Bug fix in the constructor (thanks to A.Marin)
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizer.cxx
CommitLineData
483b0559 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//_________________________________________________________________________
19// Base class for the clusterization algorithm (pure abstract)
20//*--
21//*-- Author: Yves Schutz SUBATECH
ee08edde 22//
23// Clusterization mother class. Contains common methods/data members of different
24// clusterizers. GCB 2010
483b0559 25//////////////////////////////////////////////////////////////////////////////
26
27// --- ROOT system ---
c47157cd 28#include "TClonesArray.h"
29#include "TTree.h"
ee08edde 30#include <TFile.h>
31class TFolder;
32#include <TMath.h>
33#include <TMinuit.h>
34#include <TTree.h>
35class TSystem;
36#include <TBenchmark.h>
37#include <TBrowser.h>
38#include <TROOT.h>
483b0559 39
40// --- Standard library ---
ee08edde 41#include <cassert>
483b0559 42
43// --- AliRoot header files ---
483b0559 44#include "AliEMCALClusterizer.h"
ee08edde 45#include "AliEMCALReconstructor.h"
46#include "AliRunLoader.h"
47#include "AliRun.h"
c47157cd 48#include "AliLog.h"
ee08edde 49#include "AliEMCAL.h"
50#include "AliEMCALRecPoint.h"
51#include "AliEMCALRecParam.h"
52#include "AliEMCALGeometry.h"
53#include "AliEMCALRecParam.h"
54#include "AliCDBManager.h"
55#include "AliCaloCalibPedestal.h"
56#include "AliEMCALCalibData.h"
57class AliCDBStorage;
58#include "AliCDBEntry.h"
483b0559 59
60ClassImp(AliEMCALClusterizer)
61
62//____________________________________________________________________________
c47157cd 63AliEMCALClusterizer::AliEMCALClusterizer():
ac084be7 64 fIsInputCalibrated(kFALSE),
294553d1 65 fJustClusters(kFALSE),
c47157cd 66 fDigitsArr(NULL),
67 fTreeR(NULL),
ee08edde 68 fRecPoints(NULL),
69 fGeom(NULL),
70 fCalibData(NULL),
71 fCaloPed(NULL),
e93968a8 72 fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
ee08edde 73 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
74 fDefaultInit(kFALSE),fToUnfold(kFALSE),
75 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413 76 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
77 fClusterUnfolding(NULL)
483b0559 78{
79 // ctor
ee08edde 80
81 Init();
ee08edde 82}
83
84//____________________________________________________________________________
85AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
ac084be7 86 fIsInputCalibrated(kFALSE),
294553d1 87 fJustClusters(kFALSE),
ee08edde 88 fDigitsArr(NULL),
89 fTreeR(NULL),
90 fRecPoints(NULL),
91 fGeom(geometry),
92 fCalibData(NULL),
93 fCaloPed(NULL),
e93968a8 94 fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
ee08edde 95 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
96 fDefaultInit(kFALSE),fToUnfold(kFALSE),
97 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413 98 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
99 fClusterUnfolding(NULL)
ee08edde 100{
ac084be7 101 // Ctor with the indication of the file where header Tree and digits Tree are stored.
102 // Use this contructor to avoid usage of Init() which uses runloader.
103 // Change needed by HLT - MP.
104 // Note for the future: the use on runloader should be avoided or optional at least.
105 // Another way is to make Init virtual and protected at least
106 // such that the deriving classes can overload Init();
ee08edde 107
108 if (!fGeom)
109 {
110 AliFatal("Geometry not initialized.");
111 }
65bec413 112 Int_t i=0;
113 for (i = 0; i < 8; i++)
114 fSSPars[i] = 0.;
115 for (i = 0; i < 3; i++) {
116 fPar5[i] = 0.;
117 fPar6[i] = 0.;
118 }
483b0559 119}
839828a6 120
483b0559 121//____________________________________________________________________________
ac084be7 122AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry *geometry,
123 AliEMCALCalibData *calib,
124 AliCaloCalibPedestal *caloped):
125 fIsInputCalibrated(kFALSE),
294553d1 126 fJustClusters(kFALSE),
ee08edde 127 fDigitsArr(NULL),
128 fTreeR(NULL),
129 fRecPoints(NULL),
130 fGeom(geometry),
131 fCalibData(calib),
132 fCaloPed(caloped),
e93968a8 133 fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
ee08edde 134 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
135 fDefaultInit(kFALSE),fToUnfold(kFALSE),
136 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413 137 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
138 fClusterUnfolding(NULL)
ee08edde 139{
c526514b 140 // ctor, geometry and calibration are initialized elsewhere.
141
142 if (!fGeom)
143 AliFatal("Geometry not initialized.");
ee08edde 144
65bec413 145 Int_t i=0;
146 for (i = 0; i < 8; i++)
147 fSSPars[i] = 0.;
148 for (i = 0; i < 3; i++) {
149 fPar5[i] = 0.;
150 fPar6[i] = 0.;
151 }
ee08edde 152}
153
ee08edde 154//____________________________________________________________________________
c47157cd 155AliEMCALClusterizer::~AliEMCALClusterizer()
483b0559 156{
c47157cd 157 // dtor
23ca956b 158 //Already deleted in AliEMCALReconstructor.
159
65bec413 160 if(fClusterUnfolding) delete fClusterUnfolding;
c828bc97 161
162 // make sure we delete the rec points array
163 DeleteRecPoints();
8d0210ea 164
165 //Delete digits array
166 DeleteDigits();
167
c828bc97 168}
169
170//____________________________________________________________________________
171void AliEMCALClusterizer::DeleteRecPoints()
172{
173 // free the cluster array
174 if (fRecPoints)
175 {
176 AliDebug(2, "Deleting fRecPoints.");
177 fRecPoints->Delete();
178 delete fRecPoints;
8d0210ea 179 }
180}
181
182//____________________________________________________________________________
183void AliEMCALClusterizer::DeleteDigits()
184{
185 // free the digits array
186 if (fDigitsArr)
187 {
188 AliDebug(2, "Deleting fDigitsArr.");
189 fDigitsArr->Clear("C");
190 delete fDigitsArr;
c828bc97 191 }
18a21c7c 192}
193
194//____________________________________________________________________________
783153ff 195void AliEMCALClusterizer::Calibrate(Float_t & amp, Float_t & time, const Int_t absId)
ee08edde 196{
783153ff 197 // Convert digitized amplitude into energy, calibrate time
198 // Calibration parameters are taken from OCDB : OCDB/EMCAL/Calib/Data
ac084be7 199
200 //Return energy with default parameters if calibration is not available
e93968a8 201 if (!fCalibData && !fCaloPed) {
ac084be7 202 if (fIsInputCalibrated == kTRUE)
b1324a01 203 {
204 AliDebug(10, Form("Input already calibrated!"));
783153ff 205 return ;
e93968a8 206 }
207 else{
208 AliFatal("OCDB calibration and bad map parameters are not available");
209 return;
210 }
ac084be7 211 }
212
213 if (fGeom==0)
214 AliFatal("Did not get geometry from EMCALLoader") ;
b1324a01 215
ac084be7 216 Int_t iSupMod = -1;
217 Int_t nModule = -1;
218 Int_t nIphi = -1;
219 Int_t nIeta = -1;
220 Int_t iphi = -1;
221 Int_t ieta = -1;
ee08edde 222
ac084be7 223 Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
224 if(!bCell) {
225 fGeom->PrintGeometry();
226 AliError(Form("Wrong cell id number : %i", absId));
227 //assert(0); // GCB: This aborts reconstruction of raw simulations
228 //where simulation had more SM than default geometry,
229 //change to return 0, to avoid aborting good generations.
783153ff 230 amp = 0;
231 time = 0;
232 return ;
ee08edde 233 }
b1324a01 234
ac084be7 235 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
236
237 // Check if channel is bad (dead or hot), in this case return 0.
238 // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
239 // for the moment keep it here but remember to do the selection at the sdigitizer level
240 // and remove it from here
294553d1 241 if (fCaloPed) {
242 Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
243 if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
244 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
783153ff 245 amp = 0 ;
246 time = 0 ;
247 return ;
294553d1 248 }
ac084be7 249 }
250 //Check if time is too large or too small, indication of a noisy channel, remove in this case
783153ff 251 if(time > fTimeMax || time < fTimeMin) {
252 amp = 0 ;
253 time = 0 ;
254 return ;
255 }
b1324a01 256
9c8525ff 257 if (fIsInputCalibrated||!fCalibData)
ac084be7 258 {
259 AliDebug(10, Form("Input already calibrated!"));
783153ff 260 return ;
b1324a01 261 }
ac084be7 262
263 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
264 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
783153ff 265 fTimeECA = fCalibData->GetTimeChannel(iSupMod,ieta,iphi);
266
267 time -= fTimeECA ;
268 amp = amp * fADCchannelECA - fADCpedestalECA ;
269
ee08edde 270}
271
272//____________________________________________________________________________
273void AliEMCALClusterizer::GetCalibrationParameters()
274{
275 // Set calibration parameters:
ac084be7 276 // If calibration database exists, they are read from database,
ee08edde 277 // otherwise, they are taken from digitizer.
ee08edde 278 // It is a user responsilibity to open CDB before reconstruction,
279 // for example:
280 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
ac084be7 281
282 if (fIsInputCalibrated)
283 return;
ee08edde 284
285 //Check if calibration is stored in data base
ee08edde 286 if(!fCalibData)
287 {
288 AliCDBEntry *entry = (AliCDBEntry*)
289 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
290 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
291 }
292
293 if(!fCalibData)
294 AliFatal("Calibration parameters not found in CDB!");
ee08edde 295}
296
297//____________________________________________________________________________
298void AliEMCALClusterizer::GetCaloCalibPedestal()
299{
c526514b 300 // Set calibration parameters:
301 // if calibration database exists, they are read from database,
302 // otherwise, they are taken from digitizer.
303 //
304 // It is a user responsilibity to open CDB before reconstruction,
305 // for example:
306 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
ac084be7 307
308 if (fIsInputCalibrated)
309 return;
c526514b 310
ac084be7 311 // Check if calibration is stored in data base
c526514b 312 if(!fCaloPed)
313 {
314 AliCDBEntry *entry = (AliCDBEntry*)
315 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
316 if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
317 }
318
319 if(!fCaloPed)
320 AliFatal("Pedestal info not found in CDB!");
ee08edde 321}
322
323//____________________________________________________________________________
324void AliEMCALClusterizer::Init()
325{
326 // Make all memory allocations which can not be done in default constructor.
327 // Attach the Clusterizer task to the list of EMCAL tasks
328
329 AliRunLoader *rl = AliRunLoader::Instance();
a51e676d 330 if (rl->GetAliRun()){
331 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
332 if(emcal)fGeom = emcal->GetGeometry();
333 }
334
335 if(!fGeom){
ac084be7 336 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
a51e676d 337 }
ee08edde 338
339 AliDebug(1,Form("geom %p",fGeom));
340
341 if(!gMinuit)
342 gMinuit = new TMinuit(100) ;
343
65bec413 344 Int_t i=0;
345 for (i = 0; i < 8; i++)
346 fSSPars[i] = 0.;
347 for (i = 0; i < 3; i++) {
348 fPar5[i] = 0.;
349 fPar6[i] = 0.;
350 }
ee08edde 351}
352
353//____________________________________________________________________________
354void AliEMCALClusterizer::InitParameters()
355{
ac084be7 356 // Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().
357
358 return InitParameters(AliEMCALReconstructor::GetRecParam());
359}
360
361//____________________________________________________________________________
362void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
363{
ee08edde 364 // Initializes the parameters for the Clusterizer
ac084be7 365
ee08edde 366 fNumberOfECAClusters = 0 ;
367 fCalibData = 0 ;
368 fCaloPed = 0 ;
369
ee08edde 370 if(!recParam) {
371 AliFatal("Reconstruction parameters for EMCAL not set!");
ac084be7 372 }
373
374 fECAClusteringThreshold = recParam->GetClusteringThreshold();
375 fECAW0 = recParam->GetW0();
376 fMinECut = recParam->GetMinECut();
377 fToUnfold = recParam->GetUnfold();
378 fECALocMaxCut = recParam->GetLocMaxCut();
379 fTimeCut = recParam->GetTimeCut();
380 fTimeMin = recParam->GetTimeMin();
381 fTimeMax = recParam->GetTimeMax();
65bec413 382
ac084be7 383 AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
384 "fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
385 fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
386
387 if (fToUnfold) {
388 Int_t i=0;
389 for (i = 0; i < 8; i++) {
390 fSSPars[i] = recParam->GetSSPars(i);
391 } //end of loop over parameters
392 for (i = 0; i < 3; i++) {
393 fPar5[i] = recParam->GetPar5(i);
394 fPar6[i] = recParam->GetPar6(i);
395 } //end of loop over parameters
622e10be 396
ac084be7 397 InitClusterUnfolding();
622e10be 398
ac084be7 399 for (i = 0; i < 8; i++) {
400 AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
401 }
402 for (i = 0; i < 3; i++) {
403 AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
404 AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
405 }
406 } // to unfold
ee08edde 407}
408
ee08edde 409//____________________________________________________________________________
410void AliEMCALClusterizer::Print(Option_t * /*option*/)const
411{
412 // Print clusterizer parameters
413
414 TString message("\n") ;
415
ac084be7 416 if (strcmp(GetName(),"") == 0) {
417 printf("AliEMCALClusterizer not initialized\n");
418 return;
419 }
ee08edde 420
ac084be7 421 // Print parameters
422 TString taskName(Version()) ;
ee08edde 423
ac084be7 424 printf("--------------- ");
425 printf("%s",taskName.Data()) ;
426 printf(" ");
427 printf("Clusterizing digits: ");
428 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
429 printf("\n ECA Logarithmic weight = %f", fECAW0);
430 if (fToUnfold) {
431 printf("\nUnfolding on\n");
432 printf("Unfolding parameters: fSSpars: \n");
433 Int_t i=0;
434 for (i = 0; i < 8; i++) {
435 printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
436 }
437 printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
438 for (i = 0; i < 3; i++) {
439 printf("fPar5[%d] = %f \n", i, fPar5[i]);
440 printf("fPar6[%d] = %f \n", i, fPar6[i]);
c526514b 441 }
ee08edde 442 }
443 else
ac084be7 444 printf("\nUnfolding off\n");
445
446 printf("------------------------------------------------------------------");
ee08edde 447}
448
449//____________________________________________________________________________
450void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
451{
452 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
ac084be7 453
454 if (strstr(option,"deb")) {
ee08edde 455 printf("PrintRecPoints: Clusterization result:") ;
ee08edde 456 printf(" Found %d ECA Rec Points\n ",
457 fRecPoints->GetEntriesFast()) ;
458 }
459
ac084be7 460 if (strstr(option,"all")) {
461 if (strstr(option,"deb")) {
ee08edde 462 printf("\n-----------------------------------------------------------------------\n") ;
463 printf("Clusters in ECAL section\n") ;
464 printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
465 }
466 Int_t index;
467 for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
468 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
ac084be7 469 if (!rp)
470 continue;
471
472 TVector3 globalpos;
473 //rp->GetGlobalPosition(globalpos);
474 TVector3 localpos;
475 rp->GetLocalPosition(localpos);
476 Float_t lambda[2];
477 rp->GetElipsAxis(lambda);
51e4198e 478
ac084be7 479 Int_t nprimaries=0;
480 Int_t * primaries = rp->GetPrimaries(nprimaries);
481 if(strstr(option,"deb"))
482 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
483 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
484 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
485 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
486 if(strstr(option,"deb")){
487 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
488 printf("%d ", primaries[iprimary] ) ;
ee08edde 489 }
490 }
491 }
492
493 if(strstr(option,"deb"))
494 printf("\n-----------------------------------------------------------------------\n");
495 }
496}
497
498//___________________________________________________________________
499void AliEMCALClusterizer::PrintRecoInfo()
500{
ac084be7 501 // Print reco version
502
ee08edde 503 printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
ee08edde 504}
505
506//____________________________________________________________________________
c47157cd 507void AliEMCALClusterizer::SetInput(TTree *digitsTree)
18a21c7c 508{
c47157cd 509 // Read the digits from the input tree
ac084be7 510
c47157cd 511 TBranch *branch = digitsTree->GetBranch("EMCAL");
512 if (!branch) {
513 AliError("can't get the branch with the EMCAL digits !");
514 return;
515 }
3f9ad99f 516 if (!fDigitsArr)
517 fDigitsArr = new TClonesArray("AliEMCALDigit",100);
c47157cd 518 branch->SetAddress(&fDigitsArr);
519 branch->GetEntry(0);
483b0559 520}
521
522//____________________________________________________________________________
c47157cd 523void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
483b0559 524{
c47157cd 525 // Read the digits from the input tree
ac084be7 526
c47157cd 527 AliDebug(9, "Making array for EMCAL clusters");
ac084be7 528 fRecPoints = new TObjArray(1000);
529 if (clustersTree) {
530 fTreeR = clustersTree;
531 Int_t split = 0;
532 Int_t bufsize = 32000;
533 fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
534 }
88cb7938 535}
ee08edde 536
b1324a01 537//___________________________________________________________________
ac084be7 538void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
b1324a01 539{
ac084be7 540 // Flag to indicate that input is calibrated - the case when we run already on ESD
ee08edde 541
ac084be7 542 fIsInputCalibrated = val;
543}
294553d1 544
545//___________________________________________________________________
546void AliEMCALClusterizer::SetJustClusters(Bool_t val)
547{
548 // Flag to indicate that we are running on ESDs, when calling
549 // rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes
550
551 fJustClusters = val;
552}