added run dependent corrections on T for calibration, removed old method
[u/mrichter/AliRoot.git] / PWGGA / EMCALTasks / AliAnalysisTaskEMCALClusterize.cxx
CommitLineData
6fdd30c4 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 **************************************************************************/
6fdd30c4 15
16//_________________________________________________________________________
17// This analysis provides a new list of clusters to be used in other analysis
18//
19// Author: Gustavo Conesa Balbastre,
20// Adapted from analysis class from Deepa Thomas
21//
cd231d42 22// $Id$
6fdd30c4 23//_________________________________________________________________________
24
25// --- Root ---
04fcfa08 26#include <TString.h>
27#include <TRefArray.h>
28#include <TClonesArray.h>
29#include <TTree.h>
30#include <TGeoManager.h>
31#include <TROOT.h>
32#include <TInterpreter.h>
33#include <TFile.h>
6fdd30c4 34
35// --- AliRoot Analysis Steering
36#include "AliAnalysisTask.h"
37#include "AliAnalysisManager.h"
38#include "AliESDEvent.h"
39#include "AliGeomManager.h"
40#include "AliVCaloCells.h"
41#include "AliAODCaloCluster.h"
42#include "AliCDBManager.h"
43#include "AliCDBStorage.h"
44#include "AliCDBEntry.h"
45#include "AliLog.h"
46#include "AliVEventHandler.h"
6060ed91 47#include "AliAODInputHandler.h"
3769e0cb 48#include "AliOADBContainer.h"
6fdd30c4 49
50// --- EMCAL
6fdd30c4 51#include "AliEMCALAfterBurnerUF.h"
52#include "AliEMCALGeometry.h"
53#include "AliEMCALClusterizerNxN.h"
54#include "AliEMCALClusterizerv1.h"
cfcbe5d2 55#include "AliEMCALClusterizerv2.h"
6fdd30c4 56#include "AliEMCALRecPoint.h"
57#include "AliEMCALDigit.h"
58#include "AliCaloCalibPedestal.h"
59#include "AliEMCALCalibData.h"
60
61#include "AliAnalysisTaskEMCALClusterize.h"
62
63ClassImp(AliAnalysisTaskEMCALClusterize)
64
6544055e 65//______________________________________________________________________________
6fdd30c4 66AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
6544055e 67: AliAnalysisTaskSE(name)
68, fEvent(0)
3769e0cb 69, fGeom(0), fGeomName("")
6544055e 70, fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
71, fCalibData(0), fPedestalData(0)
3769e0cb 72, fOCDBpath(""), fAccessOCDB(kFALSE)
6544055e 73, fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
74, fRecParam(0), fClusterizer(0)
75, fUnfolder(0), fJustUnfold(kFALSE)
6a797202 76, fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
77, fFillAODFile(kFALSE), fFillAODHeader(0)
6544055e 78, fFillAODCaloCells(0), fRun(-1)
79, fRecoUtils(0), fConfigName("")
80, fCellLabels(), fCellSecondLabels(), fCellTime()
adad4ea9 81, fCellMatchdEta(), fCellMatchdPhi()
3769e0cb 82, fMaxEvent(0), fDoTrackMatching(kFALSE)
83, fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
84, fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
85, fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
86, fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
6a8be5c3 87{
b43bdd14 88 // Constructor
89
e3990982 90 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
91 for(Int_t j = 0; j < 24*48*11; j++)
92 {
6a8be5c3 93 fCellLabels[j] = -1;
94 fCellSecondLabels[j] = -1;
adad4ea9 95 fCellTime[j] = 0.;
96 fCellMatchdEta[j] = -999;
97 fCellMatchdPhi[j] = -999;
6a8be5c3 98 }
a39f5b70 99
a39f5b70 100
6fdd30c4 101}
102
6544055e 103//______________________________________________________________
6fdd30c4 104AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
6544055e 105: AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
106, fEvent(0)
3769e0cb 107, fGeom(0), fGeomName("")
6544055e 108, fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
109, fCalibData(0), fPedestalData(0)
3769e0cb 110, fOCDBpath(""), fAccessOCDB(kFALSE)
6544055e 111, fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
112, fRecParam(0), fClusterizer(0)
113, fUnfolder(0), fJustUnfold(kFALSE)
6a797202 114, fOutputAODBranch(0), fOutputAODBranchName(""), fOutputAODBranchSet(0)
115, fFillAODFile(kFALSE), fFillAODHeader(0)
6544055e 116, fFillAODCaloCells(0), fRun(-1)
117, fRecoUtils(0), fConfigName("")
118, fCellLabels(), fCellSecondLabels(), fCellTime()
adad4ea9 119, fCellMatchdEta(), fCellMatchdPhi()
3769e0cb 120, fMaxEvent(0), fDoTrackMatching(kFALSE)
121, fSelectCell(kFALSE), fSelectCellMinE(0), fSelectCellMinFrac(0)
122, fRemoveLEDEvents(kTRUE), fRemoveExoticEvents(kFALSE)
9ac44255 123, fImportGeometryFromFile(kFALSE), fImportGeometryFilePath("")
3769e0cb 124, fOADBSet(kFALSE), fAccessOADB(kTRUE), fOADBFilePath("")
125
6fdd30c4 126{
127 // Constructor
b43bdd14 128
e3990982 129 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = 0;
130 for(Int_t j = 0; j < 24*48*11; j++)
131 {
37d2296c 132 fCellLabels[j] = -1;
133 fCellSecondLabels[j] = -1;
adad4ea9 134 fCellTime[j] = 0.;
135 fCellMatchdEta[j] = -999;
136 fCellMatchdPhi[j] = -999;
37d2296c 137 }
b43bdd14 138
6fdd30c4 139}
140
5994e71f 141
6544055e 142//_______________________________________________________________
6fdd30c4 143AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
144{
145 //dtor
146
b43bdd14 147 if (fDigitsArr)
148 {
6fdd30c4 149 fDigitsArr->Clear("C");
542e1379 150 delete fDigitsArr;
6fdd30c4 151 }
152
b43bdd14 153 if (fClusterArr)
154 {
6fdd30c4 155 fClusterArr->Delete();
542e1379 156 delete fClusterArr;
6fdd30c4 157 }
158
b43bdd14 159 if (fCaloClusterArr)
160 {
6fdd30c4 161 fCaloClusterArr->Delete();
542e1379 162 delete fCaloClusterArr;
6fdd30c4 163 }
6a8be5c3 164
385b7abf 165 if(fClusterizer) delete fClusterizer;
166 if(fUnfolder) delete fUnfolder;
167 if(fRecoUtils) delete fRecoUtils;
6a8be5c3 168
6fdd30c4 169}
170
3769e0cb 171//_______________________________________________
172void AliAnalysisTaskEMCALClusterize::AccessOADB()
173{
174 // Set the AODB calibration, bad channels etc. parameters at least once
175 // alignment matrices from OADB done in SetGeometryMatrices
176
177 //Set it only once
178 if(fOADBSet) return ;
179
180 Int_t runnumber = InputEvent()->GetRunNumber() ;
181 TString pass = GetPass();
182
183 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Get AODB parameters from EMCAL in %s for run %d, and <%s> \n",fOADBFilePath.Data(),runnumber,pass.Data());
184
185 Int_t nSM = fGeom->GetNumberOfSuperModules();
186
187 // Bad map
188 if(fRecoUtils->IsBadChannelsRemovalSwitchedOn())
189 {
190 AliOADBContainer *contBC=new AliOADBContainer("");
191 contBC->InitFromFile(Form("%s/EMCALBadChannels.root",fOADBFilePath.Data()),"AliEMCALBadChannels");
192
193 TObjArray *arrayBC=(TObjArray*)contBC->GetObject(runnumber);
194
195 if(arrayBC)
196 {
3769e0cb 197 fRecoUtils->SwitchOnDistToBadChannelRecalculation();
198 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Remove EMCAL bad cells \n");
199
200 for (Int_t i=0; i<nSM; ++i)
201 {
202 TH2I *hbm = fRecoUtils->GetEMCALChannelStatusMap(i);
203
204 if (hbm)
205 delete hbm;
206
17e1c80e 207 hbm=(TH2I*)arrayBC->FindObject(Form("EMCALBadChannelMap_Mod%d",i));
3769e0cb 208
209 if (!hbm)
210 {
211 AliError(Form("Can not get EMCALBadChannelMap_Mod%d",i));
212 continue;
213 }
214
215 hbm->SetDirectory(0);
216 fRecoUtils->SetEMCALChannelStatusMap(i,hbm);
217
218 } // loop
17e1c80e 219 } else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT remove EMCAL bad channels\n"); // run array
3769e0cb 220 } // Remove bad
221
222 // Energy Recalibration
223 if(fRecoUtils->IsRecalibrationOn())
224 {
225 AliOADBContainer *contRF=new AliOADBContainer("");
226
227 contRF->InitFromFile(Form("%s/EMCALRecalib.root",fOADBFilePath.Data()),"AliEMCALRecalib");
228
229 TObjArray *recal=(TObjArray*)contRF->GetObject(runnumber);
230
231 if(recal)
232 {
233 TObjArray *recalpass=(TObjArray*)recal->FindObject(pass);
234
235 if(recalpass)
236 {
237 TObjArray *recalib=(TObjArray*)recalpass->FindObject("Recalib");
238
239 if(recalib)
240 {
241 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate EMCAL \n");
242 for (Int_t i=0; i<nSM; ++i)
243 {
244 TH2F *h = fRecoUtils->GetEMCALChannelRecalibrationFactors(i);
245
246 if (h)
247 delete h;
248
249 h = (TH2F*)recalib->FindObject(Form("EMCALRecalFactors_SM%d",i));
250
251 if (!h)
252 {
253 AliError(Form("Could not load EMCALRecalFactors_SM%d",i));
254 continue;
255 }
256
257 h->SetDirectory(0);
258
259 fRecoUtils->SetEMCALChannelRecalibrationFactors(i,h);
260 } // SM loop
17e1c80e 261 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params object array \n"); // array ok
262 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for pass\n"); // array pass ok
263 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL, no params for run\n"); // run number array ok
7bf608c9 264
3769e0cb 265 } // Recalibration on
266
7bf608c9 267 // Energy Recalibration, apply on top of previous calibration factors
268 if(fRecoUtils->IsRunDepRecalibrationOn())
269 {
270 AliOADBContainer *contRFTD=new AliOADBContainer("");
271
272 contRFTD->InitFromFile(Form("%s/EMCALTemperatureCorrCalib.root",fOADBFilePath.Data()),"AliEMCALRunDepTempCalibCorrections");
273
274 TH1S *htd=(TH1S*)contRFTD->GetObject(runnumber);
275
276 if(htd)
277 {
278 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Recalibrate (Temperature) EMCAL \n");
279
280 for (Int_t ism=0; ism<nSM; ++ism)
281 {
282 for (Int_t icol=0; icol<48; ++icol)
283 {
284 for (Int_t irow=0; irow<24; ++irow)
285 {
286 Float_t factor = fRecoUtils->GetEMCALChannelRecalibrationFactor(ism,icol,irow);
287
288 Int_t absID = fGeom->GetAbsCellIdFromCellIndexes(ism, irow, icol); // original calibration factor
289 factor *= htd->GetBinContent(absID) / 10000. ; // correction dependent on T
290 //printf("\t ism %d, icol %d, irow %d,absID %d, corrA %2.3f, corrB %2.3f, corrAB %2.3f\n",ism, icol, irow, absID,
291 // GetEMCALChannelRecalibrationFactor(ism,icol,irow) , htd->GetBinContent(absID) / 10000., factor);
292 fRecoUtils->SetEMCALChannelRecalibrationFactor(ism,icol,irow,factor);
293 } // columns
294 } // rows
295 } // SM loop
296 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate EMCAL with T variations, no params TH1 \n");
297 } // Run by Run T calibration
3769e0cb 298
299 // Time Recalibration
300 if(fRecoUtils->IsTimeRecalibrationOn())
301 {
302 AliOADBContainer *contTRF=new AliOADBContainer("");
303
f46af216 304 contTRF->InitFromFile(Form("%s/EMCALTimeCalib.root",fOADBFilePath.Data()),"AliEMCALTimeCalib");
3769e0cb 305
306 TObjArray *trecal=(TObjArray*)contTRF->GetObject(runnumber);
307
308 if(trecal)
309 {
310 TObjArray *trecalpass=(TObjArray*)trecal->FindObject(pass);
311
312 if(trecalpass)
313 {
b43bdd14 314 printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Time Recalibrate EMCAL \n");
f46af216 315 for (Int_t ibc = 0; ibc < 4; ++ibc)
3769e0cb 316 {
f46af216 317 TH1F *h = fRecoUtils->GetEMCALChannelTimeRecalibrationFactors(ibc);
318
319 if (h)
320 delete h;
321
322 h = (TH1F*)trecalpass->FindObject(Form("hAllTimeAvBC%d",ibc));
323
324 if (!h)
3769e0cb 325 {
f46af216 326 AliError(Form("Could not load hAllTimeAvBC%d",ibc));
327 continue;
328 }
329
330 h->SetDirectory(0);
331
332 fRecoUtils->SetEMCALChannelTimeRecalibrationFactors(ibc,h);
333 } // bunch crossing loop
17e1c80e 334 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for pass\n"); // array pass ok
335 }else printf("AliAnalysisTaskEMCALClusterize::SetOADBParameters() - Do NOT recalibrate time EMCAL, no params for run\n"); // run number array ok
3769e0cb 336
337 } // Time recalibration on
338
339 // Parameters already set once, so do not it again
340 fOADBSet = kTRUE;
341
342}
343
6544055e 344//_________________________________________________
345Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
346{
347 //Access to OCDB stuff
1f77507b 348
6544055e 349 fEvent = InputEvent();
350 if (!fEvent)
351 {
352 Warning("AccessOCDB","Event not available!!!");
353 return kFALSE;
1f77507b 354 }
1f77507b 355
6544055e 356 if (fEvent->GetRunNumber()==fRun)
357 return kTRUE;
358 fRun = fEvent->GetRunNumber();
1f77507b 359
6544055e 360 if(DebugLevel() > 1 )
361 printf("AliAnalysisTaksEMCALClusterize::AccessODCD() - Begin");
b43bdd14 362
6544055e 363 AliCDBManager *cdb = AliCDBManager::Instance();
1f77507b 364
1f77507b 365
b43bdd14 366 if (fOCDBpath.Length())
367 {
6544055e 368 cdb->SetDefaultStorage(fOCDBpath.Data());
369 printf("AliAnalysisTaksEMCALClusterize::AccessOCDB() - Default storage %s",fOCDBpath.Data());
370 }
371
372 cdb->SetRun(fEvent->GetRunNumber());
1f77507b 373
374 //
6544055e 375 // EMCAL from RAW OCDB
376 if (fOCDBpath.Contains("alien:"))
377 {
378 cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
379 cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
380 }
1f77507b 381
6544055e 382 TString path = cdb->GetDefaultStorage()->GetBaseFolder();
1f77507b 383
6544055e 384 // init parameters:
1f77507b 385
6544055e 386 //Get calibration parameters
387 if(!fCalibData)
388 {
389 AliCDBEntry *entry = (AliCDBEntry*)
390 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
391
392 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
1f77507b 393 }
6544055e 394
395 if(!fCalibData)
396 AliFatal("Calibration parameters not found in CDB!");
397
398 //Get calibration parameters
399 if(!fPedestalData)
400 {
401 AliCDBEntry *entry = (AliCDBEntry*)
402 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
403
404 if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
1f77507b 405 }
406
6544055e 407 if(!fPedestalData)
408 AliFatal("Dead map not found in CDB!");
1f77507b 409
6544055e 410 return kTRUE;
6fdd30c4 411}
412
9ac44255 413//_____________________________________________________
6544055e 414void AliAnalysisTaskEMCALClusterize::CheckAndGetEvent()
6fdd30c4 415{
6544055e 416 // Get the input event, it can depend in embedded events what you want to get
417 // Also check if the quality of the event is good if not reject it
418
419 fEvent = 0x0;
6a8be5c3 420
220447cd 421 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
422 Int_t eventN = Entry();
423 if(aodIH) eventN = aodIH->GetReadEntry();
424
6544055e 425 if (eventN > fMaxEvent)
426 return ;
f5036bcb 427
c50e6949 428 //printf("Clusterizer --- Event %d-- \n",eventN);
b43bdd14 429
6060ed91 430 //Check if input event are embedded events
431 //If so, take output event
3769e0cb 432 if (aodIH && aodIH->GetMergeEvents())
433 {
6544055e 434 fEvent = AODEvent();
1f77507b 435
436 if(!aodIH->GetMergeEMCALCells())
437 AliFatal("Events merged but not EMCAL cells, check analysis settings!");
438
3769e0cb 439 if(DebugLevel() > 1)
440 {
441 printf("AliAnalysisTaksEMCALClusterize::UserExec() - Use embedded events\n");
442
443 printf("\t InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
444 InputEvent()->GetEMCALCells()->GetNumberOfCells());
445
446 printf("\t MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
447 aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
448
449 for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++)
450 {
451 AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
452 if(sigCluster->IsEMCAL()) printf("\t \t Signal cluster: i %d, E %f\n",icl,sigCluster->E());
453 }
454
455 printf("\t OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
456 AODEvent()->GetEMCALCells()->GetNumberOfCells());
6544055e 457 }
6060ed91 458 }
3769e0cb 459 else
460 {
6544055e 461 fEvent = InputEvent();
1f77507b 462 if(fFillAODCaloCells) FillAODCaloCells();
463 if(fFillAODHeader) FillAODHeader();
6060ed91 464 }
465
3769e0cb 466 if (!fEvent)
467 {
3fa4fb85 468 Error("UserExec","Event not available");
6544055e 469 return ;
6fdd30c4 470 }
6a8be5c3 471
cd2e4ce6 472 //-------------------------------------------------------------------------------------
6544055e 473 // Reject events if LED was firing, use only for LHC11a data
474 // Reject event if triggered by exotic cell and remove exotic cells if not triggered
cd2e4ce6 475 //-------------------------------------------------------------------------------------
6544055e 476
3769e0cb 477 if( IsLEDEvent( InputEvent()->GetRunNumber() ) ) { fEvent = 0x0 ; return ; }
6544055e 478
3769e0cb 479 if( IsExoticEvent() ) { fEvent = 0x0 ; return ; }
cd2e4ce6 480
6a797202 481 //-------------------------------------------------------------------------------------
482 // Set the cluster array in the event (output or input)
483 //-------------------------------------------------------------------------------------
484
485 if ( fFillAODFile )
486 {
487 //Magic line to write events to AOD filem put after event rejection
488 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kTRUE);
489 }
490 else if( !fOutputAODBranchSet )
491 {
492 // Create array and put it in the input event, if output AOD not selected, only once
493 InputEvent()->AddObject(fOutputAODBranch);
494 fOutputAODBranchSet = kTRUE;
495 printf("AliAnalysisTaskEMCALClusterize::UserExec() - Add AOD branch <%s> to input event\n",fOutputAODBranchName.Data());
496 }
6544055e 497
498}
499
500//____________________________________________________
501void AliAnalysisTaskEMCALClusterize::ClusterizeCells()
502{
503 // Recluster calocells, transform them into digits,
504 // feed the clusterizer with them and get new list of clusters
505
506 //In case of MC, first loop on the clusters and fill MC label to array
507 Int_t nClusters = fEvent->GetNumberOfCaloClusters();
508 Int_t nClustersOrg = 0;
b43bdd14 509
6544055e 510 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
511 if(aodIH && aodIH->GetEventToMerge()) //Embedding
512 nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
f5036bcb 513
6544055e 514 for (Int_t i = 0; i < nClusters; i++)
515 {
516 AliVCluster *clus = 0;
517 if(aodIH && aodIH->GetEventToMerge()) //Embedding
518 clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
519 else
520 clus = fEvent->GetCaloCluster(i);
5994e71f 521
6544055e 522 if(!clus) return;
6a8be5c3 523
f46af216 524 if(clus->IsEMCAL())
525 {
6544055e 526 Int_t label = clus->GetLabel();
527 Int_t label2 = -1 ;
73a4ed8f 528 //printf("Org cluster E %f, Time %e, Id = ", clus->E(), clus->GetTOF() );
6544055e 529 if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
530 UShort_t * index = clus->GetCellsAbsId() ;
f46af216 531 for(Int_t icell=0; icell < clus->GetNCells(); icell++ )
532 {
6544055e 533 fCellLabels[index[icell]] = label;
534 fCellSecondLabels[index[icell]] = label2;
73a4ed8f 535 fCellTime[index[icell]] = clus->GetTOF();
536 fCellMatchdEta[index[icell]] = clus->GetTrackDz();
537 fCellMatchdPhi[index[icell]] = clus->GetTrackDx();
538 //printf(" %d,", index[icell] );
6fdd30c4 539 }
6544055e 540 nClustersOrg++;
6fdd30c4 541 }
b43bdd14 542 // printf("\n");
6544055e 543 }
544
545 // Transform CaloCells into Digits
b43bdd14 546
6544055e 547 Int_t idigit = 0;
548 Int_t id = -1;
549 Float_t amp = -1;
550 Double_t time = -1;
551
552 AliVCaloCells *cells = fEvent->GetEMCALCells();
b43bdd14 553
a7e5a381 554 Int_t bc = InputEvent()->GetBunchCrossNumber();
6544055e 555 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
556 {
6544055e 557 // Get cell values, recalibrate and not include bad channels found in analysis, nor cells with too low energy, nor exotic cell
558 id = cells->GetCellNumber(icell);
a7e5a381 559 Bool_t accept = fRecoUtils->AcceptCalibrateCell(id,bc,amp,time,cells);
560
561 // Do not include cells with too low energy, nor exotic cell
562 if(amp < fRecParam->GetMinECut() ) accept = kFALSE;
563
564 // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
3769e0cb 565 if (time*1e9 < 1.)
566 {
a7e5a381 567 time = fCellTime[id];
f46af216 568 //printf("cell %d time org %f - ",id, time*1.e9);
a7e5a381 569 fRecoUtils->RecalibrateCellTime(id,bc,time);
f46af216 570 //printf("recal %f\n",time*1.e9);
a7e5a381 571 }
572
3769e0cb 573 if( accept && fRecoUtils->IsExoticCell(id,cells,bc))
574 {
a7e5a381 575 accept = kFALSE;
576 }
b43bdd14 577
f46af216 578 if( !accept )
579 {
6544055e 580 fCellLabels[id] =-1; //reset the entry in the array for next event
581 fCellSecondLabels[id]=-1; //reset the entry in the array for next event
adad4ea9 582 fCellTime[id] = 0.;
583 fCellMatchdEta[id] =-999;
584 fCellMatchdPhi[id] =-999;
f46af216 585 if( DebugLevel() > 2 )
586 printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - Remove channel absId %d, index %d of %d, amp %f, time %f\n",
b43bdd14 587 id,icell, cells->GetNumberOfCells(), amp, time*1.e9);
6544055e 588 continue;
589 }
b43bdd14 590
6544055e 591 //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
592 new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
b43bdd14 593
04fcfa08 594 fCellLabels[id] =-1; //reset the entry in the array for next event
6fdd30c4 595
6544055e 596 idigit++;
597 }
6a8be5c3 598
04fcfa08 599 fDigitsArr->Sort();
b43bdd14 600
6544055e 601 //-------------------------------------------------------------------------------------
602 //Do the clusterization
603 //-------------------------------------------------------------------------------------
b43bdd14 604
6544055e 605 fClusterizer->Digits2Clusters("");
606
607 //-------------------------------------------------------------------------------------
608 //Transform the recpoints into AliVClusters
609 //-------------------------------------------------------------------------------------
610
b43bdd14 611 RecPoints2Clusters();
6544055e 612
3769e0cb 613 if(!fCaloClusterArr)
614 {
6544055e 615 printf("AliAnalysisTaksEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
616 return;
617 }
618
3769e0cb 619 if( DebugLevel() > 0 )
620 {
6544055e 621 printf("AliAnalysisTaksEMCALClusterize::ClusterizeCells() - N clusters: before recluster %d, after recluster %d\n",nClustersOrg, fCaloClusterArr->GetEntriesFast());
b43bdd14 622
623 if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast())
624 {
6544055e 625 printf("\t Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
626 fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
627 fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
628 }
629 }
630
631 //Reset the array with second labels for this event
632 memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
b43bdd14 633
6544055e 634}
a39f5b70 635
73a4ed8f 636//_____________________________________________________
6544055e 637void AliAnalysisTaskEMCALClusterize::ClusterUnfolding()
638{
639 // Take the event clusters and unfold them
640
641 AliVCaloCells *cells = fEvent->GetEMCALCells();
642 Double_t cellAmplitude = 0;
643 Double_t cellTime = 0;
644 Short_t cellNumber = 0;
77e93dc2 645 Short_t cellMCLabel = 0;
646 Double_t cellEFrac = 0;
6544055e 647 Int_t nClustersOrg = 0;
648
649 // Fill the array with the EMCAL clusters, copy them
650 for (Int_t i = 0; i < fEvent->GetNumberOfCaloClusters(); i++)
651 {
652 AliVCluster *clus = fEvent->GetCaloCluster(i);
b43bdd14 653 if(clus->IsEMCAL())
654 {
6544055e 655 //recalibrate/remove bad channels/etc if requested
b43bdd14 656 if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells()))
657 {
5994e71f 658 continue;
6544055e 659 }
3bfc4732 660
b43bdd14 661 if(fRecoUtils->IsRecalibrationOn())
662 {
6544055e 663 //Calibrate cluster
664 fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
665
666 //CalibrateCells
667 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
668 {
77e93dc2 669 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime, cellMCLabel, cellEFrac) != kTRUE)
6544055e 670 break;
671
672 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
673 fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
674 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
675
676 //Do not include bad channels found in analysis?
677 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
b43bdd14 678 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
679 {
6544055e 680 fCellLabels[cellNumber] =-1; //reset the entry in the array for next event
681 fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
adad4ea9 682 fCellTime[cellNumber] = 0.;
683 fCellMatchdEta[cellNumber] =-999;
684 fCellMatchdPhi[cellNumber] =-999;
6544055e 685 continue;
686 }
687
688 cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
689
690 }// cells loop
691 }// recalibrate
5994e71f 692
6544055e 693 //Cast to ESD or AOD, needed to create the cluster array
694 AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
695 AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
b43bdd14 696
697 if (esdCluster)
698 {
6544055e 699 fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
700 }//ESD
b43bdd14 701 else if(aodCluster)
702 {
6544055e 703 fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
704 }//AOD
705 else
706 Warning("UserExec()"," - Wrong CaloCluster type?");
b43bdd14 707
6544055e 708 nClustersOrg++;
6a8be5c3 709 }
6fdd30c4 710 }
711
6544055e 712 //Do the unfolding
713 fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
6a8be5c3 714
6544055e 715 //CLEAN-UP
716 fUnfolder->Clear();
385b7abf 717
6544055e 718}
719
720//_____________________________________________________
721void AliAnalysisTaskEMCALClusterize::FillAODCaloCells()
722{
723 // Put calo cells in standard branch
724 AliVCaloCells &eventEMcells = *(fEvent->GetEMCALCells());
725 Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
6fdd30c4 726
6544055e 727 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
728 aodEMcells.CreateContainer(nEMcell);
729 aodEMcells.SetType(AliVCaloCells::kEMCALCell);
730 Double_t calibFactor = 1.;
b43bdd14 731 for (Int_t iCell = 0; iCell < nEMcell; iCell++)
732 {
6544055e 733 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
734 fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
735 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
385b7abf 736
b43bdd14 737 if(fRecoUtils->IsRecalibrationOn())
738 {
6544055e 739 calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
385b7abf 740 }
0821ba9f 741
b43bdd14 742 if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi))
743 { //Channel is not declared as bad
77e93dc2 744 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor,
745 eventEMcells.GetTime(iCell),eventEMcells.GetMCLabel(iCell),eventEMcells.GetEFraction(iCell));
6544055e 746 }
b43bdd14 747 else
748 {
77e93dc2 749 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0,-1,-1,0);
5994e71f 750 }
6fdd30c4 751 }
6544055e 752 aodEMcells.Sort();
6fdd30c4 753
6544055e 754}
6fdd30c4 755
6544055e 756//__________________________________________________
757void AliAnalysisTaskEMCALClusterize::FillAODHeader()
6fdd30c4 758{
6544055e 759 //Put event header information in standard AOD branch
e9dd2d80 760
6544055e 761 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (fEvent);
762 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (fEvent);
e9dd2d80 763
6544055e 764 Double_t pos[3] ;
765 Double_t covVtx[6];
766 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
767
768 AliAODHeader* header = AODEvent()->GetHeader();
769 header->SetRunNumber(fEvent->GetRunNumber());
770
b43bdd14 771 if(esdevent)
772 {
6544055e 773 TTree* tree = fInputHandler->GetTree();
b43bdd14 774 if (tree)
775 {
6544055e 776 TFile* file = tree->GetCurrentFile();
777 if (file) header->SetESDFileName(file->GetName());
778 }
779 }
780 else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
e9dd2d80 781
6544055e 782 header->SetBunchCrossNumber(fEvent->GetBunchCrossNumber());
783 header->SetOrbitNumber(fEvent->GetOrbitNumber());
784 header->SetPeriodNumber(fEvent->GetPeriodNumber());
785 header->SetEventType(fEvent->GetEventType());
385b7abf 786
6544055e 787 //Centrality
b43bdd14 788 if(fEvent->GetCentrality())
6544055e 789 header->SetCentrality(new AliCentrality(*(fEvent->GetCentrality())));
b43bdd14 790 else
6544055e 791 header->SetCentrality(0);
385b7abf 792
6544055e 793 //Trigger
794 header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
795 if (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
796 else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
b43bdd14 797
6544055e 798 header->SetTriggerMask(fEvent->GetTriggerMask());
799 header->SetTriggerCluster(fEvent->GetTriggerCluster());
b43bdd14 800
801 if (esdevent)
802 {
6544055e 803 header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
804 header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
805 header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
806 }
b43bdd14 807 else if (aodevent)
808 {
6544055e 809 header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
810 header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
811 header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
812 }
e9dd2d80 813
6544055e 814 header->SetMagneticField(fEvent->GetMagneticField());
815 //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
e9dd2d80 816
6544055e 817 header->SetZDCN1Energy(fEvent->GetZDCN1Energy());
818 header->SetZDCP1Energy(fEvent->GetZDCP1Energy());
819 header->SetZDCN2Energy(fEvent->GetZDCN2Energy());
820 header->SetZDCP2Energy(fEvent->GetZDCP2Energy());
821 header->SetZDCEMEnergy(fEvent->GetZDCEMEnergy(0),fEvent->GetZDCEMEnergy(1));
385b7abf 822
6544055e 823 Float_t diamxy[2]={fEvent->GetDiamondX(),fEvent->GetDiamondY()};
824 Float_t diamcov[3];
825 fEvent->GetDiamondCovXY(diamcov);
826 header->SetDiamond(diamxy,diamcov);
827 if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
828 else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
e9dd2d80 829
6544055e 830 //
831 Int_t nVertices = 1 ;/* = prim. vtx*/;
832 Int_t nCaloClus = fEvent->GetNumberOfCaloClusters();
e9dd2d80 833
6544055e 834 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
385b7abf 835
6544055e 836 // Access to the AOD container of vertices
837 TClonesArray &vertices = *(AODEvent()->GetVertices());
838 Int_t jVertices=0;
385b7abf 839
6544055e 840 // Add primary vertex. The primary tracks will be defined
841 // after the loops on the composite objects (V0, cascades, kinks)
842 fEvent->GetPrimaryVertex()->GetXYZ(pos);
843 Float_t chi = 0;
b43bdd14 844 if (esdevent)
845 {
6544055e 846 esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
847 chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
848 }
b43bdd14 849 else if (aodevent)
850 {
6544055e 851 aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
852 chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
6fdd30c4 853 }
bd9e8ebd 854
6544055e 855 AliAODVertex * primary = new(vertices[jVertices++])
856 AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
857 primary->SetName(fEvent->GetPrimaryVertex()->GetName());
858 primary->SetTitle(fEvent->GetPrimaryVertex()->GetTitle());
e9dd2d80 859
6544055e 860}
861
b43bdd14 862//_________________________________________________________
863void AliAnalysisTaskEMCALClusterize::FillCaloClusterInAOD()
864{
865 // Get the CaloClusters array, do some final calculations
866 // and put the clusters in the AOD branch
867
868 //First recalculate track-matching for the new clusters
869 if(fDoTrackMatching)
870 {
871 fRecoUtils->FindMatches(fEvent,fCaloClusterArr,fGeom);
872 }
873 //Put the new clusters in the AOD list
874
875 Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
876 for(Int_t i = 0; i < kNumberOfCaloClusters; i++)
877 {
878 AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
879
880 newCluster->SetID(i);
881
882 //Add matched track
883 if(fDoTrackMatching)
884 {
885 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
886 if(trackIndex >= 0)
887 {
888 newCluster->AddTrackMatched(fEvent->GetTrack(trackIndex));
889 if(DebugLevel() > 1)
890 printf("AliAnalysisTaksEMCALClusterize::UserExec() - Matched Track index %d to new cluster %d \n",trackIndex,i);
891 }
892
893 Float_t dR = 999., dZ = 999.;
894 fRecoUtils->GetMatchedResiduals(newCluster->GetID(),dR,dZ);
895 newCluster->SetTrackDistance(dR,dZ);
896
897 }
898 else
899 {// Assign previously assigned matched track in reco, very very rough
900 Int_t absId0 = newCluster->GetCellsAbsId()[0]; // Assign match of first cell in cluster
901 newCluster->SetTrackDistance(fCellMatchdPhi[absId0],fCellMatchdEta[absId0]);
902 }
903
904 //printf("New cluster E %f, Time %e, Id = ", newCluster->E(), newCluster->GetTOF() );
905 //for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ) printf(" %d,", newCluster->GetCellsAbsId() [icell] );
906 //printf("\n");
907
908 // Calculate distance to bad channel for new cluster. Make sure you give the list of bad channels.
909 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, fEvent->GetEMCALCells(), newCluster);
910
911 new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
912
913 if(DebugLevel() > 1 )
914 printf("AliAnalysisTaksEMCALClusterize::UserExec() - New cluster %d of %d, energy %f\n",newCluster->GetID(), kNumberOfCaloClusters, newCluster->E());
915
916 } // cluster loop
917
918 fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
919
920 // Clean up
921 fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
922
923}
924
3769e0cb 925//_______________________________________________
926TString AliAnalysisTaskEMCALClusterize::GetPass()
927{
928 // Get passx from filename.
929
930 if (!AliAnalysisManager::GetAnalysisManager()->GetTree())
931 {
932 AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Pointer to tree = 0, returning null\n");
933 return TString("");
934 }
935
936 if (!AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile())
937 {
938 AliError("AliAnalysisTaskEMCALClusterize::GetPass() - Null pointer input file, returning null\n");
939 return TString("");
940 }
941
942 TString pass(AliAnalysisManager::GetAnalysisManager()->GetTree()->GetCurrentFile()->GetName());
943 if (pass.Contains("ass1")) return TString("pass1");
944 else if (pass.Contains("ass2")) return TString("pass2");
945 else if (pass.Contains("ass3")) return TString("pass3");
946
947 // No condition fullfilled, give a default value
948 printf("AliAnalysisTaskEMCALClusterize::GetPass() - Pass number string not found \n");
949 return TString("");
950
951}
952
6544055e 953//_________________________________________
954void AliAnalysisTaskEMCALClusterize::Init()
955{
956 //Init analysis with configuration macro if available
957
3769e0cb 958 fOADBSet = kFALSE;
959 if(fOADBFilePath == "") fOADBFilePath = "$ALICE_ROOT/OADB/EMCAL" ;
960
3769e0cb 961 fCaloClusterArr = new TObjArray(10000);
962 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
963
964 if(!fRecParam) fRecParam = new AliEMCALRecParam;
965 if(!fRecoUtils) fRecoUtils = new AliEMCALRecoUtils();
966
967 if(fMaxEvent <= 0) fMaxEvent = 1000000000;
968 if(fSelectCellMinE <= 0) fSelectCellMinE = 0.005;
969 if(fSelectCellMinFrac <= 0) fSelectCellMinFrac = 0.001;
970
971 if (fOCDBpath == "") fOCDBpath = "raw://" ;
972 if (fOutputAODBranchName == "") fOutputAODBranchName = "newEMCALClusters" ;
973
974 if(gROOT->LoadMacro(fConfigName) >=0)
975 {
6544055e 976 printf("AliAnalysisTaksEMCALClusterize::Init() - Configure analysis with %s\n",fConfigName.Data());
977 AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
978 fGeomName = clus->fGeomName;
979 fLoadGeomMatrices = clus->fLoadGeomMatrices;
980 fOCDBpath = clus->fOCDBpath;
981 fAccessOCDB = clus->fAccessOCDB;
982 fRecParam = clus->fRecParam;
983 fJustUnfold = clus->fJustUnfold;
984 fFillAODFile = clus->fFillAODFile;
985 fRecoUtils = clus->fRecoUtils;
986 fConfigName = clus->fConfigName;
987 fMaxEvent = clus->fMaxEvent;
988 fDoTrackMatching = clus->fDoTrackMatching;
989 fOutputAODBranchName = clus->fOutputAODBranchName;
e3990982 990 for(Int_t i = 0; i < 12; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
6fdd30c4 991 }
e9dd2d80 992
9ac44255 993 // Init geometry, I do not like much to do it like this ...
3769e0cb 994 if(fImportGeometryFromFile && !gGeoManager)
995 {
996 if (fImportGeometryFilePath == "") fImportGeometryFilePath = "$ALICE_ROOT/PWGGA/EMCALTasks/macros/geometry.root" ; // "$ALICE_ROOT/EVE/alice-data/default_geo.root"
997 printf("AliAnalysisTaskEMCALClusterize::Init() - Import %s\n",fImportGeometryFilePath.Data());
998 TGeoManager::Import(fImportGeometryFilePath) ;
9ac44255 999 }
b43bdd14 1000
6544055e 1001}
6fdd30c4 1002
6544055e 1003//_______________________________________________________
6fdd30c4 1004void AliAnalysisTaskEMCALClusterize::InitClusterization()
1005{
1006 //Select clusterization/unfolding algorithm and set all the needed parameters
1007
04fcfa08 1008 if (fJustUnfold)
1009 {
bd9e8ebd 1010 // init the unfolding afterburner
1011 delete fUnfolder;
1759f743 1012 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut(),fRecParam->GetMinECut());
bd9e8ebd 1013 return;
6a8be5c3 1014 }
1015
bd9e8ebd 1016 //First init the clusterizer
1017 delete fClusterizer;
1018 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
1019 fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
cfcbe5d2 1020 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
1021 fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
b43bdd14 1022 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
1023 {
bd9e8ebd 1024 fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
212ac797 1025 fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
1026 fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
04fcfa08 1027 }
1028 else
1029 {
bd9e8ebd 1030 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
1031 }
6a8be5c3 1032
bd9e8ebd 1033 //Now set the parameters
1034 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
1035 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
1036 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
1037 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
1038 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
1039 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
1040 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
1041 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
1042 fClusterizer->SetInputCalibrated ( kTRUE );
2af35150 1043 fClusterizer->SetJustClusters ( kTRUE );
04fcfa08 1044
1045 // Initialize the cluster rec points and digits arrays and get them.
1046 fClusterizer->SetOutput(0);
1047 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
1048 fDigitsArr = fClusterizer->GetDigits();
1049
bd9e8ebd 1050 //In case of unfolding after clusterization is requested, set the corresponding parameters
04fcfa08 1051 if(fRecParam->GetUnfold())
1052 {
bd9e8ebd 1053 Int_t i=0;
04fcfa08 1054 for (i = 0; i < 8; i++)
1055 {
bd9e8ebd 1056 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
1057 }//end of loop over parameters
04fcfa08 1058
1059 for (i = 0; i < 3; i++)
1060 {
bd9e8ebd 1061 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
1062 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
1063 }//end of loop over parameters
6fdd30c4 1064
bd9e8ebd 1065 fClusterizer->InitClusterUnfolding();
6544055e 1066
bd9e8ebd 1067 }// to unfold
6fdd30c4 1068}
3fa4fb85 1069
6544055e 1070//_________________________________________________
1071void AliAnalysisTaskEMCALClusterize::InitGeometry()
1072{
3769e0cb 1073 // Init geometry and set the geometry matrix, for the first event, skip the rest
6544055e 1074 // Also set once the run dependent calibrations
1075
3769e0cb 1076 if(fGeomMatrixSet) return;
1077
1078 Int_t runnumber = InputEvent()->GetRunNumber() ;
1079 if (!fGeom)
1080 {
1081 if(fGeomName=="")
1082 {
1083 if (runnumber < 140000) fGeomName = "EMCAL_FIRSTYEARV1";
1084 else if(runnumber < 171000) fGeomName = "EMCAL_COMPLETEV1";
1085 else fGeomName = "EMCAL_COMPLETE12SMV1";
1086 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Set EMCAL geometry name to <%s> for run %d\n",fGeomName.Data(),runnumber);
1087 }
6544055e 1088
3769e0cb 1089 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
6544055e 1090
3769e0cb 1091 if(fDebug > 0)
1092 {
1093 printf("AliAnalysisTaskEMCALClusterize::InitGeometry(run=%d)",runnumber);
1094 if (!gGeoManager) printf(" - Careful!, gGeoManager not loaded, load misalign matrices");
1095 printf("\n");
1096 }
1097 } // geometry pointer did not exist before
1098
1099 if(fLoadGeomMatrices)
1100 {
1101 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - Load user defined EMCAL geometry matrices\n");
1102
1103 // OADB if available
1104 AliOADBContainer emcGeoMat("AliEMCALgeo");
1105 emcGeoMat.InitFromFile(Form("%s/EMCALlocal2master.root",fOADBFilePath.Data()),"AliEMCALgeo");
1106 TObjArray *matEMCAL=(TObjArray*)emcGeoMat.GetObject(runnumber,"EmcalMatrices");
1107
1108
1109 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1110 {
1111
1112 if (!fGeomMatrix[mod]) // Get it from OADB
1113 {
1114 if(fDebug > 1 )
1115 printf("AliAnalysisTaskEMCALClusterize::InitGeometry() - EMCAL matrices SM %d, %p\n",
1116 mod,((TGeoHMatrix*) matEMCAL->At(mod)));
1117 //((TGeoHMatrix*) matEMCAL->At(mod))->Print();
1118
1119 fGeomMatrix[mod] = (TGeoHMatrix*) matEMCAL->At(mod) ;
1120 }
1121
1122 if(fGeomMatrix[mod])
1123 {
1124 if(DebugLevel() > 1)
1125 fGeomMatrix[mod]->Print();
1126
1127 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
1128 }
1129
1130 fGeomMatrixSet=kTRUE;
1131
1132 }//SM loop
1133 }//Load matrices
1134 else if(!gGeoManager)
1135 {
1136 printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - Get geo matrices from data");
1137 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
1138 if(!strcmp(fEvent->GetName(),"AliAODEvent"))
1139 {
1140 if(DebugLevel() > 1)
1141 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
1142 }//AOD
1143 else
1144 {
1145 if(DebugLevel() > 1)
1146 printf("AliAnalysisTaksEMCALClusterize::InitGeometry() - AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
1147
1148 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(fEvent) ;
1149
1150 if(!esd)
1151 {
1152 Error("InitGeometry"," - This event does not contain ESDs?");
6a797202 1153 if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
3769e0cb 1154 return;
1155 }
1156
1157 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++)
1158 {
1159 if(DebugLevel() > 1)
1160 esd->GetEMCALMatrix(mod)->Print();
1161
1162 if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
1163
1164 }
1165
1166 fGeomMatrixSet=kTRUE;
1167
1168 }//ESD
1169 }//Load matrices from Data
1170
6544055e 1171}
1172
6544055e 1173//____________________________________________________
1174Bool_t AliAnalysisTaskEMCALClusterize::IsExoticEvent()
1175{
1176
1177 // Check if event is exotic, get an exotic cell and compare with triggered patch
1178 // If there is a match remove event ... to be completed, filled with something provisional
1179
1180 if(!fRemoveExoticEvents) return kFALSE;
1181
1182 // Loop on cells
1183 AliVCaloCells * cells = fEvent->GetEMCALCells();
1184 Float_t totCellE = 0;
a7e5a381 1185 Int_t bc = InputEvent()->GetBunchCrossNumber();
b43bdd14 1186 for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1187 {
6544055e 1188
1189 Float_t ecell = 0 ;
1190 Double_t tcell = 0 ;
1191
1192 Int_t absID = cells->GetCellNumber(icell);
a7e5a381 1193 Bool_t accept = fRecoUtils->AcceptCalibrateCell(absID,bc,ecell,tcell,cells);
1194 if(accept && !fRecoUtils->IsExoticCell(absID,cells,bc)) totCellE += ecell;
6544055e 1195 }
1196
1197 // TString triggerclasses = "";
1198 // if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses();
1199 // else triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
1200 // //
1201 // printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
6a797202 1202 // if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
6544055e 1203 // return;
1204 //
1205
1206 //printf("TotE cell %f\n",totCellE);
1207 if(totCellE < 1) return kTRUE;
1208
1209 return kFALSE;
1210
1211}
1212
3769e0cb 1213//________________________________________________________________
1214Bool_t AliAnalysisTaskEMCALClusterize::IsLEDEvent(const Int_t run)
6544055e 1215{
1216 //Check if event is LED
1217
1218 if(!fRemoveLEDEvents) return kFALSE;
3769e0cb 1219
1220 //check events of LHC11a period
1221 if(run < 140000 || run > 146860) return kFALSE ;
1222
6544055e 1223 // Count number of cells with energy larger than 0.1 in SM3, cut on this number
1224 Int_t ncellsSM3 = 0;
1225 Int_t ncellsSM4 = 0;
1226 AliVCaloCells * cells = fEvent->GetEMCALCells();
3769e0cb 1227 for(Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
1228 {
6544055e 1229 if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
1230 if(cells->GetAmplitude(icell) > 0.1 && cells->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
1231 }
1232
1233 TString triggerclasses = "";
1234
1235 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(fEvent);
1236 if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses();
1237 else triggerclasses = ((AliAODEvent*)fEvent)->GetFiredTriggerClasses();
1238
1239 Int_t ncellcut = 21;
1240 if(triggerclasses.Contains("EMC")) ncellcut = 35;
1241
3769e0cb 1242 if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 )
1243 {
6544055e 1244 printf("AliAnalysisTaksEMCALClusterize::IsLEDEvent() - reject event %d with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
6a797202 1245 if(fFillAODFile) AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);;
6544055e 1246 return kTRUE;
1247 }
1248
1249 return kFALSE;
1250
1251}
1252
b43bdd14 1253//_______________________________________________________
1254void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters()
6fdd30c4 1255{
1256 // Restore clusters from recPoints
b43bdd14 1257 // Cluster energy, global position, cells and their amplitude
1258 // fractions are restored
1259
6a8be5c3 1260 Int_t j = 0;
b43bdd14 1261 for(Int_t i = 0; i < fClusterArr->GetEntriesFast(); i++)
6fdd30c4 1262 {
b43bdd14 1263 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) fClusterArr->At(i);
6fdd30c4 1264
1265 const Int_t ncells = recPoint->GetMultiplicity();
f6e45fe0 1266 Int_t ncellsTrue = 0;
6fdd30c4 1267
6a8be5c3 1268 if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
1269
6fdd30c4 1270 // cells and their amplitude fractions
1271 UShort_t absIds[ncells];
1272 Double32_t ratios[ncells];
1273
f6e45fe0 1274 //For later check embedding
1275 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
f053064e 1276
1277 Float_t clusterE = 0;
f46af216 1278 for (Int_t c = 0; c < ncells; c++)
1279 {
b43bdd14 1280 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->At(recPoint->GetDigitsList()[c]);
6fdd30c4 1281
f6e45fe0 1282 absIds[ncellsTrue] = digit->GetId();
1283 ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
6544055e 1284
d390e7cb 1285 // In case of unfolding, remove digits with unfolded energy too low
b43bdd14 1286 if(fSelectCell)
1287 {
f46af216 1288 if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac)
1289 {
1290 if(DebugLevel() > 1)
1291 {
6544055e 1292 printf("AliAnalysisTaksEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
d390e7cb 1293 recPoint->GetEnergiesList()[c],digit->GetAmplitude());
1294 }
1295
1296 continue;
1297
1298 } // if cuts
1299 }// Select cells
6544055e 1300
d390e7cb 1301 //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
1302 ncellsTrue++;
1303 clusterE +=recPoint->GetEnergiesList()[c];
f053064e 1304
f6e45fe0 1305 // In case of embedding, fill ratio with amount of signal,
f46af216 1306 if (aodIH && aodIH->GetMergeEvents())
1307 {
f6e45fe0 1308 //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
1309 AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
1310 AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
1311
1312 Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1313 //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1314 Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
1315 //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
1316
1317 if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
1318 //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
1319
1320 }//Embedding
1321
f6e45fe0 1322 }// cluster cell loop
6fdd30c4 1323
f46af216 1324 if (ncellsTrue < 1)
1325 {
f053064e 1326 if (DebugLevel() > 1)
1327 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
1328 recPoint->GetEnergy(), ncells);
1329 continue;
1330 }
1331
1332 //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
1333
f46af216 1334 if(clusterE < fRecParam->GetClusteringThreshold())
1335 {
f053064e 1336 if (DebugLevel()>1)
1337 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
6fdd30c4 1338 continue;
1339 }
1340
1341 TVector3 gpos;
1342 Float_t g[3];
1343
1344 // calculate new cluster position
b43bdd14 1345 recPoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
6fdd30c4 1346 recPoint->GetGlobalPosition(gpos);
1347 gpos.GetXYZ(g);
1348
1349 // create a new cluster
b43bdd14 1350 (*fCaloClusterArr)[j] = new AliAODCaloCluster() ;
1351 AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( fCaloClusterArr->At(j) ) ;
6a8be5c3 1352 j++;
6fdd30c4 1353 clus->SetType(AliVCluster::kEMCALClusterv1);
6544055e 1354 clus->SetE(clusterE);
6fdd30c4 1355 clus->SetPosition(g);
f6e45fe0 1356 clus->SetNCells(ncellsTrue);
6fdd30c4 1357 clus->SetCellsAbsId(absIds);
1358 clus->SetCellsAmplitudeFraction(ratios);
6fdd30c4 1359 clus->SetChi2(-1); //not yet implemented
bd9e8ebd 1360 clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
6fdd30c4 1361 clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
5994e71f 1362 clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
b43bdd14 1363
f46af216 1364 if(ncells == ncellsTrue)
1365 {
f053064e 1366 Float_t elipAxis[2];
1367 recPoint->GetElipsAxis(elipAxis);
1368 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1369 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1370 clus->SetDispersion(recPoint->GetDispersion());
1371 }
f46af216 1372 else if(fSelectCell)
1373 {
afaaef51 1374 // In case some cells rejected, in unfolding case, recalculate
1375 // shower shape parameters and position
d20bc070 1376 if(DebugLevel() > 1)
f46af216 1377 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Cells removed from cluster (ncells %d, ncellsTrue %d), recalculate Shower Shape\n",ncells,ncellsTrue);
b43bdd14 1378
f053064e 1379 AliVCaloCells* cells = 0x0;
1380 if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1381 else cells = InputEvent()->GetEMCALCells();
b43bdd14 1382
f053064e 1383 fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1384 fRecoUtils->RecalculateClusterPID(clus);
afaaef51 1385 fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
f053064e 1386
f053064e 1387 }
6060ed91 1388
1389 //MC
b65291bf 1390 Int_t parentMult = 0;
6060ed91 1391 Int_t *parentList = recPoint->GetParents(parentMult);
1392 clus->SetLabel(parentList, parentMult);
1393
37d2296c 1394 //Write the second major contributor to each MC cluster.
1395 Int_t iNewLabel ;
f46af216 1396 for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ )
1397 {
37d2296c 1398
1399 Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
b43bdd14 1400 if(idCell>=0)
1401 {
6ffa213f 1402 iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1403 for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1404 {
1405 if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1406 if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1407 }
1408 if (iNewLabel == 1)
1409 {
937b6c91 1410 Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
f46af216 1411 for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ )
1412 {
6ffa213f 1413 newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1414 }
6a8be5c3 1415
6ffa213f 1416 newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1417 clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
abf0860a 1418 delete [] newLabelArray;
6ffa213f 1419 }
1420 }//positive cell number
37d2296c 1421 }
1422
6fdd30c4 1423 } // recPoints loop
6544055e 1424
1425}
1426
3769e0cb 1427
6544055e 1428//____________________________________________________________
1429void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
1430{
1431 // Init geometry, create list of output clusters
1432
6a797202 1433
1434 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1435
1436 if(fOutputAODBranchName.Length()==0)
3769e0cb 1437 {
6a797202 1438 fOutputAODBranchName = "newEMCALClustersArray";
1439 printf("Cluster branch name not set, set it to newEMCALClustersArray \n");
6544055e 1440 }
6a797202 1441
1442 fOutputAODBranch->SetName(fOutputAODBranchName);
1443
1444 if( fFillAODFile )
3769e0cb 1445 {
6a797202 1446 //fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
1447
1448
1449 //fOutputAODBranch->SetOwner(kFALSE);
1450
1451 AddAODBranch("TClonesArray", &fOutputAODBranch);
6544055e 1452 }
b43bdd14 1453
6fdd30c4 1454}
6544055e 1455
1456//_______________________________________________________
1457void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
1458{
b43bdd14 1459 // Do clusterization event by event, execute different steps
1460 // 1) Do some checks on the kind of events (ESD, AOD) or if some filtering is needed and initializations
1461 // 2) Clusterize a) just unfolding existing clusters (fJustUnfold)
1462 // b) recluster cells
1463 // + convert cells into digits (calibrating them)
1464 // + recluster digits into recPoints with the chosen clusterizer (V1, V1+Unfold,V2, NxN)
1465 // with methods in AliEMCALClusterizer
1466 // + transform recPoints into CaloClusters
1467 // 3) Do some final calculations in the found clusters (track-matching) and put them in an AOD branch
1468
1469 // Step 1
1470
1471 //Remove the contents of AOD branch output list set in the previous event
6544055e 1472 fOutputAODBranch->Clear("C");
1473
1474 LoadBranches();
1475
1476 //Get the event, do some checks and settings
1477 CheckAndGetEvent() ;
1478
a1eeb25d 1479 if (!fEvent)
1480 {
6544055e 1481 if(DebugLevel() > 0 ) printf("AliAnalysisTaksEMCALClusterize::UserExec() - Skip Event %d", (Int_t) Entry());
1482 return ;
1483 }
6a797202 1484
3769e0cb 1485 //Init pointers, geometry, clusterizer, ocdb, aodb
6544055e 1486
3769e0cb 1487 InitGeometry(); // only once, must be done before OADB, geo OADB accessed here
6544055e 1488
3769e0cb 1489 if(fAccessOCDB) AccessOCDB();
1490 if(fAccessOADB) AccessOADB(); // only once
b43bdd14 1491
3769e0cb 1492 InitClusterization();
6544055e 1493
b43bdd14 1494 // Step 2
1495
6544055e 1496 // Make clusters
1497 if (fJustUnfold) ClusterUnfolding();
1498 else ClusterizeCells() ;
6544055e 1499
b43bdd14 1500 // Step 3
6544055e 1501
b43bdd14 1502 FillCaloClusterInAOD();
b65291bf 1503
6544055e 1504}
1505