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