Analysis classes to measure inclusive pi0 spectrum in pp 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//________________________________________________________________________
2431b20f 47AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast()
48 : AliAnalysisTaskSE(),
49 fRun(0),
50 fDigitsArr(0),
51 fClusterArr(0),
52 fRecParam(0),
53 fClusterizer(0),
54 fUnfolder(0),
55 fJustUnfold(kFALSE),
56 fGeomName(),
57 fGeomMatrixSet(kFALSE),
58 fLoadGeomMatrices(kFALSE),
59 fOCDBpath(),
60 fCalibData(0),
61 fPedestalData(0),
62 fOutputAODBranch(0),
63 fOutputAODBrName(),
64 fRecoUtils(0),
65 fLoadCalib(0),
4cf9204b 66 fLoadPed(0),
b0e354fe 67 fAttachClusters(0),
95f69c66 68 fRecalibOnly(0),
69 fSubBackground(0)
2431b20f 70{
71 // Constructor
72}
73
74//________________________________________________________________________
2f7259cf 75AliAnalysisTaskEMCALClusterizeFast::AliAnalysisTaskEMCALClusterizeFast(const char *name)
76 : AliAnalysisTaskSE(name),
77 fRun(0),
78 fDigitsArr(0),
79 fClusterArr(0),
80 fRecParam(new AliEMCALRecParam),
81 fClusterizer(0),
82 fUnfolder(0),
83 fJustUnfold(kFALSE),
84 fGeomName("EMCAL_FIRSTYEARV1"),
85 fGeomMatrixSet(kFALSE),
86 fLoadGeomMatrices(kFALSE),
87 fOCDBpath(),
88 fCalibData(0),
89 fPedestalData(0),
90 fOutputAODBranch(0),
91 fOutputAODBrName(),
10d33986 92 fRecoUtils(0),
93 fLoadCalib(0),
4cf9204b 94 fLoadPed(0),
b0e354fe 95 fAttachClusters(0),
95f69c66 96 fRecalibOnly(0),
97 fSubBackground(0)
2f7259cf 98{
99 // Constructor
100
95f69c66 101 fBranchNames = "ESD:AliESDHeader.,EMCALCells. AOD:header,emcalCells";
102 for(Int_t i = 0; i < 12; ++i)
2f7259cf 103 fGeomMatrix[i] = 0;
104}
105
106//________________________________________________________________________
107AliAnalysisTaskEMCALClusterizeFast::~AliAnalysisTaskEMCALClusterizeFast()
108{
109 // Destructor.
110
111 delete fDigitsArr;
112 delete fClusterizer;
113 delete fUnfolder;
114 delete fRecoUtils;
115}
116
117//-------------------------------------------------------------------
118void AliAnalysisTaskEMCALClusterizeFast::UserCreateOutputObjects()
119{
120 // Create output objects.
121
122 if (!fOutputAODBrName.IsNull()) {
123 fOutputAODBranch = new TClonesArray("AliAODCaloCluster", 0);
124 fOutputAODBranch->SetName(fOutputAODBrName);
125 AddAODBranch("TClonesArray", &fOutputAODBranch);
126 AliInfo(Form("Created Branch: %s",fOutputAODBrName.Data()));
127 }
128}
129
130//________________________________________________________________________
131void AliAnalysisTaskEMCALClusterizeFast::UserExec(Option_t *)
132{
133 // Main loop, called for each event
134
135 // remove the contents of output list set in the previous event
136 if (fOutputAODBranch)
137 fOutputAODBranch->Clear("C");
138
139 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
140 AliAODEvent *aodevent = dynamic_cast<AliAODEvent*>(InputEvent());
141
142 if (!esdevent&&!aodevent) {
143 Error("UserExec","Event not available");
144 return;
145 }
146
147 LoadBranches();
148
149 Init();
150
151 if (fJustUnfold) {
152 AliWarning("Unfolding not implemented");
84ccad86 153 return;
154 }
155
95f69c66 156 FillDigitsArray();
84ccad86 157
95f69c66 158 if (fRecalibOnly) {
159 UpdateCells();
160 return; // not requested to run clusterizer
84ccad86 161 }
7099ec08 162
95f69c66 163 Clusterize();
164 UpdateCells();
165 UpdateClusters();
2431b20f 166
95f69c66 167 if (fOutputAODBranch)
168 RecPoints2Clusters(fOutputAODBranch);
2431b20f 169}
170
171//________________________________________________________________________
95f69c66 172void AliAnalysisTaskEMCALClusterizeFast::Clusterize()
2431b20f 173{
95f69c66 174 // Clusterize
7099ec08 175
95f69c66 176 if (fSubBackground) {
177 fClusterizer->SetInputCalibrated(kTRUE);
178 fClusterizer->SetCalibrationParameters(0);
2f7259cf 179 }
95f69c66 180 fClusterizer->Digits2Clusters("");
181 if (fSubBackground) {
182 fClusterizer->SetInputCalibrated(kFALSE);
183 fClusterizer->SetCalibrationParameters(fCalibData);
4cf9204b 184 }
4cf9204b 185}
186
187//________________________________________________________________________
95f69c66 188void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
2f7259cf 189{
190 // Fill digits from cells.
191
192 fDigitsArr->Clear("C");
95f69c66 193 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
194 Double_t avgE = 0; // for background subtraction
2f7259cf 195 Int_t ncells = cells->GetNumberOfCells();
2f7259cf 196 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
197 Double_t cellAmplitude=0, cellTime=0;
198 Short_t cellNumber=0;
199 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
200 break;
95f69c66 201 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
2f7259cf 202 digit->SetId(cellNumber);
2f7259cf 203 digit->SetTime(cellTime);
204 digit->SetTimeR(cellTime);
205 digit->SetIndexInList(idigit);
206 digit->SetType(AliEMCALDigit::kHG);
95f69c66 207 if (fRecalibOnly||fSubBackground) {
b0e354fe 208 Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
95f69c66 209 digit->SetAmplitude(energy);
210 avgE += energy;
211 } else {
212 digit->SetAmplitude(cellAmplitude);
b0e354fe 213 }
2f7259cf 214 idigit++;
215 }
2f7259cf 216
95f69c66 217 if (fSubBackground) {
218 avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
219 Int_t ndigis = fDigitsArr->GetEntries();
220 for (Int_t i = 0; i < ndigis; ++i) {
221 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
222 Double_t energy = digit->GetAmplitude() - avgE;
223 if (energy<=0.001) {
224 //fDigitsArr->RemoveAt(i);
225 digit->SetAmplitude(0);
226 } else {
227 digit->SetAmplitude(energy);
228 }
b0e354fe 229 }
2f7259cf 230 }
95f69c66 231 fDigitsArr->Compress();
84ccad86 232 fDigitsArr->Sort();
2f7259cf 233}
234
235//________________________________________________________________________________________
95f69c66 236void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
4cf9204b 237{
238 // Cluster energy, global position, cells and their amplitude fractions are restored.
239
95f69c66 240 Bool_t esdobjects = 0;
241 if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
242 esdobjects = 1;
243
4cf9204b 244 Int_t Ncls = fClusterArr->GetEntriesFast();
245 for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
246 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
247 Int_t ncells_true = 0;
248 const Int_t ncells = recpoint->GetMultiplicity();
249 UShort_t absIds[ncells];
250 Double32_t ratios[ncells];
251 Int_t *dlist = recpoint->GetDigitsList();
252 Float_t *elist = recpoint->GetEnergiesList();
253 for (Int_t c = 0; c < ncells; ++c) {
254 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
255 absIds[ncells_true] = digit->GetId();
256 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
84ccad86 257 if (ratios[ncells_true] < 0.001)
258 continue;
259 ++ncells_true;
4cf9204b 260 }
261
262 if (ncells_true < 1) {
263 AliWarning("Skipping cluster with no cells");
264 continue;
265 }
266
267 // calculate new cluster position
268 TVector3 gpos;
4cf9204b 269 recpoint->GetGlobalPosition(gpos);
84ccad86 270 Float_t g[3];
4cf9204b 271 gpos.GetXYZ(g);
272
95f69c66 273 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
4cf9204b 274 c->SetType(AliVCluster::kEMCALClusterv1);
275 c->SetE(recpoint->GetEnergy());
276 c->SetPosition(g);
277 c->SetNCells(ncells_true);
4cf9204b 278 c->SetDispersion(recpoint->GetDispersion());
84ccad86 279 c->SetEmcCpvDistance(-1); //not yet implemented
4cf9204b 280 c->SetChi2(-1); //not yet implemented
281 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
282 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
283 Float_t elipAxis[2];
284 recpoint->GetElipsAxis(elipAxis);
285 c->SetM02(elipAxis[0]*elipAxis[0]) ;
286 c->SetM20(elipAxis[1]*elipAxis[1]) ;
17cff35e 287 c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
95f69c66 288 if (esdobjects) {
289 AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
290 cesd->SetCellsAbsId(absIds);
291 cesd->SetCellsAmplitudeFraction(ratios);
292 } else {
293 AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
294 caod->SetCellsAbsId(absIds);
295 caod->SetCellsAmplitudeFraction(ratios);
296 }
84ccad86 297 }
298
299 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
300 if (!esdevent)
301 return;
302 if (!fRecoUtils)
303 return;
304
95f69c66 305 AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
84ccad86 306 fRecoUtils->FindMatches(esdevent,clus);
307
95f69c66 308 if (!esdobjects) {
309 Int_t Nclus = clus->GetEntries();
310 for(Int_t i=0; i < Nclus; ++i) {
311 AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
312 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
313 if(trackIndex >= 0) {
314 Float_t dR, dZ;
315 fRecoUtils->GetMatchedResiduals(i,dR, dZ);
316 c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
317 c->SetTrackDistance(dR,dZ); // not implemented
318 c->SetEmcCpvDistance(dR);
319 c->SetChi2(dZ);
320 if(DebugLevel() > 1)
321 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
322 }
2f7259cf 323 }
95f69c66 324 } else {
84ccad86 325 Int_t Nclus = clus->GetEntries();
326 for(Int_t i=0; i < Nclus; ++i) {
327 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
2f7259cf 328 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
329 if(trackIndex >= 0) {
84ccad86 330 Float_t dR, dZ;
331 fRecoUtils->GetMatchedResiduals(i,dR, dZ);
332 c->SetTrackDistance(dR,dZ);
333 c->SetEmcCpvDistance(dR); //to be consistent with AODs
334 c->SetChi2(dZ); //to be consistent with AODs
4cf9204b 335 TArrayI tm(1,&trackIndex);
336 c->AddTracksMatched(tm);
2f7259cf 337 if(DebugLevel() > 1)
338 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
339 }
340 }
341 }
342}
343
95f69c66 344//________________________________________________________________________
345void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
346{
347 // Update cells in case re-calibration was done.
348
349 if (!fCalibData&&!fSubBackground)
350 return;
351
352 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
353 Int_t ncells = cells->GetNumberOfCells();
354 Int_t ndigis = fDigitsArr->GetEntries();
355
356 if (ncells!=ndigis) {
357 cells->DeleteContainer();
358 cells->CreateContainer(ndigis);
359 }
360 for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
361 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
362 Double_t cellAmplitude = digit->GetCalibAmp();
363 Short_t cellNumber = digit->GetId();
364 Double_t cellTime = digit->GetTime();
365 cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
366 }
367}
368
369//________________________________________________________________________
370void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
371{
372 // Update cells in case re-calibration was done.
373
374 if (!fAttachClusters)
375 return;
376
377 TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
378 if (!clus)
379 clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
380 if(!clus)
381 return;
382
383 Int_t nents = clus->GetEntries();
384 for (Int_t i=0;i<nents;++i) {
385 AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
386 if (!c)
387 continue;
388 if (c->IsEMCAL())
389 clus->RemoveAt(i);
390 }
391 clus->Compress();
392 RecPoints2Clusters(clus);
393}
394
2f7259cf 395//________________________________________________________________________________________
396void AliAnalysisTaskEMCALClusterizeFast::Init()
397{
398 //Select clusterization/unfolding algorithm and set all the needed parameters
399
400 AliVEvent * event = InputEvent();
401 if (!event) {
402 AliWarning("Event not available!!!");
403 return;
404 }
405
406 if (event->GetRunNumber()==fRun)
407 return;
408 fRun = event->GetRunNumber();
409
410 if (fJustUnfold){
411 // init the unfolding afterburner
412 delete fUnfolder;
413 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
414 return;
415 }
416
417 AliCDBManager *cdb = AliCDBManager::Instance();
418 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
419 cdb->SetDefaultStorage(fOCDBpath);
420 if (fRun!=cdb->GetRun())
421 cdb->SetRun(fRun);
422
423 AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
424 if (!geometry) {
425 AliFatal("Geometry not available!!!");
426 return;
427 }
428
429 if (!fGeomMatrixSet) {
430 if (fLoadGeomMatrices) {
431 for(Int_t mod=0; mod < (geometry->GetEMCGeometry())->GetNumberOfSuperModules(); ++mod) {
432 if(fGeomMatrix[mod]){
433 if(DebugLevel() > 2)
434 fGeomMatrix[mod]->Print();
435 geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
436 }
437 }
29b496d5 438 } else { // get matrix from file (work around bug in aliroot)
2f7259cf 439 for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5 440 const TGeoHMatrix *gm = 0;
441 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event->GetHeader());
442 if (esdevent) {
443 gm = esdevent->GetEMCALMatrix(mod);
444 } else {
445 AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
446 if (aodheader) {
447 gm = aodheader->GetEMCALMatrix(mod);
448 }
449 }
450 if (gm) {
2f7259cf 451 if(DebugLevel() > 2)
29b496d5 452 gm->Print();
453 geometry->SetMisalMatrix(gm,mod);
2f7259cf 454 }
455 }
456 }
457 fGeomMatrixSet=kTRUE;
458 }
459
460 // setup digit array if needed
461 if (!fDigitsArr) {
462 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
463 fDigitsArr->SetOwner(1);
464 }
465
466 // then setup clusterizer
467 delete fClusterizer;
468 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
469 fClusterizer = new AliEMCALClusterizerv1(geometry);
470 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
471 fClusterizer = new AliEMCALClusterizerNxN(geometry);
472 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
473 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
474 clusterizer->SetNRowDiff(2);
475 clusterizer->SetNColDiff(2);
476 fClusterizer = clusterizer;
477 } else{
478 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
479 }
480 fClusterizer->InitParameters(fRecParam);
10d33986 481 if (!fCalibData&&fLoadCalib) {
2f7259cf 482 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
483 if (entry)
484 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
485 if (!fCalibData)
486 AliFatal("Calibration parameters not found in CDB!");
487 }
10d33986 488 if (!fPedestalData&&fLoadPed) {
2f7259cf 489 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
490 if (entry)
491 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
492 }
10d33986 493 if (fCalibData) {
494 fClusterizer->SetInputCalibrated(kFALSE);
495 fClusterizer->SetCalibrationParameters(fCalibData);
496 fClusterizer->SetCaloCalibPedestal(fPedestalData);
497 } else {
498 fClusterizer->SetInputCalibrated(kTRUE);
499 }
500 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 501 fClusterizer->SetDigitsArr(fDigitsArr);
502 fClusterizer->SetOutput(0);
503 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
504}