AliTRDCalibTask.cxx .h (running on a test task over all triggers)
[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(),
51 fRun(0),
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),
79 fRun(0),
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 }
10d33986 519 if (!fCalibData&&fLoadCalib) {
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 }
10d33986 526 if (!fPedestalData&&fLoadPed) {
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}