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