1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //_________________________________________________________________________
19 // Base class for the clusterization algorithm (pure abstract)
21 //*-- Author: Yves Schutz SUBATECH
23 // Clusterization mother class. Contains common methods/data members of different
24 // clusterizers. GCB 2010
25 //////////////////////////////////////////////////////////////////////////////
27 // --- ROOT system ---
35 #include <TBenchmark.h>
39 // --- Standard library ---
42 // --- AliRoot header files ---
43 #include "AliEMCALClusterizer.h"
44 #include "AliEMCALReconstructor.h"
45 #include "AliRunLoader.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"
57 #include "AliCDBEntry.h"
59 ClassImp(AliEMCALClusterizer)
61 //____________________________________________________________________________
62 AliEMCALClusterizer::AliEMCALClusterizer():
63 fIsInputCalibrated(kFALSE),
64 fJustClusters(kFALSE),
71 fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
72 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
73 fDefaultInit(kFALSE),fToUnfold(kFALSE),
74 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
75 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
76 fClusterUnfolding(NULL)
83 //____________________________________________________________________________
84 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
85 fIsInputCalibrated(kFALSE),
86 fJustClusters(kFALSE),
93 fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
94 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
95 fDefaultInit(kFALSE),fToUnfold(kFALSE),
96 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
97 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
98 fClusterUnfolding(NULL)
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();
109 AliFatal("Geometry not initialized.");
112 for (i = 0; i < 8; i++)
114 for (i = 0; i < 3; i++) {
120 //____________________________________________________________________________
121 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry *geometry,
122 AliEMCALCalibData *calib,
123 AliCaloCalibPedestal *caloped):
124 fIsInputCalibrated(kFALSE),
125 fJustClusters(kFALSE),
132 fADCchannelECA(0.),fADCpedestalECA(0.), fTimeECA(0.),
133 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
134 fDefaultInit(kFALSE),fToUnfold(kFALSE),
135 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
136 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
137 fClusterUnfolding(NULL)
139 // ctor, geometry and calibration are initialized elsewhere.
142 AliFatal("Geometry not initialized.");
145 for (i = 0; i < 8; i++)
147 for (i = 0; i < 3; i++) {
153 //____________________________________________________________________________
154 AliEMCALClusterizer::~AliEMCALClusterizer()
157 //Already deleted in AliEMCALReconstructor.
159 if(fClusterUnfolding) delete fClusterUnfolding;
161 // make sure we delete the rec points array
164 //Delete digits array
169 //____________________________________________________________________________
170 void AliEMCALClusterizer::DeleteRecPoints()
172 // free the cluster array
175 AliDebug(2, "Deleting fRecPoints.");
176 fRecPoints->Delete();
182 //____________________________________________________________________________
183 void AliEMCALClusterizer::DeleteDigits()
185 // free the digits array
188 AliDebug(2, "Deleting fDigitsArr.");
189 fDigitsArr->Clear("C");
195 //____________________________________________________________________________
196 void AliEMCALClusterizer::Calibrate(Float_t & amp, Float_t & time, const Int_t absId)
198 // Convert digitized amplitude into energy, calibrate time
199 // Calibration parameters are taken from OCDB : OCDB/EMCAL/Calib/Data
201 //Return energy with default parameters if calibration is not available
202 if (!fCalibData && !fCaloPed) {
203 if (fIsInputCalibrated == kTRUE)
205 AliDebug(10, Form("Input already calibrated!"));
209 AliFatal("OCDB calibration and bad map parameters are not available");
215 AliFatal("Did not get geometry from EMCALLoader") ;
224 Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
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.
236 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
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
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));
252 if (fIsInputCalibrated || !fCalibData)
254 AliDebug(10, Form("Input already calibrated!"));
258 Int_t bc = 0; // Get somehow the bunch crossing number
260 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
261 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
262 fTimeECA = fCalibData->GetTimeChannel(iSupMod,ieta,iphi, bc);
265 amp = amp * fADCchannelECA - fADCpedestalECA ;
269 //____________________________________________________________________________
270 void AliEMCALClusterizer::GetCalibrationParameters()
272 // Set calibration parameters:
273 // If calibration database exists, they are read from database,
274 // otherwise, they are taken from digitizer.
275 // It is a user responsilibity to open CDB before reconstruction,
277 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
279 if (fIsInputCalibrated)
282 //Check if calibration is stored in data base
285 AliCDBEntry *entry = (AliCDBEntry*)
286 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
287 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
291 AliFatal("Calibration parameters not found in CDB!");
294 //____________________________________________________________________________
295 void AliEMCALClusterizer::GetCaloCalibPedestal()
297 // Set calibration parameters:
298 // if calibration database exists, they are read from database,
299 // otherwise, they are taken from digitizer.
301 // It is a user responsilibity to open CDB before reconstruction,
303 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
305 if (fIsInputCalibrated)
308 // Check if calibration is stored in data base
311 AliCDBEntry *entry = (AliCDBEntry*)
312 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
313 if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
317 AliFatal("Pedestal info not found in CDB!");
320 //____________________________________________________________________________
321 void AliEMCALClusterizer::Init()
323 // Make all memory allocations which can not be done in default constructor.
324 // Attach the Clusterizer task to the list of EMCAL tasks
326 AliRunLoader *rl = AliRunLoader::Instance();
327 if (rl->GetAliRun()){
328 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
329 if(emcal)fGeom = emcal->GetGeometry();
333 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
336 AliDebug(1,Form("geom %p",fGeom));
339 gMinuit = new TMinuit(100) ;
342 for (i = 0; i < 8; i++)
344 for (i = 0; i < 3; i++) {
350 //____________________________________________________________________________
351 void AliEMCALClusterizer::InitParameters()
353 // Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().
355 return InitParameters(AliEMCALReconstructor::GetRecParam());
358 //____________________________________________________________________________
359 void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
361 // Initializes the parameters for the Clusterizer
363 fNumberOfECAClusters = 0 ;
368 AliFatal("Reconstruction parameters for EMCAL not set!");
371 fECAClusteringThreshold = recParam->GetClusteringThreshold();
372 fECAW0 = recParam->GetW0();
373 fMinECut = recParam->GetMinECut();
374 fToUnfold = recParam->GetUnfold();
375 fECALocMaxCut = recParam->GetLocMaxCut();
376 fTimeCut = recParam->GetTimeCut();
377 fTimeMin = recParam->GetTimeMin();
378 fTimeMax = recParam->GetTimeMax();
381 SetNRowDiff(recParam->GetNRowDiff());
382 SetNColDiff(recParam->GetNColDiff());
384 AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
385 "fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
386 fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
390 for (i = 0; i < 8; i++) {
391 fSSPars[i] = recParam->GetSSPars(i);
392 } //end of loop over parameters
393 for (i = 0; i < 3; i++) {
394 fPar5[i] = recParam->GetPar5(i);
395 fPar6[i] = recParam->GetPar6(i);
396 } //end of loop over parameters
398 InitClusterUnfolding();
400 for (i = 0; i < 8; i++) {
401 AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
403 for (i = 0; i < 3; i++) {
404 AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
405 AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
410 //____________________________________________________________________________
411 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
413 // Print clusterizer parameters
415 TString message("\n") ;
417 if (strcmp(GetName(),"") == 0) {
418 printf("AliEMCALClusterizer not initialized\n");
423 TString taskName(Version()) ;
425 printf("--------------- ");
426 printf("%s",taskName.Data()) ;
428 printf("Clusterizing digits: ");
429 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
430 printf("\n ECA Logarithmic weight = %f", fECAW0);
432 printf("\nUnfolding on\n");
433 printf("Unfolding parameters: fSSpars: \n");
435 for (i = 0; i < 8; i++) {
436 printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
438 printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
439 for (i = 0; i < 3; i++) {
440 printf("fPar5[%d] = %f \n", i, fPar5[i]);
441 printf("fPar6[%d] = %f \n", i, fPar6[i]);
445 printf("\nUnfolding off\n");
447 printf("------------------------------------------------------------------");
450 //____________________________________________________________________________
451 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
453 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
455 if (strstr(option,"deb")) {
456 printf("PrintRecPoints: Clusterization result:") ;
457 printf(" Found %d ECA Rec Points\n ",
458 fRecPoints->GetEntriesFast()) ;
461 if (strstr(option,"all")) {
462 if (strstr(option,"deb")) {
463 printf("\n-----------------------------------------------------------------------\n") ;
464 printf("Clusters in ECAL section\n") ;
465 printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
468 for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
469 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
474 //rp->GetGlobalPosition(globalpos);
476 rp->GetLocalPosition(localpos);
478 rp->GetElipsAxis(lambda);
481 Int_t * primaries = rp->GetPrimaries(nprimaries);
482 if(strstr(option,"deb"))
483 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
484 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
485 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
486 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
487 if(strstr(option,"deb")){
488 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
489 printf("%d ", primaries[iprimary] ) ;
494 if(strstr(option,"deb"))
495 printf("\n-----------------------------------------------------------------------\n");
499 //___________________________________________________________________
500 void AliEMCALClusterizer::PrintRecoInfo()
502 // Print reco version
504 printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
507 //____________________________________________________________________________
508 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
510 // Read the digits from the input tree
512 TBranch *branch = digitsTree->GetBranch("EMCAL");
514 AliError("can't get the branch with the EMCAL digits !");
518 fDigitsArr = new TClonesArray("AliEMCALDigit",100);
519 branch->SetAddress(&fDigitsArr);
523 //____________________________________________________________________________
524 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
526 // Read the digits from the input tree
528 AliDebug(9, "Making array for EMCAL clusters");
529 fRecPoints = new TObjArray(1000);
531 fTreeR = clustersTree;
533 Int_t bufsize = 32000;
534 fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
538 //___________________________________________________________________
539 void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
541 // Flag to indicate that input is calibrated - the case when we run already on ESD
543 fIsInputCalibrated = val;
546 //___________________________________________________________________
547 void AliEMCALClusterizer::SetJustClusters(Bool_t val)
549 // Flag to indicate that we are running on ESDs, when calling
550 // rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes