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