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