check if calo ped is zero
[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
01bc1ce8 101 fBranchNames = "ESD:AliESDHeader.,AliESDRun.,EMCALCells. AOD:header,emcalCells";
95f69c66 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) {
4889ed8b 182 if (fCalibData) {
183 fClusterizer->SetInputCalibrated(kFALSE);
184 fClusterizer->SetCalibrationParameters(fCalibData);
185 }
4cf9204b 186 }
4cf9204b 187}
188
189//________________________________________________________________________
95f69c66 190void AliAnalysisTaskEMCALClusterizeFast::FillDigitsArray()
2f7259cf 191{
192 // Fill digits from cells.
193
194 fDigitsArr->Clear("C");
95f69c66 195 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
196 Double_t avgE = 0; // for background subtraction
2f7259cf 197 Int_t ncells = cells->GetNumberOfCells();
2f7259cf 198 for (Int_t icell = 0, idigit = 0; icell < ncells; ++icell) {
199 Double_t cellAmplitude=0, cellTime=0;
200 Short_t cellNumber=0;
201 if (cells->GetCell(icell, cellNumber, cellAmplitude, cellTime) != kTRUE)
202 break;
95f69c66 203 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->New(idigit));
2f7259cf 204 digit->SetId(cellNumber);
2f7259cf 205 digit->SetTime(cellTime);
206 digit->SetTimeR(cellTime);
207 digit->SetIndexInList(idigit);
208 digit->SetType(AliEMCALDigit::kHG);
95f69c66 209 if (fRecalibOnly||fSubBackground) {
b0e354fe 210 Double_t energy = fClusterizer->Calibrate(cellAmplitude,cellTime,cellNumber);
95f69c66 211 digit->SetAmplitude(energy);
212 avgE += energy;
213 } else {
214 digit->SetAmplitude(cellAmplitude);
b0e354fe 215 }
2f7259cf 216 idigit++;
217 }
2f7259cf 218
95f69c66 219 if (fSubBackground) {
220 avgE /= AliEMCALGeometry::GetInstance(fGeomName)->GetNumberOfSuperModules()*48*24;
221 Int_t ndigis = fDigitsArr->GetEntries();
222 for (Int_t i = 0; i < ndigis; ++i) {
223 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(i));
224 Double_t energy = digit->GetAmplitude() - avgE;
225 if (energy<=0.001) {
226 //fDigitsArr->RemoveAt(i);
227 digit->SetAmplitude(0);
228 } else {
229 digit->SetAmplitude(energy);
230 }
b0e354fe 231 }
2f7259cf 232 }
95f69c66 233 fDigitsArr->Compress();
84ccad86 234 fDigitsArr->Sort();
2f7259cf 235}
236
237//________________________________________________________________________________________
95f69c66 238void AliAnalysisTaskEMCALClusterizeFast::RecPoints2Clusters(TClonesArray *clus)
4cf9204b 239{
240 // Cluster energy, global position, cells and their amplitude fractions are restored.
241
95f69c66 242 Bool_t esdobjects = 0;
243 if (strcmp(clus->GetClass()->GetName(),"AliESDCaloCluster")==0)
244 esdobjects = 1;
245
4cf9204b 246 Int_t Ncls = fClusterArr->GetEntriesFast();
247 for(Int_t i=0, nout=clus->GetEntries(); i < Ncls; ++i) {
248 AliEMCALRecPoint *recpoint = static_cast<AliEMCALRecPoint*>(fClusterArr->At(i));
249 Int_t ncells_true = 0;
250 const Int_t ncells = recpoint->GetMultiplicity();
251 UShort_t absIds[ncells];
252 Double32_t ratios[ncells];
253 Int_t *dlist = recpoint->GetDigitsList();
254 Float_t *elist = recpoint->GetEnergiesList();
255 for (Int_t c = 0; c < ncells; ++c) {
256 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(dlist[c]));
257 absIds[ncells_true] = digit->GetId();
258 ratios[ncells_true] = elist[c]/digit->GetAmplitude();
84ccad86 259 if (ratios[ncells_true] < 0.001)
260 continue;
261 ++ncells_true;
4cf9204b 262 }
263
264 if (ncells_true < 1) {
265 AliWarning("Skipping cluster with no cells");
266 continue;
267 }
268
269 // calculate new cluster position
270 TVector3 gpos;
4cf9204b 271 recpoint->GetGlobalPosition(gpos);
84ccad86 272 Float_t g[3];
4cf9204b 273 gpos.GetXYZ(g);
274
95f69c66 275 AliVCluster *c = static_cast<AliVCluster*>(clus->New(nout++));
4cf9204b 276 c->SetType(AliVCluster::kEMCALClusterv1);
277 c->SetE(recpoint->GetEnergy());
278 c->SetPosition(g);
279 c->SetNCells(ncells_true);
4cf9204b 280 c->SetDispersion(recpoint->GetDispersion());
84ccad86 281 c->SetEmcCpvDistance(-1); //not yet implemented
4cf9204b 282 c->SetChi2(-1); //not yet implemented
283 c->SetTOF(recpoint->GetTime()) ; //time-of-flight
284 c->SetNExMax(recpoint->GetNExMax()); //number of local maxima
285 Float_t elipAxis[2];
286 recpoint->GetElipsAxis(elipAxis);
287 c->SetM02(elipAxis[0]*elipAxis[0]) ;
288 c->SetM20(elipAxis[1]*elipAxis[1]) ;
17cff35e 289 c->SetDistanceToBadChannel(recpoint->GetDistanceToBadTower());
95f69c66 290 if (esdobjects) {
291 AliESDCaloCluster *cesd = static_cast<AliESDCaloCluster*>(c);
292 cesd->SetCellsAbsId(absIds);
293 cesd->SetCellsAmplitudeFraction(ratios);
294 } else {
295 AliAODCaloCluster *caod = static_cast<AliAODCaloCluster*>(c);
296 caod->SetCellsAbsId(absIds);
297 caod->SetCellsAmplitudeFraction(ratios);
298 }
84ccad86 299 }
300
301 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(InputEvent());
302 if (!esdevent)
303 return;
304 if (!fRecoUtils)
305 return;
306
95f69c66 307 AliAnalysisManager::GetAnalysisManager()->LoadBranch("Tracks");
84ccad86 308 fRecoUtils->FindMatches(esdevent,clus);
309
95f69c66 310 if (!esdobjects) {
311 Int_t Nclus = clus->GetEntries();
312 for(Int_t i=0; i < Nclus; ++i) {
313 AliAODCaloCluster *c = static_cast<AliAODCaloCluster*>(clus->At(i));
314 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
315 if(trackIndex >= 0) {
316 Float_t dR, dZ;
317 fRecoUtils->GetMatchedResiduals(i,dR, dZ);
318 c->AddTrackMatched(0x0); //esdevent->GetTrack(trackIndex));
319 c->SetTrackDistance(dR,dZ); // not implemented
320 c->SetEmcCpvDistance(dR);
321 c->SetChi2(dZ);
322 if(DebugLevel() > 1)
323 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
324 }
2f7259cf 325 }
95f69c66 326 } else {
84ccad86 327 Int_t Nclus = clus->GetEntries();
328 for(Int_t i=0; i < Nclus; ++i) {
329 AliESDCaloCluster *c = static_cast<AliESDCaloCluster*>(clus->At(i));
2f7259cf 330 Int_t trackIndex = fRecoUtils->GetMatchedTrackIndex(i);
331 if(trackIndex >= 0) {
84ccad86 332 Float_t dR, dZ;
333 fRecoUtils->GetMatchedResiduals(i,dR, dZ);
334 c->SetTrackDistance(dR,dZ);
335 c->SetEmcCpvDistance(dR); //to be consistent with AODs
336 c->SetChi2(dZ); //to be consistent with AODs
4cf9204b 337 TArrayI tm(1,&trackIndex);
338 c->AddTracksMatched(tm);
2f7259cf 339 if(DebugLevel() > 1)
340 AliInfo(Form("Matched Track index %d to new cluster %d\n",trackIndex,i));
341 }
342 }
343 }
344}
345
95f69c66 346//________________________________________________________________________
347void AliAnalysisTaskEMCALClusterizeFast::UpdateCells()
348{
349 // Update cells in case re-calibration was done.
350
351 if (!fCalibData&&!fSubBackground)
352 return;
353
354 AliVCaloCells *cells = InputEvent()->GetEMCALCells();
355 Int_t ncells = cells->GetNumberOfCells();
356 Int_t ndigis = fDigitsArr->GetEntries();
357
358 if (ncells!=ndigis) {
359 cells->DeleteContainer();
360 cells->CreateContainer(ndigis);
361 }
362 for (Int_t idigit = 0; idigit < ndigis; ++idigit) {
363 AliEMCALDigit *digit = static_cast<AliEMCALDigit*>(fDigitsArr->At(idigit));
364 Double_t cellAmplitude = digit->GetCalibAmp();
365 Short_t cellNumber = digit->GetId();
366 Double_t cellTime = digit->GetTime();
367 cells->SetCell(idigit, cellNumber, cellAmplitude, cellTime);
368 }
369}
370
371//________________________________________________________________________
372void AliAnalysisTaskEMCALClusterizeFast::UpdateClusters()
373{
374 // Update cells in case re-calibration was done.
375
376 if (!fAttachClusters)
377 return;
378
379 TClonesArray *clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("caloClusters"));
380 if (!clus)
381 clus = dynamic_cast<TClonesArray*>(InputEvent()->FindListObject("CaloClusters"));
382 if(!clus)
383 return;
384
385 Int_t nents = clus->GetEntries();
386 for (Int_t i=0;i<nents;++i) {
387 AliVCluster *c = static_cast<AliVCluster*>(clus->At(i));
388 if (!c)
389 continue;
390 if (c->IsEMCAL())
391 clus->RemoveAt(i);
392 }
393 clus->Compress();
394 RecPoints2Clusters(clus);
395}
396
2f7259cf 397//________________________________________________________________________________________
398void AliAnalysisTaskEMCALClusterizeFast::Init()
399{
400 //Select clusterization/unfolding algorithm and set all the needed parameters
01bc1ce8 401
2f7259cf 402 AliVEvent * event = InputEvent();
403 if (!event) {
404 AliWarning("Event not available!!!");
405 return;
406 }
407
408 if (event->GetRunNumber()==fRun)
409 return;
410 fRun = event->GetRunNumber();
411
412 if (fJustUnfold){
413 // init the unfolding afterburner
414 delete fUnfolder;
415 fUnfolder = new AliEMCALAfterBurnerUF(fRecParam->GetW0(),fRecParam->GetLocMaxCut());
416 return;
417 }
418
2f7259cf 419 AliEMCALGeometry *geometry = AliEMCALGeometry::GetInstance(fGeomName);
420 if (!geometry) {
421 AliFatal("Geometry not available!!!");
422 return;
423 }
424
425 if (!fGeomMatrixSet) {
426 if (fLoadGeomMatrices) {
e51bc3f5 427 for(Int_t mod=0; mod < geometry->GetNumberOfSuperModules(); ++mod) {
2f7259cf 428 if(fGeomMatrix[mod]){
429 if(DebugLevel() > 2)
430 fGeomMatrix[mod]->Print();
431 geometry->SetMisalMatrix(fGeomMatrix[mod],mod);
432 }
433 }
29b496d5 434 } else { // get matrix from file (work around bug in aliroot)
2f7259cf 435 for(Int_t mod=0; mod < geometry->GetEMCGeometry()->GetNumberOfSuperModules(); ++mod) {
29b496d5 436 const TGeoHMatrix *gm = 0;
437 AliESDEvent *esdevent = dynamic_cast<AliESDEvent*>(event->GetHeader());
438 if (esdevent) {
439 gm = esdevent->GetEMCALMatrix(mod);
440 } else {
441 AliAODHeader *aodheader = dynamic_cast<AliAODHeader*>(event->GetHeader());
442 if (aodheader) {
443 gm = aodheader->GetEMCALMatrix(mod);
444 }
445 }
446 if (gm) {
2f7259cf 447 if(DebugLevel() > 2)
29b496d5 448 gm->Print();
449 geometry->SetMisalMatrix(gm,mod);
2f7259cf 450 }
451 }
452 }
453 fGeomMatrixSet=kTRUE;
454 }
455
456 // setup digit array if needed
457 if (!fDigitsArr) {
458 fDigitsArr = new TClonesArray("AliEMCALDigit", 1000);
459 fDigitsArr->SetOwner(1);
460 }
461
462 // then setup clusterizer
463 delete fClusterizer;
464 if (fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerv1)
465 fClusterizer = new AliEMCALClusterizerv1(geometry);
466 else if(fRecParam->GetClusterizerFlag() == AliEMCALRecParam::kClusterizerNxN)
467 fClusterizer = new AliEMCALClusterizerNxN(geometry);
468 else if(fRecParam->GetClusterizerFlag() > AliEMCALRecParam::kClusterizerNxN) {
469 AliEMCALClusterizerNxN *clusterizer = new AliEMCALClusterizerNxN(geometry);
470 clusterizer->SetNRowDiff(2);
471 clusterizer->SetNColDiff(2);
472 fClusterizer = clusterizer;
473 } else{
474 AliFatal(Form("Clusterizer < %d > not available", fRecParam->GetClusterizerFlag()));
475 }
476 fClusterizer->InitParameters(fRecParam);
01bc1ce8 477
478 if ((!fCalibData&&fLoadCalib) || (!fPedestalData&&fLoadPed)) {
479 AliCDBManager *cdb = AliCDBManager::Instance();
480 if (!cdb->IsDefaultStorageSet() && !fOCDBpath.IsNull())
481 cdb->SetDefaultStorage(fOCDBpath);
482 if (fRun!=cdb->GetRun())
483 cdb->SetRun(fRun);
484 }
485
10d33986 486 if (!fCalibData&&fLoadCalib) {
2f7259cf 487 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Data"));
488 if (entry)
489 fCalibData = static_cast<AliEMCALCalibData*>(entry->GetObject());
490 if (!fCalibData)
491 AliFatal("Calibration parameters not found in CDB!");
492 }
10d33986 493 if (!fPedestalData&&fLoadPed) {
2f7259cf 494 AliCDBEntry *entry = static_cast<AliCDBEntry*>(AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals"));
495 if (entry)
496 fPedestalData = static_cast<AliCaloCalibPedestal*>(entry->GetObject());
497 }
10d33986 498 if (fCalibData) {
499 fClusterizer->SetInputCalibrated(kFALSE);
500 fClusterizer->SetCalibrationParameters(fCalibData);
501 fClusterizer->SetCaloCalibPedestal(fPedestalData);
502 } else {
503 fClusterizer->SetInputCalibrated(kTRUE);
504 }
505 fClusterizer->SetJustClusters(kTRUE);
2f7259cf 506 fClusterizer->SetDigitsArr(fDigitsArr);
507 fClusterizer->SetOutput(0);
508 fClusterArr = const_cast<TObjArray *>(fClusterizer->GetRecPoints());
509}