Implementing centrality bins for track finding efficiency for lead-lead collisions
[u/mrichter/AliRoot.git] / PWG4 / CaloCalib / AliAnalysisTaskEMCALClusterizeFast.cxx
CommitLineData
2f7259cf 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
16/* $Id$ */
17
18// --- Root ---
19#include <TClonesArray.h>
20#include <TGeoManager.h>
21#include <TObjArray.h>
22#include <TString.h>
23#include <TTree.h>
24#include "AliAODCaloCluster.h"
25#include "AliAODEvent.h"
26#include "AliAnalysisManager.h"
27#include "AliCDBEntry.h"
28#include "AliCDBManager.h"
29#include "AliCaloCalibPedestal.h"
30#include "AliEMCALAfterBurnerUF.h"
31#include "AliEMCALCalibData.h"
32#include "AliEMCALClusterizerNxN.h"
33#include "AliEMCALClusterizerv1.h"
34#include "AliEMCALDigit.h"
35#include "AliEMCALGeometry.h"
36#include "AliEMCALRecParam.h"
37#include "AliEMCALRecPoint.h"
38#include "AliEMCALRecoUtils.h"
39#include "AliESDEvent.h"
40#include "AliLog.h"
41
42#include "AliAnalysisTaskEMCALClusterizeFast.h"
43
44ClassImp(AliAnalysisTaskEMCALClusterizeFast)
45
46//________________________________________________________________________
47AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
48 : AliAnalysisTaskSE(name),
49 fRun(0),
50 fDigitsArr(0),
51 fClusterArr(0),
52 fRecParam(new AliEMCALRecParam),
53 fClusterizer(0),
54 fUnfolder(0),
55 fJustUnfold(kFALSE),
56 fGeomName("EMCAL_FIRSTYEARV1"),
57 fGeomMatrixSet(kFALSE),
58 fLoadGeomMatrices(kFALSE),
59 fOCDBpath(),
60 fCalibData(0),
61 fPedestalData(0),
62 fOutputAODBranch(0),
63 fOutputAODBrName(),
10d33986 64 fRecoUtils(0),
65 fLoadCalib(0),
66 fLoadPed(0)
2f7259cf 67{
68 // Constructor
69
10922c99 70 fBranchNames = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
2f7259cf 71 for(Int_t i = 0; i < 10; ++i)
72 fGeomMatrix[i] = 0;
73}
74
75//________________________________________________________________________
76AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
77{
78 // Destructor.
79
80 delete fDigitsArr;
81 delete fClusterizer;
82 delete fUnfolder;
83 delete fRecoUtils;
84}
85
86//-------------------------------------------------------------------
87void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
88{
89 // Create output objects.
90
91 if (!fOutputAODBrName.IsNull()) {
92 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
93 fOutputAODBranch->SetName(fOutputAODBrName);
94 AddAODBranch("TClonesArray", &fOutputAODBranch);
95 AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
96 }
97}
98
99//________________________________________________________________________
100void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
101{
102 // Main loop, called for each event
103
104 // remove the contents of output list set in the previous event
105 if (fOutputAODBranch)
106 fOutputAODBranch->Clear("C");
107
108 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
109 AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
110
111 if (!esdevent&&!aodevent) {
112 Error("UserExec","Event not available");
113 return;
114 }
115
116 LoadBranches();
117
118 Init();
119
120 if (fJustUnfold) {
121 AliWarning("Unfolding not implemented");
122 } else {
123 if (esdevent)
124 FillDigitsArray(esdevent);
125 else
126 FillDigitsArray(aodevent);
127 fClusterizer->Digits2Clusters("");
128 if (esdevent &&fRecoUtils)
129 fRecoUtils->FindMatches(esdevent,fClusterArr);
130 if (fOutputAODBranch)
131 RecPoints2Clusters();
132 }
133}
134
135//________________________________________________________________________
136void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliAODEvent *event)
137{
138 // Fill digits from cells.
139
140 fDigitsArr->Clear("C");
141 AliAODCaloCells *cells = event->GetEMCALCells();
142 Int_t ncells = cells->GetNumberOfCells();
143 if (ncells>fDigitsArr->GetSize())
144 fDigitsArr->Expand(2*ncells);
145 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
146 Double_t cellAmplitude=0, cellTime=0;
147 Short_t cellNumber=0;
148 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
149 break;
150 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
151 digit->SetId(cellNumber);
152 digit->SetAmplitude(cellAmplitude);
153 digit->SetTime(cellTime);
154 digit->SetTimeR(cellTime);
155 digit->SetIndexInList(idigit);
156 digit->SetType(AliEMCALDigit::kHG);
157 idigit++;
158 }
159}
160
161//________________________________________________________________________
162void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray(AliESDEvent *event)
163{
164 // Fill digits from cells.
165
166 fDigitsArr->Clear("C");
167 AliESDCaloCells *cells = event->GetEMCALCells();
168 Int_t ncells = cells->GetNumberOfCells();
169 if (ncells>fDigitsArr->GetSize())
170 fDigitsArr->Expand(2*ncells);
171 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
172 Double_t cellAmplitude=0, cellTime=0;
173 Short_t cellNumber=0;
174 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
175 break;
176 AliEMCALDigit *digit = (AliEMCALDigit*) fDigitsArr->New(idigit);
177 digit->SetId(cellNumber);
178 digit->SetAmplitude(cellAmplitude);
179 digit->SetTime(cellTime);
180 digit->SetTimeR(cellTime);
181 digit->SetIndexInList(idigit);
182 digit->SetType(AliEMCALDigit::kHG);
183 idigit++;
184 }
185}
186
187//________________________________________________________________________________________
188void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters()
189{
190 // Cluster energy, global position, cells and their amplitude fractions are restored.
191
192 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
193
194 Int_t Ncls = fClusterArr->GetEntriesFast();
195 for(Int_t i=0, nout=0; i < Ncls; ++i) {
196 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
197 Int_t ncells_true = 0;
198 const Int_t ncells = recpoint->GetMultiplicity();
199 UShort_t absIds[ncells];
200 Double32_t ratios[ncells];
201 Int_t *dlist = recpoint->GetDigitsList();
202 Float_t *elist = recpoint->GetEnergiesList();
203
204 for (Int_t c = 0; c < ncells; ++c) {
205 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
206 absIds[ncells_true] = digit->GetId();
207 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
208 if (ratios[ncells_true] > 0.001)
209 ++ncells_true;
210 }
211
212 if (ncells_true < 1) {
213 AliWarning("Skipping cluster with no cells");
214 continue;
215 }
216
217 // calculate new cluster position
218 TVector3 gpos;
219 Float_t g[3];
220
221 recpoint->EvalGlobalPosition(fRecParam->GetW0(), fDigitsArr);
222 recpoint->GetGlobalPosition(gpos);
223 gpos.GetXYZ(g);
224
225 AliAODCaloCluster *clus = static_cast<AliAODCaloCluster*>(fOutputAODBranch->New(nout++));
226 clus->SetType(AliVCluster::kEMCALClusterv1);
227 clus->SetE(recpoint->GetEnergy());
228 clus->SetPosition(g);
229 clus->SetNCells(ncells_true);
230 clus->SetCellsAbsId(absIds);
231 clus->SetCellsAmplitudeFraction(ratios);
232 clus->SetDispersion(recpoint->GetDispersion());
233 clus->SetChi2(-1); //not yet implemented
234 clus->SetTOF(recpoint->GetTime()) ; //time-of-flight
235 clus->SetNExMax(recpoint->GetNExMax()); //number of local maxima
236 Float_t elipAxis[2];
237 recpoint->GetElipsAxis(elipAxis);
238 clus->SetM02(elipAxis[0]*elipAxis[0]) ;
239 clus->SetM20(elipAxis[1]*elipAxis[1]) ;
240 clus->SetDistToBadChannel(recpoint->GetDistanceToBadTower());
241 if (esdevent && fRecoUtils) {
242 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
243 if(trackIndex >= 0) {
244 clus->AddTrackMatched(esdevent->GetTrack(trackIndex));
245 if(DebugLevel() > 1)
246 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
247 }
248 }
249 }
250}
251
252//________________________________________________________________________________________
253void AliAnalysisTaskEMCALClusterizeFast::Init()
254{
255 //Select clusterization/unfolding algorithm and set all the needed parameters
256
257 AliVEvent * event = InputEvent();
258 if (!event) {
259 AliWarning("Event not available!!!");
260 return;
261 }
262
263 if (event->GetRunNumber()==fRun)
264 return;
265 fRun = event->GetRunNumber();
266
267 if (fJustUnfold){
268 // init the unfolding afterburner
269 delete fUnfolder;
270 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
271 return;
272 }
273
274 AliCDBManager *cdb = AliCDBManager::Instance();
275 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
276 cdb->SetDefaultStorage(fOCDBpath);
277 if (fRun!=cdb->GetRun())
278 cdb->SetRun(fRun);
279
280 AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
281 if (!geometry) {
282 AliFatal("Geometry not available!!!");
283 return;
284 }
285
286 if (!fGeomMatrixSet) {
287 if (fLoadGeomMatrices) {
288 for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
289 if(fGeomMatrix[mod]){
290 if(DebugLevel() > 2)
291 fGeomMatrix[mod]->Print();
292 geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
293 }
294 }
295 } else {
296 for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
297 if(event->GetEMCALMatrix(mod)) {
298 if(DebugLevel() > 2)
299 event->GetEMCALMatrix(mod)->Print();
300 geometry->SetMisalMatrix(event->GetEMCALMatrix(mod),mod);
301 }
302 }
303 }
304 fGeomMatrixSet=kTRUE;
305 }
306
307 // setup digit array if needed
308 if (!fDigitsArr) {
309 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
310 fDigitsArr->SetOwner(1);
311 }
312
313 // then setup clusterizer
314 delete fClusterizer;
315 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
316 fClusterizer = new AliEMCALClusterizerv1(geometry);
317 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
318 fClusterizer = new AliEMCALClusterizerNxN(geometry);
319 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
320 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
321 clusterizer->SetNRowDiff(2);
322 clusterizer->SetNColDiff(2);
323 fClusterizer = clusterizer;
324 } else{
325 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
326 }
327 fClusterizer->InitParameters(fRecParam);
10d33986 328 if (!fCalibData&&fLoadCalib) {
2f7259cf 329 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
330 if (entry)
331 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
332 if (!fCalibData)
333 AliFatal("Calibration parameters not found in CDB!");
334 }
10d33986 335 if (!fPedestalData&&fLoadPed) {
2f7259cf 336 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
337 if (entry)
338 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
339 }
10d33986 340 if (fCalibData) {
341 fClusterizer->SetInputCalibrated(kFALSE);
342 fClusterizer->SetCalibrationParameters(fCalibData);
343 fClusterizer->SetCaloCalibPedestal(fPedestalData);
344 } else {
345 fClusterizer->SetInputCalibrated(kTRUE);
346 }
347 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 348 fClusterizer->SetDigitsArr(fDigitsArr);
349 fClusterizer->SetOutput(0);
350 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
351}