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"
35 // --- AliRoot Analysis Steering
36 #include "AliAnalysisTask.h"
37 #include "AliAnalysisManager.h"
38 #include "AliESDEvent.h"
39 #include "AliGeomManager.h"
40 #include "AliVCaloCells.h"
41 #include "AliAODCaloCluster.h"
42 #include "AliCDBManager.h"
43 #include "AliCDBStorage.h"
44 #include "AliCDBEntry.h"
46 #include "AliVEventHandler.h"
49 #include "AliEMCALRecParam.h"
50 #include "AliEMCALAfterBurnerUF.h"
51 #include "AliEMCALGeometry.h"
52 #include "AliEMCALClusterizerNxN.h"
53 #include "AliEMCALClusterizerv1.h"
54 #include "AliEMCALRecPoint.h"
55 #include "AliEMCALDigit.h"
56 #include "AliCaloCalibPedestal.h"
57 #include "AliEMCALCalibData.h"
58 #include "AliEMCALRecoUtils.h"
60 #include "AliAnalysisTaskEMCALClusterize.h"
62 ClassImp(AliAnalysisTaskEMCALClusterize)
64 //________________________________________________________________________
65 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize(const char *name)
66 : AliAnalysisTaskSE(name)
67 , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
68 , fCalibData(0), fPedestalData(0), fOCDBpath("raw://")
69 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
70 , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
71 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kTRUE)
72 , fRun(-1), fRecoUtils(0), fConfigName("")
76 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
77 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
78 fClusterArr = new TObjArray(100);
79 fCaloClusterArr = new TObjArray(100);
80 fRecParam = new AliEMCALRecParam;
81 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
82 fRecoUtils = new AliEMCALRecoUtils();
85 //________________________________________________________________________
86 AliAnalysisTaskEMCALClusterize::AliAnalysisTaskEMCALClusterize()
87 : AliAnalysisTaskSE("DefaultAnalysis_AliAnalysisTaskEMCALClusterize")
88 , fGeom(0), fGeomName("EMCAL_FIRSTYEARV1"), fGeomMatrixSet(kFALSE), fLoadGeomMatrices(kFALSE)
89 , fCalibData(0), fPedestalData(0), fOCDBpath("raw://")
90 , fDigitsArr(0), fClusterArr(0), fCaloClusterArr(0)
91 , fRecParam(0), fClusterizer(0), fUnfolder(0), fJustUnfold(kFALSE)
92 , fOutputAODBranch(0), fOutputAODBranchName("newEMCALClusters"), fFillAODFile(kFALSE)
93 , fRun(-1), fRecoUtils(0), fConfigName("")
96 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = 0 ;
97 fDigitsArr = new TClonesArray("AliEMCALDigit",200);
98 fClusterArr = new TObjArray(100);
99 fCaloClusterArr = new TObjArray(100);
100 fRecParam = new AliEMCALRecParam;
101 fBranchNames = "ESD:AliESDHeader.,EMCALCells.";
102 fRecoUtils = new AliEMCALRecoUtils();
106 //________________________________________________________________________
107 AliAnalysisTaskEMCALClusterize::~AliAnalysisTaskEMCALClusterize()
112 fDigitsArr->Clear("C");
117 fClusterArr->Delete();
121 if (fCaloClusterArr){
122 fCaloClusterArr->Delete();
123 delete fCaloClusterArr;
126 if(fClusterizer) delete fClusterizer;
127 if(fUnfolder) delete fUnfolder;
128 if(fRecoUtils) delete fRecoUtils;
132 //-------------------------------------------------------------------
133 void AliAnalysisTaskEMCALClusterize::Init()
135 //Init analysis with configuration macro if available
137 if(gROOT->LoadMacro(fConfigName) >=0){
138 printf("Configure analysis with %s\n",fConfigName.Data());
139 AliAnalysisTaskEMCALClusterize *clus = (AliAnalysisTaskEMCALClusterize*)gInterpreter->ProcessLine("ConfigEMCALClusterize()");
140 fGeomName = clus->fGeomName;
141 fLoadGeomMatrices = clus->fLoadGeomMatrices;
142 fOCDBpath = clus->fOCDBpath;
143 fRecParam = clus->fRecParam;
144 fJustUnfold = clus->fJustUnfold;
145 fFillAODFile = clus->fFillAODFile;
146 fRecoUtils = clus->fRecoUtils;
147 fConfigName = clus->fConfigName;
148 fOutputAODBranchName = clus->fOutputAODBranchName;
149 for(Int_t i = 0; i < 10; i++) fGeomMatrix[i] = clus->fGeomMatrix[i] ;
155 //-------------------------------------------------------------------
156 void AliAnalysisTaskEMCALClusterize::UserCreateOutputObjects()
158 // Init geometry, create list of output clusters
160 fGeom = AliEMCALGeometry::GetInstance(fGeomName) ;
162 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
163 fOutputAODBranch->SetName(fOutputAODBranchName);
164 AddAODBranch("TClonesArray", &fOutputAODBranch);
168 //________________________________________________________________________
169 void AliAnalysisTaskEMCALClusterize::UserExec(Option_t *)
172 // Called for each event
174 //printf("--- Event %d -- \n",(Int_t)Entry());
175 //Remove the contents of output list set in the previous event
176 fOutputAODBranch->Clear("C");
178 AliVEvent * event = InputEvent();
179 AliESDEvent * esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
182 Error("UserExec","Event not available");
186 //Magic line to write events to AOD file
187 AliAnalysisManager::GetAnalysisManager()->GetOutputEventHandler()->SetFillAOD(fFillAODFile);
192 //-------------------------------------------------------------------------------------
193 //Set the geometry matrix, for the first event, skip the rest
194 //-------------------------------------------------------------------------------------
196 if(fLoadGeomMatrices){
197 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
198 if(fGeomMatrix[mod]){
200 fGeomMatrix[mod]->Print();
201 fGeom->SetMisalMatrix(fGeomMatrix[mod],mod) ;
203 fGeomMatrixSet=kTRUE;
206 else if(!gGeoManager){
207 Info("UserExec","Get geo matrices from data");
208 //Still not implemented in AOD, just a workaround to be able to work at least with ESDs
209 if(!strcmp(event->GetName(),"AliAODEvent")) {
211 Warning("UserExec","Use ideal geometry, values geometry matrix not kept in AODs.");
215 Info("UserExec","AliAnalysisTaskEMCALClusterize Load Misaligned matrices.");
216 AliESDEvent* esd = dynamic_cast<AliESDEvent*>(event) ;
218 Error("UserExec","This event does not contain ESDs?");
221 for(Int_t mod=0; mod < (fGeom->GetEMCGeometry())->GetNumberOfSuperModules(); mod++){
223 esd->GetEMCALMatrix(mod)->Print();
224 if(esd->GetEMCALMatrix(mod)) fGeom->SetMisalMatrix(esd->GetEMCALMatrix(mod),mod) ;
226 fGeomMatrixSet=kTRUE;
228 }//Load matrices from Data
230 //Recover time dependent corrections, put then in recalibration histograms. Do it once
231 fRecoUtils->SetTimeDependentCorrections(InputEvent()->GetRunNumber());
235 //Get the list of cells needed for unfolding and reclustering
236 AliVCaloCells *cells= event->GetEMCALCells();
237 Int_t nClustersOrg = 0;
238 //-------------------------------------------
239 //---------Unfolding clusters----------------
240 //-------------------------------------------
243 //Fill the array with the EMCAL clusters, copy them
244 for (Int_t i = 0; i < event->GetNumberOfCaloClusters(); i++)
246 AliVCluster *clus = event->GetCaloCluster(i);
248 //printf("Org Cluster %d, E %f\n",i, clus->E());
250 //recalibrate/remove bad channels/etc if requested
251 if(fRecoUtils->ClusterContainsBadChannel(fGeom,clus->GetCellsAbsId(), clus->GetNCells())){
252 //printf("Remove cluster\n");
256 if(fRecoUtils->IsRecalibrationOn()){
257 //printf("Energy before %f ",clus->E());
258 fRecoUtils->RecalibrateClusterEnergy(fGeom, clus, cells);
259 //printf("; after %f\n",clus->E());
261 //Cast to ESD or AOD, needed to create the cluster array
262 AliESDCaloCluster * esdCluster = dynamic_cast<AliESDCaloCluster*> (clus);
263 AliAODCaloCluster * aodCluster = dynamic_cast<AliAODCaloCluster*> (clus);
265 fCaloClusterArr->Add( new AliESDCaloCluster(*esdCluster) );
268 fCaloClusterArr->Add( new AliAODCaloCluster(*aodCluster) );
271 Warning("UserExec()"," - Wrong CaloCluster type?");
277 fUnfolder->UnfoldClusters(fCaloClusterArr, cells);
283 //printf("nClustersOrg %d\n",nClustersOrg);
285 //-------------------------------------------
286 //---------- Recluster cells ----------------
287 //-------------------------------------------
290 //-------------------------------------------------------------------------------------
291 //Transform CaloCells into Digits
292 //-------------------------------------------------------------------------------------
298 Double_t cellAmplitude = 0;
299 Double_t cellTime = 0;
300 Short_t cellNumber = 0;
302 TTree *digitsTree = new TTree("digitstree","digitstree");
303 digitsTree->Branch("EMCAL","TClonesArray", &fDigitsArr, 32000);
305 for (Int_t icell = 0; icell < cells->GetNumberOfCells(); icell++)
307 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
314 Int_t imod = -1, iphi =-1, ieta=-1,iTower = -1, iIphi = -1, iIeta = -1;
315 fGeom->GetCellIndex(id,imod,iTower,iIphi,iIeta);
316 fGeom->GetCellPhiEtaIndexInSModule(imod,iTower,iIphi, iIeta,iphi,ieta);
318 //Do not include bad channels found in analysis?
319 if( fRecoUtils->IsBadChannelsRemovalSwitchedOn() &&
320 fRecoUtils->GetEMCALChannelStatus(imod, ieta, iphi)){
321 //printf("Remove channel %d\n",id);
326 if(fRecoUtils->IsRecalibrationOn()){
327 //printf("CalibFactor %f times %f for id %d\n",fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi),amp,id);
328 amp *=fRecoUtils->GetEMCALChannelRecalibrationFactor(imod,ieta,iphi);
331 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
333 digit->SetAmplitude(amp);
334 digit->SetTime(time);
335 digit->SetTimeR(time);
336 digit->SetIndexInList(idigit);
337 digit->SetType(AliEMCALDigit::kHG);
339 //printf("Digit: Id %d, amp %f, time %e, index %d\n",id, amp,time,idigit);
343 //Fill the tree with digits
346 //-------------------------------------------------------------------------------------
347 //Do the clusterization
348 //-------------------------------------------------------------------------------------
349 TTree *clustersTree = new TTree("clustertree","clustertree");
351 fClusterizer->SetInput(digitsTree);
352 fClusterizer->SetOutput(clustersTree);
353 fClusterizer->Digits2Clusters("");
355 //-------------------------------------------------------------------------------------
356 //Transform the recpoints into AliVClusters
357 //-------------------------------------------------------------------------------------
359 clustersTree->SetBranchStatus("*",0); //disable all branches
360 clustersTree->SetBranchStatus("EMCALECARP",1); //Enable only the branch we need
362 TBranch *branch = clustersTree->GetBranch("EMCALECARP");
363 branch->SetAddress(&fClusterArr);
366 RecPoints2Clusters(fDigitsArr, fClusterArr, fCaloClusterArr);
369 fClusterizer->Clear();
370 fDigitsArr ->Clear("C");
371 fClusterArr ->Delete(); // Do not Clear(), it leaks, why?
373 clustersTree->Delete("all");
374 digitsTree ->Delete("all");
377 //Recalculate track-matching for the new clusters, only with ESDs
378 if(esdevent)fRecoUtils->FindMatches(esdevent,fCaloClusterArr);
381 //-------------------------------------------------------------------------------------
382 //Put the new clusters in the AOD list
383 //-------------------------------------------------------------------------------------
385 Int_t kNumberOfCaloClusters = fCaloClusterArr->GetEntries();
386 //Info("UserExec","New clusters %d",kNumberOfCaloClusters);
387 //if(nClustersOrg!=kNumberOfCaloClusters) Info("UserExec","Different number: Org %d, New %d\n",nClustersOrg,kNumberOfCaloClusters);
388 for(Int_t i = 0; i < kNumberOfCaloClusters; i++){
389 AliAODCaloCluster *newCluster = (AliAODCaloCluster *) fCaloClusterArr->At(i);
390 //if(Entry()==0) Info("UserExec","newCluster E %f\n", newCluster->E());
392 //Add matched track, if any, only with ESDs
394 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
396 newCluster->AddTrackMatched(event->GetTrack(trackIndex));
398 Info("UserExec","Matched Track index %d to new cluster %d \n",trackIndex,i);
402 //In case of new bad channels, recalculate distance to bad channels
403 if(fRecoUtils->IsBadChannelsRemovalSwitchedOn()){
404 //printf("DistToBad before %f ",newCluster->GetDistanceToBadChannel());
405 fRecoUtils->RecalculateClusterDistanceToBadChannel(fGeom, event->GetEMCALCells(), newCluster);
406 //printf("; after %f \n",newCluster->GetDistanceToBadChannel());
409 newCluster->SetID(i);
410 //printf("Cluster %d, energy %f\n",newCluster->GetID(),newCluster->E());
411 new((*fOutputAODBranch)[i]) AliAODCaloCluster(*newCluster);
415 fCaloClusterArr->Delete(); // Do not Clear(), it leaks, why?
418 //_____________________________________________________________________
419 Bool_t AliAnalysisTaskEMCALClusterize::AccessOCDB()
421 //Access to OCDB stuff
423 AliVEvent * event = InputEvent();
426 Warning("AccessOCDB","Event not available!!!");
430 if (event->GetRunNumber()==fRun)
432 fRun = event->GetRunNumber();
434 if(DebugLevel() > 1 )
435 Info("AccessODCD()"," Begin");
437 fGeom = AliEMCALGeometry::GetInstance(fGeomName);
440 AliCDBManager *cdb = AliCDBManager::Instance();
443 if (fOCDBpath.Length()){
444 cdb->SetDefaultStorage(fOCDBpath.Data());
445 Info("AccessOCDB","Default storage %s",fOCDBpath.Data());
448 cdb->SetRun(event->GetRunNumber());
451 // EMCAL from RAW OCDB
452 if (fOCDBpath.Contains("alien:"))
454 cdb->SetSpecificStorage("EMCAL/Calib/Data","alien://Folder=/alice/data/2010/OCDB");
455 cdb->SetSpecificStorage("EMCAL/Calib/Pedestals","alien://Folder=/alice/data/2010/OCDB");
458 TString path = cdb->GetDefaultStorage()->GetBaseFolder();
462 //Get calibration parameters
465 AliCDBEntry *entry = (AliCDBEntry*)
466 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
468 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
472 AliFatal("Calibration parameters not found in CDB!");
474 //Get calibration parameters
477 AliCDBEntry *entry = (AliCDBEntry*)
478 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
480 if (entry) fPedestalData = (AliCaloCalibPedestal*) entry->GetObject();
484 AliFatal("Dead map not found in CDB!");
486 InitClusterization();
491 //________________________________________________________________________________________
492 void AliAnalysisTaskEMCALClusterize::InitClusterization()
494 //Select clusterization/unfolding algorithm and set all the needed parameters
497 // init the unfolding afterburner
499 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
503 //First init the clusterizer
505 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
506 fClusterizer = new AliEMCALClusterizerv1 (fGeom, fCalibData, fPedestalData);
507 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
508 fClusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
509 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
510 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(fGeom, fCalibData, fPedestalData);
511 clusterizer->SetNRowDiff(2);
512 clusterizer->SetNColDiff(2);
513 fClusterizer = clusterizer;
515 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
518 //Now set the parameters
519 fClusterizer->SetECAClusteringThreshold( fRecParam->GetClusteringThreshold() );
520 fClusterizer->SetECALogWeight ( fRecParam->GetW0() );
521 fClusterizer->SetMinECut ( fRecParam->GetMinECut() );
522 fClusterizer->SetUnfolding ( fRecParam->GetUnfold() );
523 fClusterizer->SetECALocalMaxCut ( fRecParam->GetLocMaxCut() );
524 fClusterizer->SetTimeCut ( fRecParam->GetTimeCut() );
525 fClusterizer->SetTimeMin ( fRecParam->GetTimeMin() );
526 fClusterizer->SetTimeMax ( fRecParam->GetTimeMax() );
527 fClusterizer->SetInputCalibrated ( kTRUE );
529 //In case of unfolding after clusterization is requested, set the corresponding parameters
530 if(fRecParam->GetUnfold()){
532 for (i = 0; i < 8; i++) {
533 fClusterizer->SetSSPars(i, fRecParam->GetSSPars(i));
534 }//end of loop over parameters
535 for (i = 0; i < 3; i++) {
536 fClusterizer->SetPar5 (i, fRecParam->GetPar5(i));
537 fClusterizer->SetPar6 (i, fRecParam->GetPar6(i));
538 }//end of loop over parameters
540 fClusterizer->InitClusterUnfolding();
544 //________________________________________________________________________________________
545 void AliAnalysisTaskEMCALClusterize::RecPoints2Clusters(TClonesArray *digitsArr, TObjArray *recPoints, TObjArray *clusArray)
547 // Restore clusters from recPoints
548 // Cluster energy, global position, cells and their amplitude fractions are restored
550 for(Int_t i = 0; i < recPoints->GetEntriesFast(); i++)
552 AliEMCALRecPoint *recPoint = (AliEMCALRecPoint*) recPoints->At(i);
554 const Int_t ncells = recPoint->GetMultiplicity();
555 Int_t ncells_true = 0;
557 // cells and their amplitude fractions
558 UShort_t absIds[ncells];
559 Double32_t ratios[ncells];
561 for (Int_t c = 0; c < ncells; c++) {
562 AliEMCALDigit *digit = (AliEMCALDigit*) digitsArr->At(recPoint->GetDigitsList()[c]);
564 absIds[ncells_true] = digit->GetId();
565 ratios[ncells_true] = recPoint->GetEnergiesList()[c]/digit->GetAmplitude();
567 if (ratios[ncells_true] > 0.001) ncells_true++;
570 if (ncells_true < 1) {
571 Warning("AliAnalysisTaskEMCALClusterize::RecPoints2Clusters", "skipping cluster with no cells");
578 // calculate new cluster position
579 recPoint->EvalGlobalPosition(fRecParam->GetW0(), digitsArr);
580 recPoint->GetGlobalPosition(gpos);
583 // create a new cluster
584 AliAODCaloCluster *clus = new AliAODCaloCluster();
585 clus->SetType(AliVCluster::kEMCALClusterv1);
586 clus->SetE(recPoint->GetEnergy());
587 clus->SetPosition(g);
588 clus->SetNCells(ncells_true);
589 clus->SetCellsAbsId(absIds);
590 clus->SetCellsAmplitudeFraction(ratios);
591 clus->SetDispersion(recPoint->GetDispersion());
592 clus->SetChi2(-1); //not yet implemented
593 clus->SetTOF(recPoint->GetTime()) ; //time-of-flight
594 clus->SetNExMax(recPoint->GetNExMax()); //number of local maxima
596 recPoint->GetElipsAxis(elipAxis);
597 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
598 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
599 clus->SetDistanceToBadChannel(recPoint->GetDistanceToBadTower());
600 clusArray->Add(clus);