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