memory leak fixed
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / 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 **************************************************************************/
15/* $Id: $ */
16
17//_________________________________________________________________________
18// This analysis provides a new list of clusters to be used in other analysis
19//
20// Author: Gustavo Conesa Balbastre,
21// Adapted from analysis class from Deepa Thomas
22//
23//
24//_________________________________________________________________________
25
26// --- Root ---
27#include "TString.h"
28#include "TRefArray.h"
29#include "TClonesArray.h"
30#include "TTree.h"
31#include "TGeoManager.h"
5994e71f 32#include "TROOT.h"
33#include "TInterpreter.h"
1f77507b 34#include "TFile.h"
01313069 35//#include "string.h"
6fdd30c4 36
37// --- AliRoot Analysis Steering
38#include "AliAnalysisTask.h"
39#include "AliAnalysisManager.h"
40#include "AliESDEvent.h"
41#include "AliGeomManager.h"
42#include "AliVCaloCells.h"
43#include "AliAODCaloCluster.h"
44#include "AliCDBManager.h"
45#include "AliCDBStorage.h"
46#include "AliCDBEntry.h"
47#include "AliLog.h"
48#include "AliVEventHandler.h"
6060ed91 49#include "AliAODInputHandler.h"
6fdd30c4 50
51// --- EMCAL
52#include "AliEMCALRecParam.h"
53#include "AliEMCALAfterBurnerUF.h"
54#include "AliEMCALGeometry.h"
55#include "AliEMCALClusterizerNxN.h"
56#include "AliEMCALClusterizerv1.h"
cfcbe5d2 57#include "AliEMCALClusterizerv2.h"
6fdd30c4 58#include "AliEMCALRecPoint.h"
59#include "AliEMCALDigit.h"
60#include "AliCaloCalibPedestal.h"
61#include "AliEMCALCalibData.h"
385b7abf 62#include "AliEMCALRecoUtils.h"
6fdd30c4 63
64#include "AliAnalysisTaskEMCALClusterize.h"
65
66ClassImp(AliAnalysisTaskEMCALClusterize)
67
68//________________________________________________________________________
69AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
70 : AliAnalysisTaskSE(name)
dd9a9f63 71 , fGeom(0), fGeomName("EMCAL_COMPLETEV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
e9dd2d80 72 , fCalibData(0), fPedestalData(0), fOCDBpath("raw://"), fAccessOCDB(kFALSE)
6fdd30c4 73 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
74 , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
1f77507b 75 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
76 , fFillAODFile(kTRUE), fFillAODHeader(0), fFillAODCaloCells(0)
37d2296c 77 , fRun(-1), fRecoUtils(0), fConfigName("")
cb3001b9 78 , fCellLabels(), fCellSecondLabels(), fMaxEvent(1000000000), fDoTrackMatching(kFALSE)
6fdd30c4 79
6a8be5c3 80{
6fdd30c4 81 //ctor
6060ed91 82 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
6a8be5c3 83 for(Int_t j = 0; j < 24*48*11; j++) {
84 fCellLabels[j] = -1;
85 fCellSecondLabels[j] = -1;
86 }
6fdd30c4 87 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
88 fClusterArr = new TObjArray(100);
c0feb73e 89 fCaloClusterArr = new TObjArray(1000);
6fdd30c4 90 fRecParam = new AliEMCALRecParam;
385b7abf 91 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
92 fRecoUtils = new AliEMCALRecoUtils();
6fdd30c4 93}
94
95//________________________________________________________________________
96AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
97 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
dd9a9f63 98 , fGeom(0), fGeomName("EMCAL_COMPLETEV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
e9dd2d80 99 , fCalibData(0), fPedestalData(0), fOCDBpath("raw://"), fAccessOCDB(kFALSE)
1f77507b 100 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
101 , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
102 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
103 , fFillAODFile(kFALSE), fFillAODHeader(0), fFillAODCaloCells(0)
104 , fRun(-1), fRecoUtils(0), fConfigName("")
cb3001b9 105 , fCellLabels(), fCellSecondLabels(), fMaxEvent(1000000000), fDoTrackMatching(kFALSE)
6fdd30c4 106{
107 // Constructor
6060ed91 108 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
37d2296c 109 for(Int_t j = 0; j < 24*48*11; j++) {
110 fCellLabels[j] = -1;
111 fCellSecondLabels[j] = -1;
112 }
6fdd30c4 113 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
114 fClusterArr = new TObjArray(100);
115 fCaloClusterArr = new TObjArray(100);
116 fRecParam = new AliEMCALRecParam;
385b7abf 117 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
118 fRecoUtils = new AliEMCALRecoUtils();
6fdd30c4 119}
120
5994e71f 121
6fdd30c4 122//________________________________________________________________________
123AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
124{
125 //dtor
126
127 if (fDigitsArr){
128 fDigitsArr->Clear("C");
542e1379 129 delete fDigitsArr;
6fdd30c4 130 }
131
132 if (fClusterArr){
133 fClusterArr->Delete();
542e1379 134 delete fClusterArr;
6fdd30c4 135 }
136
137 if (fCaloClusterArr){
138 fCaloClusterArr->Delete();
542e1379 139 delete fCaloClusterArr;
6fdd30c4 140 }
6a8be5c3 141
385b7abf 142 if(fClusterizer) delete fClusterizer;
143 if(fUnfolder) delete fUnfolder;
144 if(fRecoUtils) delete fRecoUtils;
6a8be5c3 145
6fdd30c4 146}
147
148//-------------------------------------------------------------------
5994e71f 149void AliAnalysisTaskEMCALClusterize::Init()
150{
151 //Init analysis with configuration macro if available
152
153 if(gROOT->LoadMacro(fConfigName) >=0){
154 printf("Configure analysis with %s\n",fConfigName.Data());
155 AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
156 fGeomName = clus->fGeomName;
157 fLoadGeomMatrices = clus->fLoadGeomMatrices;
2ef5608f 158 fOCDBpath = clus->fOCDBpath;
159 fAccessOCDB = clus->fAccessOCDB;
5994e71f 160 fRecParam = clus->fRecParam;
161 fJustUnfold = clus->fJustUnfold;
162 fFillAODFile = clus->fFillAODFile;
163 fRecoUtils = clus->fRecoUtils;
164 fConfigName = clus->fConfigName;
2ef5608f 165 fMaxEvent = clus->fMaxEvent;
166 fDoTrackMatching = clus->fDoTrackMatching;
5994e71f 167 fOutputAODBranchName = clus->fOutputAODBranchName;
168 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
6a8be5c3 169
5994e71f 170 }
6a8be5c3 171
5994e71f 172}
173
174//-------------------------------------------------------------------
6fdd30c4 175void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
176{
177 // Init geometry, create list of output clusters
6a8be5c3 178
6fdd30c4 179 fGeom = AliEMCALGeometry::GetInstance(fGeomName) ;
1f77507b 180 if(fOutputAODBranchName.Length()!=0){
181 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
182 fOutputAODBranch->SetName(fOutputAODBranchName);
183 AddAODBranch("TClonesArray", &fOutputAODBranch);
184 }
185 else {
186 AliFatal("fOutputAODBranchName not set\n");
187 }
188}
189
190//________________________________________________________________________
191void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
192 // Put calo cells in standard branch
193 AliVEvent * event = InputEvent();
194 AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
195 Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
196
197 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
198 aodEMcells.CreateContainer(nEMcell);
199 aodEMcells.SetType(AliVCaloCells::kEMCALCell);
200 Double_t calibFactor = 1.;
201 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
202 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
203 fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
204 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
205
206 if(fRecoUtils->IsRecalibrationOn()){
207 calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
208 }
209
210 if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
211 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
1f77507b 212 }
213 else {
214 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
1f77507b 215 }
216 }
217 aodEMcells.Sort();
218
219}
6fdd30c4 220
1f77507b 221//________________________________________________________________________
222void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
223 //Put event header information in standard AOD branch
224
225 AliVEvent* event = InputEvent();
226 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
227 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
6a8be5c3 228
1f77507b 229 Double_t pos[3] ;
230 Double_t covVtx[6];
231 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
232
233 AliAODHeader* header = AODEvent()->GetHeader();
234 header->SetRunNumber(event->GetRunNumber());
235
236 if(esdevent){
237 TTree* tree = fInputHandler->GetTree();
238 if (tree) {
239 TFile* file = tree->GetCurrentFile();
240 if (file) header->SetESDFileName(file->GetName());
241 }
242 }
9a42b255 243 else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
1f77507b 244
245 header->SetBunchCrossNumber(event->GetBunchCrossNumber());
246 header->SetOrbitNumber(event->GetOrbitNumber());
247 header->SetPeriodNumber(event->GetPeriodNumber());
248 header->SetEventType(event->GetEventType());
249
250 //Centrality
251 if(event->GetCentrality()){
252 header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
253 }
254 else{
255 header->SetCentrality(0);
256 }
257
258 //Trigger
259 header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
9a42b255 260 if (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
261 else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
1f77507b 262 header->SetTriggerMask(event->GetTriggerMask());
263 header->SetTriggerCluster(event->GetTriggerCluster());
264 if(esdevent){
265 header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
266 header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
267 header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
268 }
9a42b255 269 else if (aodevent){
1f77507b 270 header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
271 header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
272 header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
273 }
274
275 header->SetMagneticField(event->GetMagneticField());
276 //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
277
278 header->SetZDCN1Energy(event->GetZDCN1Energy());
279 header->SetZDCP1Energy(event->GetZDCP1Energy());
280 header->SetZDCN2Energy(event->GetZDCN2Energy());
281 header->SetZDCP2Energy(event->GetZDCP2Energy());
282 header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
283
284 Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
285 Float_t diamcov[3];
286 event->GetDiamondCovXY(diamcov);
287 header->SetDiamond(diamxy,diamcov);
9a42b255 288 if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
289 else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
1f77507b 290
291 //
292 //
293 Int_t nVertices = 1 ;/* = prim. vtx*/;
294 Int_t nCaloClus = event->GetNumberOfCaloClusters();
295
296 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
297
298 // Access to the AOD container of vertices
299 TClonesArray &vertices = *(AODEvent()->GetVertices());
300 Int_t jVertices=0;
301
302 // Add primary vertex. The primary tracks will be defined
303 // after the loops on the composite objects (V0, cascades, kinks)
304 event->GetPrimaryVertex()->GetXYZ(pos);
305 Float_t chi = 0;
306 if (esdevent){
307 esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
308 chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
309 }
310 else if (aodevent){
311 aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
312 chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
313 }
314
315 AliAODVertex * primary = new(vertices[jVertices++])
316 AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
317 primary->SetName(event->GetPrimaryVertex()->GetName());
318 primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
319
6fdd30c4 320}
321
322//________________________________________________________________________
323void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
324{
325 // Main loop
326 // Called for each event
6a8be5c3 327
220447cd 328 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
329 Int_t eventN = Entry();
330 if(aodIH) eventN = aodIH->GetReadEntry();
331
1f77507b 332 if (eventN > fMaxEvent) return;
c50e6949 333 //printf("Clusterizer --- Event %d-- \n",eventN);
6a8be5c3 334
6fdd30c4 335 //Remove the contents of output list set in the previous event
336 fOutputAODBranch->Clear("C");
337
1f77507b 338 //Magic line to write events to AOD file
339 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
340 LoadBranches();
341
e9dd2d80 342 //Init pointers, clusterizer, ocdb
343 if(fAccessOCDB) AccessOCDB();
344 InitClusterization();
6a8be5c3 345
6060ed91 346 //Get the event
347 AliVEvent * event = 0;
348 AliESDEvent * esdevent = 0;
6a8be5c3 349
350 //Fill output event with header
351
6060ed91 352 //Check if input event are embedded events
353 //If so, take output event
6060ed91 354 if (aodIH && aodIH->GetMergeEvents()) {
e9dd2d80 355 //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
6060ed91 356 event = AODEvent();
1f77507b 357
358 if(!aodIH->GetMergeEMCALCells())
359 AliFatal("Events merged but not EMCAL cells, check analysis settings!");
360
6a8be5c3 361 // printf("InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
362 // InputEvent()->GetEMCALCells()->GetNumberOfCells());
363 // printf("MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
364 // aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
365 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
366 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
367 // if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E %f\n",icl,sigCluster->E());
368 // }
369 // printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
370 // AODEvent()->GetEMCALCells()->GetNumberOfCells());
1f77507b 371
6060ed91 372 }
373 else {
374 event = InputEvent();
375 esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
1f77507b 376 if(fFillAODCaloCells) FillAODCaloCells();
377 if(fFillAODHeader) FillAODHeader();
6060ed91 378 }
379
6fdd30c4 380 if (!event) {
3fa4fb85 381 Error("UserExec","Event not available");
6fdd30c4 382 return;
383 }
6a8be5c3 384
6fdd30c4 385 //-------------------------------------------------------------------------------------
386 //Set the geometry matrix, for the first event, skip the rest
387 //-------------------------------------------------------------------------------------
388 if(!fGeomMatrixSet){
389 if(fLoadGeomMatrices){
390 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
391 if(fGeomMatrix[mod]){
392 if(DebugLevel() > 1)
393 fGeomMatrix[mod]->Print();
394 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
395 }
396 fGeomMatrixSet=kTRUE;
397 }//SM loop
398 }//Load matrices
399 else if(!gGeoManager){
3fa4fb85 400 Info("UserExec","Get geo matrices from data");
6fdd30c4 401 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
402 if(!strcmp(event->GetName(),"AliAODEvent")) {
403 if(DebugLevel() > 1)
3fa4fb85 404 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
6fdd30c4 405 }//AOD
406 else {
3fa4fb85 407 if(DebugLevel() > 1)
408 Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
6fdd30c4 409 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
410 if(!esd) {
3fa4fb85 411 Error("UserExec","This event does not contain ESDs?");
6fdd30c4 412 return;
413 }
414 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
415 if(DebugLevel() > 1)
416 esd->GetEMCALMatrix(mod)->Print();
417 if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
418 }
419 fGeomMatrixSet=kTRUE;
420 }//ESD
421 }//Load matrices from Data
5994e71f 422
423 //Recover time dependent corrections, put then in recalibration histograms. Do it once
424 fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
6a8be5c3 425
6fdd30c4 426 }//first event
6a8be5c3 427
6fdd30c4 428 //Get the list of cells needed for unfolding and reclustering
6a8be5c3 429 AliVCaloCells *cells = event->GetEMCALCells();
430 Int_t nClustersOrg = 0;
431 Double_t cellAmplitude = 0;
432 Double_t cellTime = 0;
433 Short_t cellNumber = 0;
434
6fdd30c4 435 //-------------------------------------------
436 //---------Unfolding clusters----------------
437 //-------------------------------------------
438 if (fJustUnfold) {
439
440 //Fill the array with the EMCAL clusters, copy them
441 for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
442 {
443 AliVCluster *clus = event->GetCaloCluster(i);
444 if(clus->IsEMCAL()){
5994e71f 445
446 //recalibrate/remove bad channels/etc if requested
447 if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
5994e71f 448 continue;
449 }
450
451 if(fRecoUtils->IsRecalibrationOn()){
6a8be5c3 452
453 //Calibrate cluster
5994e71f 454 //printf("Energy before %f ",clus->E());
455 fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
456 //printf("; after %f\n",clus->E());
6a8be5c3 457
458 //CalibrateCells
459 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
460 {
461 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
462 break;
463
464
465 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
466 fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
467 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
468
469 //Do not include bad channels found in analysis?
470 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
471 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
472 fCellLabels[cellNumber] =-1; //reset the entry in the array for next event
473 fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
474 continue;
475 }
476
477 cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
478
479 }// cells loop
480 }// recalibrate
481
5994e71f 482 //Cast to ESD or AOD, needed to create the cluster array
e615d952 483 AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
484 AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
485 if (esdCluster){
486 fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
6fdd30c4 487 }//ESD
e615d952 488 else if(aodCluster){
489 fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
6fdd30c4 490 }//AOD
3fa4fb85 491 else
492 Warning("UserExec()"," - Wrong CaloCluster type?");
6fdd30c4 493 nClustersOrg++;
494 }
495 }
496
497 //Do the unfolding
498 fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
499
500 //CLEAN-UP
501 fUnfolder->Clear();
502
6a8be5c3 503 } // just unfold ESD/AOD cluster
504
6fdd30c4 505 //-------------------------------------------
506 //---------- Recluster cells ----------------
507 //-------------------------------------------
508
509 else{
510 //-------------------------------------------------------------------------------------
511 //Transform CaloCells into Digits
512 //-------------------------------------------------------------------------------------
6060ed91 513
514 //In case of MC, first loop on the clusters and fill MC label to array
515 //.....................................................................
516
6a8be5c3 517 // for(Int_t j = 0; j < 24*48*11; j++) {
518 // if(fCellLabels[j] !=-1) printf("Not well initialized cell %d, label1 %d\n",j,fCellLabels[j] ) ;
519 // if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
520 // }
6060ed91 521
1f77507b 522 Int_t nClusters = event->GetNumberOfCaloClusters();
006e9636 523 if(aodIH && aodIH->GetEventToMerge()) //Embedding
1f77507b 524 nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
525 for (Int_t i = 0; i < nClusters; i++)
6060ed91 526 {
1f77507b 527 AliVCluster *clus = 0;
006e9636 528 if(aodIH && aodIH->GetEventToMerge()) //Embedding
529 clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
530 else
531 clus = event->GetCaloCluster(i);
6a8be5c3 532
1f77507b 533 if(!clus) {
534 printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
535 return;
536 }
537
6060ed91 538 if(clus->IsEMCAL()){
1f77507b 539 //printf("Cluster Signal %d, energy %f\n",clus->GetID(),clus->E());
6060ed91 540 Int_t label = clus->GetLabel();
1f77507b 541 Int_t label2 = -1 ;
542 if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
6060ed91 543 UShort_t * index = clus->GetCellsAbsId() ;
544 for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
220447cd 545 fCellLabels[index[icell]] =label;
1f77507b 546 fCellSecondLabels[index[icell]]=label2;
abf0860a 547 //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
6060ed91 548 }
549 nClustersOrg++;
550 }
551 }
6a8be5c3 552
6060ed91 553 // Create digits
554 //.................
6fdd30c4 555 Int_t idigit = 0;
556 Int_t id = -1;
557 Float_t amp = -1;
558 Float_t time = -1;
559
6fdd30c4 560 TTree *digitsTree = new TTree("digitstree","digitstree");
561 digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
562
563 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
564 {
565 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
566 break;
567
568 time = cellTime;
569 amp = cellAmplitude;
570 id = cellNumber;
571
5994e71f 572 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
573 fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
574 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
575
576 //Do not include bad channels found in analysis?
577 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
6a8be5c3 578 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
579 fCellLabels[id] =-1; //reset the entry in the array for next event
580 fCellSecondLabels[id]=-1; //reset the entry in the array for next event
5994e71f 581 //printf("Remove channel %d\n",id);
582 continue;
583 }
6a8be5c3 584
5994e71f 585 //Recalibrate?
586 if(fRecoUtils->IsRecalibrationOn()){
587 //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
588 amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
589 }
6a8be5c3 590
6060ed91 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);
593 //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
594 //else printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
220447cd 595 fCellLabels[id] =-1; //reset the entry in the array for next event
5994e71f 596
597 //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
6fdd30c4 598 idigit++;
599 }
600
601 //Fill the tree with digits
602 digitsTree->Fill();
603
604 //-------------------------------------------------------------------------------------
605 //Do the clusterization
606 //-------------------------------------------------------------------------------------
607 TTree *clustersTree = new TTree("clustertree","clustertree");
6a8be5c3 608
6fdd30c4 609 fClusterizer->SetInput(digitsTree);
610 fClusterizer->SetOutput(clustersTree);
611 fClusterizer->Digits2Clusters("");
6a8be5c3 612
6fdd30c4 613 //-------------------------------------------------------------------------------------
614 //Transform the recpoints into AliVClusters
615 //-------------------------------------------------------------------------------------
616
617 clustersTree->SetBranchStatus("*",0); //disable all branches
618 clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
619
620 TBranch *branch = clustersTree->GetBranch("EMCALECARP");
621 branch->SetAddress(&fClusterArr);
622 branch->GetEntry(0);
6a8be5c3 623
6fdd30c4 624 RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
625
6a8be5c3 626 if(!fCaloClusterArr){
627 printf("AliAnalisysTaskEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
628 return;
629 }
630
631 if(fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
9f15c23c 632 printf("AliAnalisysTaskEMCALClusterize::UserExec() - Some RecRoints not transformed into CaloClusters: Input entries %d - Output entries %d - %d (not fast)\n",
6a8be5c3 633 fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
634 }
635
01313069 636 //Reset the array with second labels for this event
637 memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
638
6fdd30c4 639 //---CLEAN UP-----
640 fClusterizer->Clear();
641 fDigitsArr ->Clear("C");
642 fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
6a8be5c3 643
6fdd30c4 644 clustersTree->Delete("all");
645 digitsTree ->Delete("all");
6fdd30c4 646 }
647
385b7abf 648 //Recalculate track-matching for the new clusters, only with ESDs
cb3001b9 649 if(fDoTrackMatching) fRecoUtils->FindMatches(event,fCaloClusterArr,fGeom);
6a8be5c3 650
385b7abf 651
6fdd30c4 652 //-------------------------------------------------------------------------------------
653 //Put the new clusters in the AOD list
654 //-------------------------------------------------------------------------------------
655
dd9a9f63 656 Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
6060ed91 657 //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
6fdd30c4 658 for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
659 AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
3fa4fb85 660 //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
385b7abf 661
662 //Add matched track, if any, only with ESDs
cb3001b9 663 if(esdevent && fDoTrackMatching){
385b7abf 664 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
665 if(trackIndex >= 0){
666 newCluster->AddTrackMatched(event->GetTrack(trackIndex));
667 if(DebugLevel() > 1)
668 Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
669 }
670 }
0821ba9f 671
5994e71f 672 //In case of new bad channels, recalculate distance to bad channels
673 if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
674 //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
675 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
676 //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
677 }
678
6a8be5c3 679 // if(newCluster->GetNLabels()>0){
680 // printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
681 // UShort_t * newindex = newCluster->GetCellsAbsId() ;
682 // for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
683 // printf("\t absID %d\n",newindex[icell]);
684 // }
685 // }
1f77507b 686
6a8be5c3 687 // printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
688 // printf("labels: ");
689 // for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
690 // printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab));
691 // printf("\n ");
6060ed91 692
0821ba9f 693 newCluster->SetID(i);
6fdd30c4 694 new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
695 }
696
1f77507b 697 //if(fOutputAODBranchName.Length()!=0)
698 fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
699
6fdd30c4 700 //---CLEAN UP-----
701 fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
3fa4fb85 702}
6fdd30c4 703
704//_____________________________________________________________________
bd9e8ebd 705Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
6fdd30c4 706{
707 //Access to OCDB stuff
e9dd2d80 708
6fdd30c4 709 AliVEvent * event = InputEvent();
bd9e8ebd 710 if (!event)
6fdd30c4 711 {
385b7abf 712 Warning("AccessOCDB","Event not available!!!");
bd9e8ebd 713 return kFALSE;
714 }
e9dd2d80 715
bd9e8ebd 716 if (event->GetRunNumber()==fRun)
717 return kTRUE;
718 fRun = event->GetRunNumber();
e9dd2d80 719
bd9e8ebd 720 if(DebugLevel() > 1 )
721 Info("AccessODCD()"," Begin");
385b7abf 722
e9dd2d80 723 //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
385b7abf 724
bd9e8ebd 725 AliCDBManager *cdb = AliCDBManager::Instance();
e9dd2d80 726
727
385b7abf 728 if (fOCDBpath.Length()){
6fdd30c4 729 cdb->SetDefaultStorage(fOCDBpath.Data());
385b7abf 730 Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
731 }
732
bd9e8ebd 733 cdb->SetRun(event->GetRunNumber());
e9dd2d80 734
bd9e8ebd 735 //
736 // EMCAL from RAW OCDB
737 if (fOCDBpath.Contains("alien:"))
738 {
739 cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
740 cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
741 }
e9dd2d80 742
bd9e8ebd 743 TString path = cdb->GetDefaultStorage()->GetBaseFolder();
385b7abf 744
bd9e8ebd 745 // init parameters:
385b7abf 746
bd9e8ebd 747 //Get calibration parameters
748 if(!fCalibData)
749 {
750 AliCDBEntry *entry = (AliCDBEntry*)
e9dd2d80 751 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
385b7abf 752
bd9e8ebd 753 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
6fdd30c4 754 }
bd9e8ebd 755
756 if(!fCalibData)
757 AliFatal("Calibration parameters not found in CDB!");
e9dd2d80 758
bd9e8ebd 759 //Get calibration parameters
760 if(!fPedestalData)
6fdd30c4 761 {
bd9e8ebd 762 AliCDBEntry *entry = (AliCDBEntry*)
e9dd2d80 763 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
385b7abf 764
bd9e8ebd 765 if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
6fdd30c4 766 }
e9dd2d80 767
bd9e8ebd 768 if(!fPedestalData)
769 AliFatal("Dead map not found in CDB!");
385b7abf 770
6fdd30c4 771 return kTRUE;
6fdd30c4 772}
773
774//________________________________________________________________________________________
775void AliAnalysisTaskEMCALClusterize::InitClusterization()
776{
777 //Select clusterization/unfolding algorithm and set all the needed parameters
778
bd9e8ebd 779 if (fJustUnfold){
780 // init the unfolding afterburner
781 delete fUnfolder;
782 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
783 return;
6a8be5c3 784 }
785
bd9e8ebd 786 //First init the clusterizer
787 delete fClusterizer;
788 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
789 fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
cfcbe5d2 790 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
791 fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
bd9e8ebd 792 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
793 fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
cfcbe5d2 794 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerv2) { //FIX this other way.
6a8be5c3 795 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
796 clusterizer->SetNRowDiff(2);
797 clusterizer->SetNColDiff(2);
798 fClusterizer = clusterizer;
bd9e8ebd 799 } else{
800 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
801 }
6a8be5c3 802
bd9e8ebd 803 //Now set the parameters
804 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
805 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
806 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
807 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
808 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
809 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
810 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
811 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
812 fClusterizer->SetInputCalibrated ( kTRUE );
2af35150 813 fClusterizer->SetJustClusters ( kTRUE );
bd9e8ebd 814 //In case of unfolding after clusterization is requested, set the corresponding parameters
815 if(fRecParam->GetUnfold()){
816 Int_t i=0;
817 for (i = 0; i < 8; i++) {
818 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
819 }//end of loop over parameters
820 for (i = 0; i < 3; i++) {
821 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
822 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
823 }//end of loop over parameters
6fdd30c4 824
bd9e8ebd 825 fClusterizer->InitClusterUnfolding();
826 }// to unfold
6fdd30c4 827}
3fa4fb85 828
6fdd30c4 829//________________________________________________________________________________________
830void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
831{
832 // Restore clusters from recPoints
833 // Cluster energy, global position, cells and their amplitude fractions are restored
6a8be5c3 834 Int_t j = 0;
6fdd30c4 835 for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
836 {
837 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
838
839 const Int_t ncells = recPoint->GetMultiplicity();
840 Int_t ncells_true = 0;
841
6a8be5c3 842 if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
843
6fdd30c4 844 // cells and their amplitude fractions
845 UShort_t absIds[ncells];
846 Double32_t ratios[ncells];
847
848 for (Int_t c = 0; c < ncells; c++) {
849 AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
850
851 absIds[ncells_true] = digit->GetId();
852 ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
853
854 if (ratios[ncells_true] > 0.001) ncells_true++;
855 }
856
857 if (ncells_true < 1) {
6a8be5c3 858 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",recPoint->GetEnergy(), ncells);
6fdd30c4 859 continue;
860 }
861
862 TVector3 gpos;
863 Float_t g[3];
864
865 // calculate new cluster position
866 recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
867 recPoint->GetGlobalPosition(gpos);
868 gpos.GetXYZ(g);
869
870 // create a new cluster
937b6c91 871 //AliAODCaloCluster *clus = new AliAODCaloCluster();
6a8be5c3 872 (*clusArray)[j] = new AliAODCaloCluster() ;
873 AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
874 j++;
6fdd30c4 875 clus->SetType(AliVCluster::kEMCALClusterv1);
876 clus->SetE(recPoint->GetEnergy());
877 clus->SetPosition(g);
878 clus->SetNCells(ncells_true);
879 clus->SetCellsAbsId(absIds);
880 clus->SetCellsAmplitudeFraction(ratios);
881 clus->SetDispersion(recPoint->GetDispersion());
882 clus->SetChi2(-1); //not yet implemented
bd9e8ebd 883 clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
6fdd30c4 884 clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
885 Float_t elipAxis[2];
886 recPoint->GetElipsAxis(elipAxis);
887 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
888 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
5994e71f 889 clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
6060ed91 890
891 //MC
892 Int_t parentMult = 0;
893 Int_t *parentList = recPoint->GetParents(parentMult);
894 clus->SetLabel(parentList, parentMult);
895
6a8be5c3 896 // if(parentMult)printf("RecToESD: Labels %d",parentMult);
897 // for (Int_t iii = 0; iii < parentMult; iii++) {
898 // printf("\t label %d\n",parentList[iii]);
899 // }
900
37d2296c 901 //Write the second major contributor to each MC cluster.
902 Int_t iNewLabel ;
903 for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
904
905 Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
6ffa213f 906 if(idCell>=0){
907 iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
908 for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
909 {
910 if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
911 if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
912 }
913 if (iNewLabel == 1)
914 {
937b6c91 915 //Int_t * newLabelArray = (Int_t*)malloc((clus->GetNLabels()+1)*sizeof(Int_t)) ;
916 Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
6ffa213f 917 for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
918 newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
919 }
6a8be5c3 920
6ffa213f 921 newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
220447cd 922 //if(fCellSecondLabels[idCell]>-1)printf("\t new label %d\n",fCellSecondLabels[idCell]);
6ffa213f 923 clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
abf0860a 924 delete [] newLabelArray;
6ffa213f 925 }
926 }//positive cell number
37d2296c 927 }
928
6fdd30c4 929 } // recPoints loop
6a8be5c3 930
6fdd30c4 931}