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 // --- AliRoot Analysis Steering
34 #include "AliAnalysisTask.h"
35 #include "AliAnalysisManager.h"
36 #include "AliESDEvent.h"
37 #include "AliGeomManager.h"
38 #include "AliVCaloCells.h"
39 #include "AliAODCaloCluster.h"
40 #include "AliCDBManager.h"
41 #include "AliCDBStorage.h"
42 #include "AliCDBEntry.h"
44 #include "AliVEventHandler.h"
47 #include "AliEMCALRecParam.h"
48 #include "AliEMCALAfterBurnerUF.h"
49 #include "AliEMCALGeometry.h"
50 #include "AliEMCALClusterizerNxN.h"
51 #include "AliEMCALClusterizerv1.h"
52 #include "AliEMCALRecPoint.h"
53 #include "AliEMCALDigit.h"
54 #include "AliCaloCalibPedestal.h"
55 #include "AliEMCALCalibData.h"
56 #include "AliEMCALRecoUtils.h"
58 #include "AliAnalysisTaskEMCALClusterize.h"
60 ClassImp(AliAnalysisTaskEMCALClusterize)
62 //________________________________________________________________________
63 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
64 : AliAnalysisTaskSE(name)
65 , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
66 , fCalibData(0), fPedestalData(0), fOCDBpath("raw://")
67 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
68 , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
69 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kTRUE)
70 , fRun(-1), fRecoUtils(0)
74 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
75 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
76 fClusterArr = new TObjArray(100);
77 fCaloClusterArr = new TObjArray(100);
78 fRecParam = new AliEMCALRecParam;
79 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
80 fRecoUtils = new AliEMCALRecoUtils();
83 //________________________________________________________________________
84 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
85 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
86 , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
87 , fCalibData(0), fPedestalData(0), fOCDBpath("raw://")
88 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
89 , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
90 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kFALSE)
91 , fRun(-1), fRecoUtils(0)
94 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
95 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
96 fClusterArr = new TObjArray(100);
97 fCaloClusterArr = new TObjArray(100);
98 fRecParam = new AliEMCALRecParam;
99 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
100 fRecoUtils = new AliEMCALRecoUtils();
103 //________________________________________________________________________
104 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
109 fDigitsArr->Clear("C");
114 fClusterArr->Delete();
118 if (fCaloClusterArr){
119 fCaloClusterArr->Delete();
120 delete fCaloClusterArr;
123 if(fClusterizer) delete fClusterizer;
124 if(fUnfolder) delete fUnfolder;
125 if(fRecoUtils) delete fRecoUtils;
129 //-------------------------------------------------------------------
130 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
132 // Init geometry, create list of output clusters
134 fGeom = AliEMCALGeometry::GetInstance(fGeomName) ;
136 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
137 fOutputAODBranch->SetName(fOutputAODBranchName);
138 AddAODBranch("TClonesArray", &fOutputAODBranch);
139 Info("UserCreateOutputObjects","Create Branch: %s",fOutputAODBranchName.Data());
142 //________________________________________________________________________
143 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
146 // Called for each event
148 //Remove the contents of output list set in the previous event
149 fOutputAODBranch->Clear("C");
151 AliVEvent * event = InputEvent();
152 AliESDEvent * esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
155 Error("UserExec","Event not available");
159 //Magic line to write events to AOD file
160 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
165 //-------------------------------------------------------------------------------------
166 //Set the geometry matrix, for the first event, skip the rest
167 //-------------------------------------------------------------------------------------
169 if(fLoadGeomMatrices){
170 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
171 if(fGeomMatrix[mod]){
173 fGeomMatrix[mod]->Print();
174 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
176 fGeomMatrixSet=kTRUE;
179 else if(!gGeoManager){
180 Info("UserExec","Get geo matrices from data");
181 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
182 if(!strcmp(event->GetName(),"AliAODEvent")) {
184 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
188 Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
189 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
191 Error("UserExec","This event does not contain ESDs?");
194 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
196 esd->GetEMCALMatrix(mod)->Print();
197 if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
199 fGeomMatrixSet=kTRUE;
201 }//Load matrices from Data
204 //Get the list of cells needed for unfolding and reclustering
205 AliVCaloCells *cells= event->GetEMCALCells();
206 Int_t nClustersOrg = 0;
207 //-------------------------------------------
208 //---------Unfolding clusters----------------
209 //-------------------------------------------
212 //Fill the array with the EMCAL clusters, copy them
213 for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
215 AliVCluster *clus = event->GetCaloCluster(i);
217 //printf("Org Cluster %d, E %f\n",i, clus->E());
218 AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
219 AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
221 fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
224 fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
227 Warning("UserExec()"," - Wrong CaloCluster type?");
233 fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
238 //printf("nClustersOrg %d\n",nClustersOrg);
240 //-------------------------------------------
241 //---------- Recluster cells ----------------
242 //-------------------------------------------
245 //-------------------------------------------------------------------------------------
246 //Transform CaloCells into Digits
247 //-------------------------------------------------------------------------------------
253 Double_t cellAmplitude = 0;
254 Double_t cellTime = 0;
255 Short_t cellNumber = 0;
257 TTree *digitsTree = new TTree("digitstree","digitstree");
258 digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
260 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
262 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
269 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
271 digit->SetAmplitude(amp);
272 digit->SetTime(time);
273 digit->SetTimeR(time);
274 digit->SetIndexInList(idigit);
275 digit->SetType(AliEMCALDigit::kHG);
276 //if(Entry()==0) printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
280 //Fill the tree with digits
283 //-------------------------------------------------------------------------------------
284 //Do the clusterization
285 //-------------------------------------------------------------------------------------
286 TTree *clustersTree = new TTree("clustertree","clustertree");
288 fClusterizer->SetInput(digitsTree);
289 fClusterizer->SetOutput(clustersTree);
290 fClusterizer->Digits2Clusters("");
292 //-------------------------------------------------------------------------------------
293 //Transform the recpoints into AliVClusters
294 //-------------------------------------------------------------------------------------
296 clustersTree->SetBranchStatus("*",0); //disable all branches
297 clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
299 TBranch *branch = clustersTree->GetBranch("EMCALECARP");
300 branch->SetAddress(&fClusterArr);
303 RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
306 fClusterizer->Clear();
307 fDigitsArr ->Clear("C");
308 fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
310 clustersTree->Delete("all");
311 digitsTree ->Delete("all");
314 //Recalculate track-matching for the new clusters, only with ESDs
315 if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);
318 //-------------------------------------------------------------------------------------
319 //Put the new clusters in the AOD list
320 //-------------------------------------------------------------------------------------
322 Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntries();
323 //Info("UserExec","New clusters %d",kNumberOfCaloClusters);
324 //if(nClustersOrg!=kNumberOfCaloClusters) Info("UserExec","Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
325 for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
326 AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
327 //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
329 //Add matched track, if any, only with ESDs
331 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
333 newCluster->AddTrackMatched(event->GetTrack(trackIndex));
335 Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
338 new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
342 fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
345 //_____________________________________________________________________
346 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
348 //Access to OCDB stuff
350 AliVEvent * event = InputEvent();
353 Warning("AccessOCDB","Event not available!!!");
357 if (event->GetRunNumber()==fRun)
359 fRun = event->GetRunNumber();
361 if(DebugLevel() > 1 )
362 Info("AccessODCD()"," Begin");
364 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
367 AliCDBManager *cdb = AliCDBManager::Instance();
370 if (fOCDBpath.Length()){
371 cdb->SetDefaultStorage(fOCDBpath.Data());
372 Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
375 cdb->SetRun(event->GetRunNumber());
378 // EMCAL from RAW OCDB
379 if (fOCDBpath.Contains("alien:"))
381 cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
382 cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
385 TString path = cdb->GetDefaultStorage()->GetBaseFolder();
389 //Get calibration parameters
392 AliCDBEntry *entry = (AliCDBEntry*)
393 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
395 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
399 AliFatal("Calibration parameters not found in CDB!");
401 //Get calibration parameters
404 AliCDBEntry *entry = (AliCDBEntry*)
405 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
407 if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
411 AliFatal("Dead map not found in CDB!");
413 InitClusterization();
418 //________________________________________________________________________________________
419 void AliAnalysisTaskEMCALClusterize::InitClusterization()
421 //Select clusterization/unfolding algorithm and set all the needed parameters
424 // init the unfolding afterburner
426 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
430 //First init the clusterizer
432 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
433 fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
434 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
435 fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
436 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
437 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
438 clusterizer->SetNRowDiff(2);
439 clusterizer->SetNColDiff(2);
440 fClusterizer = clusterizer;
442 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
445 //Now set the parameters
446 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
447 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
448 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
449 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
450 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
451 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
452 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
453 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
454 fClusterizer->SetInputCalibrated ( kTRUE );
456 //In case of unfolding after clusterization is requested, set the corresponding parameters
457 if(fRecParam->GetUnfold()){
459 for (i = 0; i < 8; i++) {
460 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
461 }//end of loop over parameters
462 for (i = 0; i < 3; i++) {
463 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
464 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
465 }//end of loop over parameters
467 fClusterizer->InitClusterUnfolding();
471 //________________________________________________________________________________________
472 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
474 // Restore clusters from recPoints
475 // Cluster energy, global position, cells and their amplitude fractions are restored
477 for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
479 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
481 const Int_t ncells = recPoint->GetMultiplicity();
482 Int_t ncells_true = 0;
484 // cells and their amplitude fractions
485 UShort_t absIds[ncells];
486 Double32_t ratios[ncells];
488 for (Int_t c = 0; c < ncells; c++) {
489 AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
491 absIds[ncells_true] = digit->GetId();
492 ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
494 if (ratios[ncells_true] > 0.001) ncells_true++;
497 if (ncells_true < 1) {
498 Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
505 // calculate new cluster position
506 recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
507 recPoint->GetGlobalPosition(gpos);
510 // create a new cluster
511 AliAODCaloCluster *clus = new AliAODCaloCluster();
512 clus->SetType(AliVCluster::kEMCALClusterv1);
513 clus->SetE(recPoint->GetEnergy());
514 clus->SetPosition(g);
515 clus->SetNCells(ncells_true);
516 clus->SetCellsAbsId(absIds);
517 clus->SetCellsAmplitudeFraction(ratios);
518 clus->SetDispersion(recPoint->GetDispersion());
519 clus->SetChi2(-1); //not yet implemented
520 clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
521 clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
523 recPoint->GetElipsAxis(elipAxis);
524 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
525 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
526 clus->SetDistToBadChannel(recPoint->GetDistanceToBadTower());
527 clusArray->Add(clus);