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