Fix (Yifei)
[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"
65bec413 59#include "AliEMCALUnfolding.h"
483b0559 60
61ClassImp(AliEMCALClusterizer)
62
63//____________________________________________________________________________
c47157cd 64AliEMCALClusterizer::AliEMCALClusterizer():
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();
81
82}
83
84//____________________________________________________________________________
85AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
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{
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
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 such that the deriving classes can overload
105 // Init() ;
106 //
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
ee08edde 121//____________________________________________________________________________
122AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped):
123 fDigitsArr(NULL),
124 fTreeR(NULL),
125 fRecPoints(NULL),
126 fGeom(geometry),
127 fCalibData(calib),
128 fCaloPed(caloped),
129 fADCchannelECA(0.),fADCpedestalECA(0.),
130 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
131 fDefaultInit(kFALSE),fToUnfold(kFALSE),
132 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
65bec413 133 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
134 fClusterUnfolding(NULL)
ee08edde 135{
c526514b 136 // ctor, geometry and calibration are initialized elsewhere.
137
138 if (!fGeom)
139 AliFatal("Geometry not initialized.");
ee08edde 140
65bec413 141 Int_t i=0;
142 for (i = 0; i < 8; i++)
143 fSSPars[i] = 0.;
144 for (i = 0; i < 3; i++) {
145 fPar5[i] = 0.;
146 fPar6[i] = 0.;
147 }
148
ee08edde 149}
150
151
483b0559 152//____________________________________________________________________________
c47157cd 153AliEMCALClusterizer::~AliEMCALClusterizer()
483b0559 154{
c47157cd 155 // dtor
23ca956b 156 //Already deleted in AliEMCALReconstructor.
157
158// if(fGeom) delete fGeom;
159// if(fCalibData) delete fCalibData;
160// if(fCaloPed) delete fCaloPed;
161
162// if (fDigitsArr) {
163// fDigitsArr->Clear("C");
164// delete fDigitsArr;
165// }
166// if (fRecPoints) {
167// fRecPoints->Delete();
168// delete fRecPoints;
169// }
65bec413 170
171 if(fClusterUnfolding) delete fClusterUnfolding;
18a21c7c 172}
173
ee08edde 174//____________________________________________________________________________
175Float_t AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId)
176{
177
178 // Convert digitized amplitude into energy.
179 // Calibration parameters are taken from calibration data base for raw data,
180 // or from digitizer parameters for simulated data.
181 if(fCalibData){
182
183 if (fGeom==0)
184 AliFatal("Did not get geometry from EMCALLoader") ;
185
186 Int_t iSupMod = -1;
187 Int_t nModule = -1;
188 Int_t nIphi = -1;
189 Int_t nIeta = -1;
190 Int_t iphi = -1;
191 Int_t ieta = -1;
192
193 Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
194 if(!bCell) {
195 fGeom->PrintGeometry();
196 Error("Calibrate()"," Wrong cell id number : %i", absId);
197 assert(0);
198 }
199
200 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
201
202 // Check if channel is bad (dead or hot), in this case return 0.
203 // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
204 // for the moment keep it here but remember to do the selection at the sdigitizer level
205 // and remove it from here
206 Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
207 if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
208 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
209 return 0;
210 }
211 //Check if time is too large or too small, indication of a noisy channel, remove in this case
212 if(time > fTimeMax || time < fTimeMin) return 0;
213
214 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
215 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
216
217 return -fADCpedestalECA + amp * fADCchannelECA ;
218
219 }
220 else //Return energy with default parameters if calibration is not available
221 return -fADCpedestalECA + amp * fADCchannelECA ;
222
223}
224
225//____________________________________________________________________________
226void AliEMCALClusterizer::GetCalibrationParameters()
227{
228 // Set calibration parameters:
229 // if calibration database exists, they are read from database,
230 // otherwise, they are taken from digitizer.
231 //
232 // It is a user responsilibity to open CDB before reconstruction,
233 // for example:
234 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
235
236 //Check if calibration is stored in data base
237
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!");
247
248}
249
250//____________________________________________________________________________
251void AliEMCALClusterizer::GetCaloCalibPedestal()
252{
c526514b 253 // Set calibration parameters:
254 // if calibration database exists, they are read from database,
255 // otherwise, they are taken from digitizer.
256 //
257 // It is a user responsilibity to open CDB before reconstruction,
258 // for example:
259 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
260
261 //Check if calibration is stored in data base
262
263 if(!fCaloPed)
264 {
265 AliCDBEntry *entry = (AliCDBEntry*)
266 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
267 if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
268 }
269
270 if(!fCaloPed)
271 AliFatal("Pedestal info not found in CDB!");
272
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){
ee08edde 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 }
303
ee08edde 304}
305
306//____________________________________________________________________________
307void AliEMCALClusterizer::InitParameters()
308{
309 // Initializes the parameters for the Clusterizer
310 fNumberOfECAClusters = 0 ;
311 fCalibData = 0 ;
312 fCaloPed = 0 ;
313
314 const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
315 if(!recParam) {
316 AliFatal("Reconstruction parameters for EMCAL not set!");
317 } else {
318 fECAClusteringThreshold = recParam->GetClusteringThreshold();
319 fECAW0 = recParam->GetW0();
320 fMinECut = recParam->GetMinECut();
321 fToUnfold = recParam->GetUnfold();
ee08edde 322 fECALocMaxCut = recParam->GetLocMaxCut();
323 fTimeCut = recParam->GetTimeCut();
324 fTimeMin = recParam->GetTimeMin();
325 fTimeMax = recParam->GetTimeMax();
326
327 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",
328 fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
65bec413 329
622e10be 330 if(fToUnfold){
331
332 Int_t i=0;
333 for (i = 0; i < 8; i++) {
334 fSSPars[i] = recParam->GetSSPars(i);
335 }//end of loop over parameters
336 for (i = 0; i < 3; i++) {
337 fPar5[i] = recParam->GetPar5(i);
338 fPar6[i] = recParam->GetPar6(i);
339 }//end of loop over parameters
340
341 fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6);
342
343 for (i = 0; i < 8; i++) {
344 AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
345 }
346 for (i = 0; i < 3; i++) {
347 AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
348 AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
349 }
350
351 }// to unfold
352 }// recparam not null
353
ee08edde 354}
355
356
357//____________________________________________________________________________
358void AliEMCALClusterizer::Print(Option_t * /*option*/)const
359{
360 // Print clusterizer parameters
361
362 TString message("\n") ;
363
364 if( strcmp(GetName(), "") !=0 ){
365
366 // Print parameters
367
368 TString taskName(Version()) ;
369
370 printf("--------------- ");
371 printf("%s",taskName.Data()) ;
372 printf(" ");
373 printf("Clusterizing digits: ");
374 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
375 printf("\n ECA Logarithmic weight = %f", fECAW0);
c526514b 376 if(fToUnfold){
ee08edde 377 printf("\nUnfolding on\n");
c526514b 378 printf("Unfolding parameters: fSSpars: \n");
379 Int_t i=0;
380 for (i = 0; i < 8; i++) {
381 printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
382 }
383 printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
384 for (i = 0; i < 3; i++) {
385 printf("fPar5[%d] = %f \n", i, fPar5[i]);
386 printf("fPar6[%d] = %f \n", i, fPar6[i]);
387 }
388 }
ee08edde 389 else
390 printf("\nUnfolding off\n");
391
392 printf("------------------------------------------------------------------");
393 }
394 else
395 printf("AliEMCALClusterizer not initialized ") ;
396}
397
398//____________________________________________________________________________
399void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
400{
401 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
402 if(strstr(option,"deb")) {
403 printf("PrintRecPoints: Clusterization result:") ;
404
405 printf(" Found %d ECA Rec Points\n ",
406 fRecPoints->GetEntriesFast()) ;
407 }
408
409 if(strstr(option,"all")) {
410 if(strstr(option,"deb")) {
411 printf("\n-----------------------------------------------------------------------\n") ;
412 printf("Clusters in ECAL section\n") ;
413 printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
414 }
415 Int_t index;
416 for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
417 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
51e4198e 418 if(rp){
419 TVector3 globalpos;
420 //rp->GetGlobalPosition(globalpos);
421 TVector3 localpos;
422 rp->GetLocalPosition(localpos);
423 Float_t lambda[2];
424 rp->GetElipsAxis(lambda);
425
426 Int_t nprimaries=0;
427 Int_t * primaries = rp->GetPrimaries(nprimaries);
428
429 if(strstr(option,"deb"))
430 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
431 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
432 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
433 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
434 if(strstr(option,"deb")){
435 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
436 printf("%d ", primaries[iprimary] ) ;
437 }
ee08edde 438 }
439 }
440 }
441
442 if(strstr(option,"deb"))
443 printf("\n-----------------------------------------------------------------------\n");
444 }
445}
446
447//___________________________________________________________________
448void AliEMCALClusterizer::PrintRecoInfo()
449{
450 //print reco version
451 printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
452
453}
454
18a21c7c 455//____________________________________________________________________________
c47157cd 456void AliEMCALClusterizer::SetInput(TTree *digitsTree)
18a21c7c 457{
c47157cd 458 // Read the digits from the input tree
459 TBranch *branch = digitsTree->GetBranch("EMCAL");
460 if (!branch) {
461 AliError("can't get the branch with the EMCAL digits !");
462 return;
463 }
3f9ad99f 464 if (!fDigitsArr)
465 fDigitsArr = new TClonesArray("AliEMCALDigit",100);
c47157cd 466 branch->SetAddress(&fDigitsArr);
467 branch->GetEntry(0);
483b0559 468}
469
470//____________________________________________________________________________
c47157cd 471void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
483b0559 472{
c47157cd 473 // Read the digits from the input tree
474 fTreeR = clustersTree;
475
476 AliDebug(9, "Making array for EMCAL clusters");
477 fRecPoints = new TObjArray(100) ;
478 Int_t split = 0;
479 Int_t bufsize = 32000;
480 fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
88cb7938 481}
ee08edde 482
483
484