1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
17 //_________________________________________________________________________
18 // This analysis provides a new list of clusters to be used in other analysis
20 // Author: Gustavo Conesa Balbastre,
21 // Adapted from analysis class from Deepa Thomas
24 //_________________________________________________________________________
28 #include "TRefArray.h"
29 #include "TClonesArray.h"
31 #include "TGeoManager.h"
33 #include "TInterpreter.h"
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"
48 #include "AliVEventHandler.h"
49 #include "AliAODInputHandler.h"
52 #include "AliEMCALRecParam.h"
53 #include "AliEMCALAfterBurnerUF.h"
54 #include "AliEMCALGeometry.h"
55 #include "AliEMCALClusterizerNxN.h"
56 #include "AliEMCALClusterizerv1.h"
57 #include "AliEMCALClusterizerv2.h"
58 #include "AliEMCALRecPoint.h"
59 #include "AliEMCALDigit.h"
60 #include "AliCaloCalibPedestal.h"
61 #include "AliEMCALCalibData.h"
62 #include "AliEMCALRecoUtils.h"
64 #include "AliAnalysisTaskEMCALClusterize.h"
66 ClassImp(AliAnalysisTaskEMCALClusterize)
68 //________________________________________________________________________
69 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
70 : AliAnalysisTaskSE(name)
71 , fGeom(0), fGeomName("EMCAL_COMPLETEV1")
72 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
73 , fCalibData(0), fPedestalData(0)
74 , fOCDBpath("raw://"), fAccessOCDB(kFALSE)
75 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
76 , fRecParam(0), fClusterizer(0)
77 , fUnfolder(0), fJustUnfold(kFALSE)
78 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
79 , fFillAODFile(kTRUE), fFillAODHeader(0)
80 , fFillAODCaloCells(0), fRun(-1)
81 , fRecoUtils(0), fConfigName("")
82 , fCellLabels(), fCellSecondLabels(), fCellTime()
83 , fMaxEvent(1000000000), fDoTrackMatching(kFALSE)
84 , fSelectCell(kFALSE), fSelectCellMinE(0.005), fSelectCellMinFrac(0.001)
85 , fRemoveLEDEvents(kFALSE)
88 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
89 for(Int_t j = 0; j < 24*48*11; j++) {
91 fCellSecondLabels[j] = -1;
95 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
96 fClusterArr = new TObjArray(100);
97 fCaloClusterArr = new TObjArray(1000);
98 fRecParam = new AliEMCALRecParam;
99 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
100 fRecoUtils = new AliEMCALRecoUtils();
104 //________________________________________________________________________
105 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
106 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
107 , fGeom(0), fGeomName("EMCAL_COMPLETEV1")
108 , fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
109 , fCalibData(0), fPedestalData(0)
110 , fOCDBpath("raw://"), fAccessOCDB(kFALSE)
111 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
112 , fRecParam(0), fClusterizer(0)
113 , fUnfolder(0), fJustUnfold(kFALSE)
114 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters")
115 , fFillAODFile(kTRUE), fFillAODHeader(0)
116 , fFillAODCaloCells(0), fRun(-1)
117 , fRecoUtils(0), fConfigName("")
118 , fCellLabels(), fCellSecondLabels(), fCellTime()
119 , fMaxEvent(1000000000), fDoTrackMatching(kFALSE)
120 , fSelectCell(kFALSE), fSelectCellMinE(0.005), fSelectCellMinFrac(0.001)
121 , fRemoveLEDEvents(kFALSE)
125 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0;
126 for(Int_t j = 0; j < 24*48*11; j++) {
128 fCellSecondLabels[j] = -1;
131 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
132 fClusterArr = new TObjArray(100);
133 fCaloClusterArr = new TObjArray(100);
134 fRecParam = new AliEMCALRecParam;
135 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
136 fRecoUtils = new AliEMCALRecoUtils();
140 //________________________________________________________________________
141 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
146 fDigitsArr->Clear("C");
151 fClusterArr->Delete();
155 if (fCaloClusterArr){
156 fCaloClusterArr->Delete();
157 delete fCaloClusterArr;
160 if(fClusterizer) delete fClusterizer;
161 if(fUnfolder) delete fUnfolder;
162 if(fRecoUtils) delete fRecoUtils;
166 //-------------------------------------------------------------------
167 void AliAnalysisTaskEMCALClusterize::Init()
169 //Init analysis with configuration macro if available
171 if(gROOT->LoadMacro(fConfigName) >=0){
172 printf("Configure analysis with %s\n",fConfigName.Data());
173 AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
174 fGeomName = clus->fGeomName;
175 fLoadGeomMatrices = clus->fLoadGeomMatrices;
176 fOCDBpath = clus->fOCDBpath;
177 fAccessOCDB = clus->fAccessOCDB;
178 fRecParam = clus->fRecParam;
179 fJustUnfold = clus->fJustUnfold;
180 fFillAODFile = clus->fFillAODFile;
181 fRecoUtils = clus->fRecoUtils;
182 fConfigName = clus->fConfigName;
183 fMaxEvent = clus->fMaxEvent;
184 fDoTrackMatching = clus->fDoTrackMatching;
185 fOutputAODBranchName = clus->fOutputAODBranchName;
186 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
192 //-------------------------------------------------------------------
193 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
195 // Init geometry, create list of output clusters
197 fGeom = AliEMCALGeometry::GetInstance(fGeomName) ;
198 if(fOutputAODBranchName.Length()!=0){
199 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
200 fOutputAODBranch->SetName(fOutputAODBranchName);
201 AddAODBranch("TClonesArray", &fOutputAODBranch);
204 AliFatal("fOutputAODBranchName not set\n");
208 //________________________________________________________________________
209 void AliAnalysisTaskEMCALClusterize::FillAODCaloCells() {
210 // Put calo cells in standard branch
211 AliVEvent * event = InputEvent();
212 AliVCaloCells &eventEMcells = *(event->GetEMCALCells());
213 Int_t nEMcell = eventEMcells.GetNumberOfCells() ;
215 AliAODCaloCells &aodEMcells = *(AODEvent()->GetEMCALCells());
216 aodEMcells.CreateContainer(nEMcell);
217 aodEMcells.SetType(AliVCaloCells::kEMCALCell);
218 Double_t calibFactor = 1.;
219 for (Int_t iCell = 0; iCell < nEMcell; iCell++) {
220 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
221 fGeom->GetCellIndex(eventEMcells.GetCellNumber(iCell),imod,iTower,iIphi,iIeta);
222 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
224 if(fRecoUtils->IsRecalibrationOn()){
225 calibFactor = fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
228 if(!fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){ //Channel is not declared as bad
229 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),eventEMcells.GetAmplitude(iCell)*calibFactor);
232 aodEMcells.SetCell(iCell,eventEMcells.GetCellNumber(iCell),0);
239 //________________________________________________________________________
240 void AliAnalysisTaskEMCALClusterize::FillAODHeader() {
241 //Put event header information in standard AOD branch
243 AliVEvent* event = InputEvent();
244 AliESDEvent* esdevent = dynamic_cast<AliESDEvent*> (event);
245 AliAODEvent* aodevent = dynamic_cast<AliAODEvent*> (event);
249 for (Int_t i = 0; i < 6; i++) covVtx[i] = 0.;
251 AliAODHeader* header = AODEvent()->GetHeader();
252 header->SetRunNumber(event->GetRunNumber());
255 TTree* tree = fInputHandler->GetTree();
257 TFile* file = tree->GetCurrentFile();
258 if (file) header->SetESDFileName(file->GetName());
261 else if (aodevent) header->SetESDFileName(aodevent->GetHeader()->GetESDFileName());
263 header->SetBunchCrossNumber(event->GetBunchCrossNumber());
264 header->SetOrbitNumber(event->GetOrbitNumber());
265 header->SetPeriodNumber(event->GetPeriodNumber());
266 header->SetEventType(event->GetEventType());
269 if(event->GetCentrality()){
270 header->SetCentrality(new AliCentrality(*(event->GetCentrality())));
273 header->SetCentrality(0);
277 header->SetOfflineTrigger(fInputHandler->IsEventSelected()); // propagate the decision of the physics selection
278 if (esdevent) header->SetFiredTriggerClasses(esdevent->GetFiredTriggerClasses());
279 else if (aodevent) header->SetFiredTriggerClasses(aodevent->GetFiredTriggerClasses());
280 header->SetTriggerMask(event->GetTriggerMask());
281 header->SetTriggerCluster(event->GetTriggerCluster());
283 header->SetL0TriggerInputs(esdevent->GetHeader()->GetL0TriggerInputs());
284 header->SetL1TriggerInputs(esdevent->GetHeader()->GetL1TriggerInputs());
285 header->SetL2TriggerInputs(esdevent->GetHeader()->GetL2TriggerInputs());
288 header->SetL0TriggerInputs(aodevent->GetHeader()->GetL0TriggerInputs());
289 header->SetL1TriggerInputs(aodevent->GetHeader()->GetL1TriggerInputs());
290 header->SetL2TriggerInputs(aodevent->GetHeader()->GetL2TriggerInputs());
293 header->SetMagneticField(event->GetMagneticField());
294 //header->SetMuonMagFieldScale(esdevent->GetCurrentDip()/6000.);
296 header->SetZDCN1Energy(event->GetZDCN1Energy());
297 header->SetZDCP1Energy(event->GetZDCP1Energy());
298 header->SetZDCN2Energy(event->GetZDCN2Energy());
299 header->SetZDCP2Energy(event->GetZDCP2Energy());
300 header->SetZDCEMEnergy(event->GetZDCEMEnergy(0),event->GetZDCEMEnergy(1));
302 Float_t diamxy[2]={event->GetDiamondX(),event->GetDiamondY()};
304 event->GetDiamondCovXY(diamcov);
305 header->SetDiamond(diamxy,diamcov);
306 if (esdevent) header->SetDiamondZ(esdevent->GetDiamondZ(),esdevent->GetSigma2DiamondZ());
307 else if (aodevent) header->SetDiamondZ(aodevent->GetDiamondZ(),aodevent->GetSigma2DiamondZ());
311 Int_t nVertices = 1 ;/* = prim. vtx*/;
312 Int_t nCaloClus = event->GetNumberOfCaloClusters();
314 AODEvent()->ResetStd(0, nVertices, 0, 0, 0, nCaloClus, 0, 0);
316 // Access to the AOD container of vertices
317 TClonesArray &vertices = *(AODEvent()->GetVertices());
320 // Add primary vertex. The primary tracks will be defined
321 // after the loops on the composite objects (V0, cascades, kinks)
322 event->GetPrimaryVertex()->GetXYZ(pos);
325 esdevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
326 chi = esdevent->GetPrimaryVertex()->GetChi2toNDF();
329 aodevent->GetPrimaryVertex()->GetCovMatrix(covVtx);
330 chi = aodevent->GetPrimaryVertex()->GetChi2perNDF();//Different from ESD?
333 AliAODVertex * primary = new(vertices[jVertices++])
334 AliAODVertex(pos, covVtx, chi, NULL, -1, AliAODVertex::kPrimary);
335 primary->SetName(event->GetPrimaryVertex()->GetName());
336 primary->SetTitle(event->GetPrimaryVertex()->GetTitle());
340 //________________________________________________________________________
341 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
344 // Called for each event
346 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
347 Int_t eventN = Entry();
348 if(aodIH) eventN = aodIH->GetReadEntry();
350 if (eventN > fMaxEvent) {
351 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
355 //printf("Clusterizer --- Event %d-- \n",eventN);
357 //Remove the contents of output list set in the previous event
358 fOutputAODBranch->Clear("C");
362 //Init pointers, clusterizer, ocdb
363 if(fAccessOCDB) AccessOCDB();
364 InitClusterization();
367 AliVEvent * event = 0;
368 AliESDEvent * esdevent = 0;
370 //Fill output event with header
372 //Check if input event are embedded events
373 //If so, take output event
374 if (aodIH && aodIH->GetMergeEvents()) {
375 //printf("AliAnalysisTaskEMCALClusterize::UserExec() - Use embedded events\n");
378 if(!aodIH->GetMergeEMCALCells())
379 AliFatal("Events merged but not EMCAL cells, check analysis settings!");
381 // printf("InputEvent N Clusters %d, N Cells %d\n",InputEvent()->GetNumberOfCaloClusters(),
382 // InputEvent()->GetEMCALCells()->GetNumberOfCells());
383 // printf("MergedEvent N Clusters %d, N Cells %d\n",aodIH->GetEventToMerge()->GetNumberOfCaloClusters(),
384 // aodIH->GetEventToMerge()->GetEMCALCells()->GetNumberOfCells());
385 // for (Int_t icl=0; icl < aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); icl++) {
386 // AliAODCaloCluster *sigCluster = aodIH->GetEventToMerge()->GetCaloCluster(icl);
387 // if(sigCluster->IsEMCAL()) printf("Signal cluster: i %d, E %f\n",icl,sigCluster->E());
389 // printf("OutputEvent N Clusters %d, N Cells %d\n", AODEvent()->GetNumberOfCaloClusters(),
390 // AODEvent()->GetEMCALCells()->GetNumberOfCells());
394 event = InputEvent();
395 esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
396 if(fFillAODCaloCells) FillAODCaloCells();
397 if(fFillAODHeader) FillAODHeader();
401 Error("UserExec","Event not available");
402 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
406 //-------------------------------------------------------------------------------------
407 // Reject event if large clusters with large energy
408 // Use only for LHC11a data for the moment, and if input is clusterizer V1 or V1+unfolding
409 // If clusterzer NxN or V2 it does not help
410 //-------------------------------------------------------------------------------------
411 if(fRemoveLEDEvents){
412 for (Int_t i = 0; i < InputEvent()->GetNumberOfCaloClusters(); i++)
414 AliVCluster *clus = InputEvent()->GetCaloCluster(i);
417 if ((clus->E() > 500 && clus->GetNCells() > 200 ) || clus->GetNCells() > 200) {
418 Int_t absID = clus->GetCellsAbsId()[0];
419 Int_t sm = fGeom->GetSuperModuleNumber(absID);
420 printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster : E %f, ncells %d, absId(0) %d, SM %d\n",(Int_t)Entry(), clus->E(), clus->GetNCells(),absID, sm);
421 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
427 // Count number of cells with energy larger than 0.1 in SM3, cut on this number
430 for(Int_t icell = 0; icell < event->GetEMCALCells()->GetNumberOfCells(); icell++){
431 if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==3) ncellsSM3++;
432 if(event->GetEMCALCells()->GetAmplitude(icell) > 0.1 && event->GetEMCALCells()->GetCellNumber(icell)/(24*48)==4) ncellsSM4++;
435 TString triggerclasses = "";
436 if(esdevent) triggerclasses = esdevent ->GetFiredTriggerClasses();
437 else triggerclasses = ((AliAODEvent*)event)->GetFiredTriggerClasses();
440 if(triggerclasses.Contains("EMC")) ncellcut = 35;
442 if( ncellsSM3 >= ncellcut || ncellsSM4 >= 100 ){
443 printf("AliAnalysisTaskEMCALClusterize - reject event %d with cluster - reject event with ncells in SM3 %d and SM4 %d\n",(Int_t)Entry(),ncellsSM3, ncellsSM4);
444 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
448 }// Remove LED events
450 //Magic line to write events to AOD filem put after event rejection
451 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
454 //-------------------------------------------------------------------------------------
455 //Set the geometry matrix, for the first event, skip the rest
456 //-------------------------------------------------------------------------------------
458 if(fLoadGeomMatrices){
459 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
460 if(fGeomMatrix[mod]){
462 fGeomMatrix[mod]->Print();
463 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
465 fGeomMatrixSet=kTRUE;
468 else if(!gGeoManager){
469 Info("UserExec","Get geo matrices from data");
470 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
471 if(!strcmp(event->GetName(),"AliAODEvent")) {
473 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
477 Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
478 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
480 Error("UserExec","This event does not contain ESDs?");
481 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
484 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
486 esd->GetEMCALMatrix(mod)->Print();
487 if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
489 fGeomMatrixSet=kTRUE;
491 }//Load matrices from Data
493 //Recover time dependent corrections, put then in recalibration histograms. Do it once
494 fRecoUtils->SetRunDependentCorrections(InputEvent()->GetRunNumber());
498 //Get the list of cells needed for unfolding and reclustering
499 AliVCaloCells *cells = event->GetEMCALCells();
500 Int_t nClustersOrg = 0;
501 Double_t cellAmplitude = 0;
502 Double_t cellTime = 0;
503 Short_t cellNumber = 0;
505 //-------------------------------------------
506 //---------Unfolding clusters----------------
507 //-------------------------------------------
510 //Fill the array with the EMCAL clusters, copy them
511 for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
513 AliVCluster *clus = event->GetCaloCluster(i);
516 //recalibrate/remove bad channels/etc if requested
517 if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
521 if(fRecoUtils->IsRecalibrationOn()){
524 //printf("Energy before %f ",clus->E());
525 fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
526 //printf("; after %f\n",clus->E());
529 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
531 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
535 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
536 fGeom->GetCellIndex(cellNumber,imod,iTower,iIphi,iIeta);
537 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
539 //Do not include bad channels found in analysis?
540 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
541 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
542 fCellLabels[cellNumber] =-1; //reset the entry in the array for next event
543 fCellSecondLabels[cellNumber]=-1; //reset the entry in the array for next event
544 fCellTime[cellNumber] = 0.;
548 cells->SetCell(icell, cellNumber, cellAmplitude*fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),cellTime);
553 //Cast to ESD or AOD, needed to create the cluster array
554 AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
555 AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
557 fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
560 fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
563 Warning("UserExec()"," - Wrong CaloCluster type?");
569 fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
574 } // just unfold ESD/AOD cluster
576 //-------------------------------------------
577 //---------- Recluster cells ----------------
578 //-------------------------------------------
581 //-------------------------------------------------------------------------------------
582 //Transform CaloCells into Digits
583 //-------------------------------------------------------------------------------------
585 //In case of MC, first loop on the clusters and fill MC label to array
586 //.....................................................................
588 // for(Int_t j = 0; j < 24*48*11; j++) {
589 // if(fCellLabels[j] !=-1) printf("Not well initialized cell %d, label1 %d\n",j,fCellLabels[j] ) ;
590 // if(fCellSecondLabels[j]!=-1) printf("Not well initialized cell %d, label2 %d\n",j,fCellSecondLabels[j]) ;
593 Int_t nClusters = event->GetNumberOfCaloClusters();
594 if(aodIH && aodIH->GetEventToMerge()) //Embedding
595 nClusters = aodIH->GetEventToMerge()->GetNumberOfCaloClusters(); //Get clusters directly from embedded signal
596 for (Int_t i = 0; i < nClusters; i++)
598 AliVCluster *clus = 0;
599 if(aodIH && aodIH->GetEventToMerge()) //Embedding
600 clus = aodIH->GetEventToMerge()->GetCaloCluster(i); //Get clusters directly from embedded signal
602 clus = event->GetCaloCluster(i);
605 printf("AliEMCALReclusterize::UserExec() - No Clusters\n");
606 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
611 //printf("Cluster Signal %d, energy %f\n",clus->GetID(),clus->E());
612 Int_t label = clus->GetLabel();
614 if (clus->GetNLabels()>=2) label2 = clus->GetLabelAt(1) ;
615 UShort_t * index = clus->GetCellsAbsId() ;
616 for(Int_t icell=0; icell < clus->GetNCells(); icell++ ){
617 fCellLabels[index[icell]] = label;
618 fCellSecondLabels[index[icell]] = label2;
619 //printf("Clusterizer in : TOF %g\n",clus->GetTOF()*1.e9);
621 fCellTime[icell] = clus->GetTOF();
623 //printf("1) absID %d, label[0] %d label[1] %d\n",index[icell], fCellLabels[index[icell]],fCellSecondLabels[index[icell]]);
636 TTree *digitsTree = new TTree("digitstree","digitstree");
637 digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
639 Int_t bc = InputEvent()->GetBunchCrossNumber();
641 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
643 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
650 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
651 fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
652 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
654 //Do not include bad channels found in analysis?
655 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
656 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
657 fCellLabels[id] =-1; //reset the entry in the array for next event
658 fCellSecondLabels[id]=-1; //reset the entry in the array for next event
660 //printf("Remove channel %d\n",id);
665 if(fRecoUtils->IsRecalibrationOn()){
666 //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
667 amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
670 // In case of AOD analysis cell time is 0, approximate replacing by time of the cluster the digit belongs.
671 if (time*1e9 < 1.) time = fCellTime[id];
674 fRecoUtils->RecalibrateCellTime(id,bc,time);
676 // printf("Clusterizer: Id %d, Time org %e, Time new %e; Amp org %f, Amp new %f\n",
677 // id, cells->GetTime(icell),time, cells->GetAmplitude(icell),amp);
680 //Create the digit, put a fake primary deposited energy to trick the clusterizer when checking the most likely primary
681 new((*fDigitsArr)[idigit]) AliEMCALDigit( fCellLabels[id], fCellLabels[id],id, amp, time,AliEMCALDigit::kHG,idigit, 0, 0, 1);
682 //if(fCellLabels[id]>=0)printf("2) Digit cell %d, label %d\n",id,fCellLabels[id]) ;
683 //else printf("2) Digit cell %d, no label, amp %f \n",id,amp) ;
684 fCellLabels[id] =-1; //reset the entry in the array for next event
686 //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
690 //Fill the tree with digits
693 //-------------------------------------------------------------------------------------
694 //Do the clusterization
695 //-------------------------------------------------------------------------------------
696 TTree *clustersTree = new TTree("clustertree","clustertree");
698 fClusterizer->SetInput(digitsTree);
699 fClusterizer->SetOutput(clustersTree);
700 fClusterizer->Digits2Clusters("");
702 //-------------------------------------------------------------------------------------
703 //Transform the recpoints into AliVClusters
704 //-------------------------------------------------------------------------------------
706 clustersTree->SetBranchStatus("*",0); //disable all branches
707 clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
709 TBranch *branch = clustersTree->GetBranch("EMCALECARP");
710 branch->SetAddress(&fClusterArr);
713 RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
715 if(!fCaloClusterArr){
716 printf("AliAnalisysTaskEMCALClusterize::UserExec() - No array with CaloClusters, input RecPoints entries %d\n",fClusterArr->GetEntriesFast());
717 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(kFALSE);
721 if( DebugLevel() > 0 && fCaloClusterArr->GetEntriesFast() != fClusterArr->GetEntriesFast()){
722 printf("AliAnalisysTaskEMCALClusterize::UserExec() - Some RecRoints not transformed into CaloClusters (clusterizer %d, unfold %d): Input entries %d - Output entries %d - %d (not fast)\n",
723 fRecParam->GetClusterizerFlag(),fRecParam->GetUnfold(),
724 fClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntriesFast(), fCaloClusterArr->GetEntries());
727 //Reset the array with second labels for this event
728 memset(fCellSecondLabels, -1, sizeof(fCellSecondLabels));
731 fClusterizer->Clear();
732 fDigitsArr ->Clear("C");
733 fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
735 clustersTree->Delete("all");
736 digitsTree ->Delete("all");
739 //Recalculate track-matching for the new clusters, only with ESDs
740 if(fDoTrackMatching) fRecoUtils->FindMatches(event,fCaloClusterArr,fGeom);
743 //-------------------------------------------------------------------------------------
744 //Put the new clusters in the AOD list
745 //-------------------------------------------------------------------------------------
747 Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntriesFast();
748 //printf("New clusters %d, Org clusters %d\n",kNumberOfCaloClusters, nClustersOrg);
749 for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
750 AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
751 //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
754 if(fDoTrackMatching){
755 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
757 newCluster->AddTrackMatched(event->GetTrack(trackIndex));
759 Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
763 //In case of new bad channels, recalculate distance to bad channels
764 if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
765 //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
766 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
767 //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
770 // if(newCluster->GetNLabels()>0){
771 // printf("3) MC: N labels %d, label %d ; ", newCluster->GetNLabels(), newCluster->GetLabel() );
772 // UShort_t * newindex = newCluster->GetCellsAbsId() ;
773 // for(Int_t icell=0; icell < newCluster->GetNCells(); icell++ ){
774 // printf("\t absID %d\n",newindex[icell]);
778 // printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
779 // printf("labels: ");
780 // for(UInt_t ilab = 0; ilab < newCluster->GetNLabels(); ilab++)
781 // printf("Label[%d]: %d ",ilab, newCluster->GetLabelAt(ilab));
784 newCluster->SetID(i);
785 new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
788 //if(fOutputAODBranchName.Length()!=0)
789 fOutputAODBranch->Expand(kNumberOfCaloClusters); // resize TObjArray to 'remove' slots
792 fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
795 //_____________________________________________________________________
796 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
798 //Access to OCDB stuff
800 AliVEvent * event = InputEvent();
803 Warning("AccessOCDB","Event not available!!!");
807 if (event->GetRunNumber()==fRun)
809 fRun = event->GetRunNumber();
811 if(DebugLevel() > 1 )
812 Info("AccessODCD()"," Begin");
814 //fGeom = AliEMCALGeometry::GetInstance(fGeomName);
816 AliCDBManager *cdb = AliCDBManager::Instance();
819 if (fOCDBpath.Length()){
820 cdb->SetDefaultStorage(fOCDBpath.Data());
821 Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
824 cdb->SetRun(event->GetRunNumber());
827 // EMCAL from RAW OCDB
828 if (fOCDBpath.Contains("alien:"))
830 cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
831 cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
834 TString path = cdb->GetDefaultStorage()->GetBaseFolder();
838 //Get calibration parameters
841 AliCDBEntry *entry = (AliCDBEntry*)
842 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
844 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
848 AliFatal("Calibration parameters not found in CDB!");
850 //Get calibration parameters
853 AliCDBEntry *entry = (AliCDBEntry*)
854 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
856 if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
860 AliFatal("Dead map not found in CDB!");
865 //________________________________________________________________________________________
866 void AliAnalysisTaskEMCALClusterize::InitClusterization()
868 //Select clusterization/unfolding algorithm and set all the needed parameters
871 // init the unfolding afterburner
873 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
877 //First init the clusterizer
879 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
880 fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
881 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv2)
882 fClusterizer = new AliEMCALClusterizerv2(fGeom, fCalibData, fPedestalData);
883 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN){
884 fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
885 fClusterizer->SetNRowDiff(fRecParam->GetNRowDiff());
886 fClusterizer->SetNColDiff(fRecParam->GetNColDiff());
888 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
891 //Now set the parameters
892 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
893 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
894 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
895 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
896 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
897 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
898 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
899 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
900 fClusterizer->SetInputCalibrated ( kTRUE );
901 fClusterizer->SetJustClusters ( kTRUE );
902 //In case of unfolding after clusterization is requested, set the corresponding parameters
903 if(fRecParam->GetUnfold()){
905 for (i = 0; i < 8; i++) {
906 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
907 }//end of loop over parameters
908 for (i = 0; i < 3; i++) {
909 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
910 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
911 }//end of loop over parameters
913 fClusterizer->InitClusterUnfolding();
918 //________________________________________________________________________________________
919 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
921 // Restore clusters from recPoints
922 // Cluster energy, global position, cells and their amplitude fractions are restored
924 for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
926 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
928 const Int_t ncells = recPoint->GetMultiplicity();
929 Int_t ncellsTrue = 0;
931 if(recPoint->GetEnergy() < fRecParam->GetClusteringThreshold()) continue;
933 // cells and their amplitude fractions
934 UShort_t absIds[ncells];
935 Double32_t ratios[ncells];
937 //For later check embedding
938 AliAODInputHandler* aodIH = dynamic_cast<AliAODInputHandler*>((AliAnalysisManager::GetAnalysisManager())->GetInputEventHandler());
940 Float_t clusterE = 0;
941 for (Int_t c = 0; c < ncells; c++) {
942 AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
944 absIds[ncellsTrue] = digit->GetId();
945 ratios[ncellsTrue] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
947 // In case of unfolding, remove digits with unfolded energy too low
949 if (recPoint->GetEnergiesList()[c] < fSelectCellMinE || ratios[ncellsTrue] < fSelectCellMinFrac) {
951 if(DebugLevel() > 1) {
952 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Too small energy in cell of cluster: cluster cell %f, digit %f\n",
953 recPoint->GetEnergiesList()[c],digit->GetAmplitude());
961 //Recalculate cluster energy and number of cluster cells in case any of the cells was rejected
963 clusterE +=recPoint->GetEnergiesList()[c];
965 // In case of embedding, fill ratio with amount of signal,
966 if (aodIH && aodIH->GetMergeEvents()) {
968 //AliVCaloCells* inEMCALCells = InputEvent()->GetEMCALCells();
969 AliVCaloCells* meEMCALCells = aodIH->GetEventToMerge()->GetEMCALCells();
970 AliVCaloCells* ouEMCALCells = AODEvent()->GetEMCALCells();
972 Float_t sigAmplitude = meEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
973 //Float_t bkgAmplitude = inEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
974 Float_t sumAmplitude = ouEMCALCells->GetCellAmplitude(absIds[ncellsTrue]);
975 //printf("\t AbsID %d, amplitude : bkg %f, sigAmplitude %f, summed %f - %f\n",absIds[ncellsTrue], bkgAmplitude, sigAmplitude, sumAmplitude, digit->GetAmplitude());
977 if(sumAmplitude > 0) ratios[ncellsTrue] = sigAmplitude/sumAmplitude;
978 //printf("\t \t ratio %f\n",ratios[ncellsTrue]);
982 }// cluster cell loop
984 if (ncellsTrue < 1) {
985 if (DebugLevel() > 1)
986 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Skipping cluster with no cells avobe threshold E = %f, ncells %d\n",
987 recPoint->GetEnergy(), ncells);
991 //if(ncellsTrue != ncells) printf("Old E %f, ncells %d; New E %f, ncells %d\n",recPoint->GetEnergy(),ncells,clusterE,ncellsTrue);
993 if(clusterE < fRecParam->GetClusteringThreshold()) {
995 printf("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters() - Remove cluster with energy below seed threshold %f\n",clusterE);
1002 // calculate new cluster position
1003 recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
1004 recPoint->GetGlobalPosition(gpos);
1007 // create a new cluster
1008 //AliAODCaloCluster *clus = new AliAODCaloCluster();
1009 (*clusArray)[j] = new AliAODCaloCluster() ;
1010 AliAODCaloCluster *clus = dynamic_cast<AliAODCaloCluster *>( clusArray->At(j) ) ;
1012 clus->SetType(AliVCluster::kEMCALClusterv1);
1013 clus->SetE(clusterE);//recPoint->GetEnergy());
1014 clus->SetPosition(g);
1015 clus->SetNCells(ncellsTrue);
1016 clus->SetCellsAbsId(absIds);
1017 clus->SetCellsAmplitudeFraction(ratios);
1018 clus->SetChi2(-1); //not yet implemented
1019 clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
1020 clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
1021 clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
1023 if(ncells == ncellsTrue){
1024 Float_t elipAxis[2];
1025 recPoint->GetElipsAxis(elipAxis);
1026 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
1027 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
1028 clus->SetDispersion(recPoint->GetDispersion());
1030 else if(fSelectCell){
1031 // In case some cells rejected, in unfolding case, recalculate
1032 // shower shape parameters and position
1033 AliVCaloCells* cells = 0x0;
1034 if (aodIH && aodIH->GetMergeEvents()) cells = AODEvent() ->GetEMCALCells();
1035 else cells = InputEvent()->GetEMCALCells();
1036 fRecoUtils->RecalculateClusterShowerShapeParameters(fGeom,cells,clus);
1037 fRecoUtils->RecalculateClusterPID(clus);
1038 fRecoUtils->RecalculateClusterPosition(fGeom,cells,clus);
1040 // Float_t elipAxis[2];
1041 // recPoint->GetElipsAxis(elipAxis);
1042 // if(ncellsTrue > 2)
1043 // printf("SS, old: l0 %f, l1 %f, D %f; New l0 %f, l1 %f, D %f\n",
1044 // elipAxis[0]*elipAxis[0],elipAxis[1]*elipAxis[1], recPoint->GetDispersion(),
1045 // clus->GetM02(),clus->GetM20(),clus->GetDispersion());
1049 Int_t parentMult = 0;
1050 Int_t *parentList = recPoint->GetParents(parentMult);
1051 clus->SetLabel(parentList, parentMult);
1053 // if(parentMult)printf("RecToESD: Labels %d",parentMult);
1054 // for (Int_t iii = 0; iii < parentMult; iii++) {
1055 // printf("\t label %d\n",parentList[iii]);
1058 //Write the second major contributor to each MC cluster.
1060 for ( Int_t iLoopCell = 0 ; iLoopCell < clus->GetNCells() ; iLoopCell++ ){
1062 Int_t idCell = clus->GetCellAbsId(iLoopCell) ;
1064 iNewLabel = 1 ; //iNewLabel makes sure we don't write twice the same label.
1065 for ( UInt_t iLoopLabels = 0 ; iLoopLabels < clus->GetNLabels() ; iLoopLabels++ )
1067 if ( fCellSecondLabels[idCell] == -1 ) iNewLabel = 0; // -1 is never a good second label.
1068 if ( fCellSecondLabels[idCell] == clus->GetLabelAt(iLoopLabels) ) iNewLabel = 0;
1072 //Int_t * newLabelArray = (Int_t*)malloc((clus->GetNLabels()+1)*sizeof(Int_t)) ;
1073 Int_t * newLabelArray = new Int_t[clus->GetNLabels()+1] ;
1074 for ( UInt_t iLoopNewLabels = 0 ; iLoopNewLabels < clus->GetNLabels() ; iLoopNewLabels++ ){
1075 newLabelArray[iLoopNewLabels] = clus->GetLabelAt(iLoopNewLabels) ;
1078 newLabelArray[clus->GetNLabels()] = fCellSecondLabels[idCell] ;
1079 //if(fCellSecondLabels[idCell]>-1)printf("\t new label %d\n",fCellSecondLabels[idCell]);
1080 clus->SetLabel(newLabelArray,clus->GetNLabels()+1) ;
1081 delete [] newLabelArray;
1083 }//positive cell number