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 ---
28 #include "TClonesArray.h"
36 #include <TBenchmark.h>
40 // --- Standard library ---
43 // --- AliRoot header files ---
44 #include "AliEMCALClusterizer.h"
45 #include "AliEMCALReconstructor.h"
46 #include "AliRunLoader.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"
58 #include "AliCDBEntry.h"
60 ClassImp(AliEMCALClusterizer)
62 //____________________________________________________________________________
63 AliEMCALClusterizer::AliEMCALClusterizer():
64 fIsInputCalibrated(kFALSE),
65 fJustClusters(kFALSE),
72 fADCchannelECA(0.),fADCpedestalECA(0.),
73 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
74 fDefaultInit(kFALSE),fToUnfold(kFALSE),
75 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
76 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
77 fClusterUnfolding(NULL)
84 //____________________________________________________________________________
85 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
86 fIsInputCalibrated(kFALSE),
87 fJustClusters(kFALSE),
94 fADCchannelECA(0.),fADCpedestalECA(0.),
95 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
96 fDefaultInit(kFALSE),fToUnfold(kFALSE),
97 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
98 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
99 fClusterUnfolding(NULL)
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 // Note for the future: the use on runloader should be avoided or optional at least.
105 // Another way is to make Init virtual and protected at least
106 // such that the deriving classes can overload Init();
110 AliFatal("Geometry not initialized.");
113 for (i = 0; i < 8; i++)
115 for (i = 0; i < 3; i++) {
121 //____________________________________________________________________________
122 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry *geometry,
123 AliEMCALCalibData *calib,
124 AliCaloCalibPedestal *caloped):
125 fIsInputCalibrated(kFALSE),
126 fJustClusters(kFALSE),
133 fADCchannelECA(0.),fADCpedestalECA(0.),
134 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
135 fDefaultInit(kFALSE),fToUnfold(kFALSE),
136 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
137 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
138 fClusterUnfolding(NULL)
140 // ctor, geometry and calibration are initialized elsewhere.
143 AliFatal("Geometry not initialized.");
146 for (i = 0; i < 8; i++)
148 for (i = 0; i < 3; i++) {
154 //____________________________________________________________________________
155 AliEMCALClusterizer::~AliEMCALClusterizer()
158 //Already deleted in AliEMCALReconstructor.
160 if(fClusterUnfolding) delete fClusterUnfolding;
163 //____________________________________________________________________________
164 Float_t AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId)
166 // Convert digitized amplitude into energy.
167 // Calibration parameters are taken from calibration data base for raw data,
168 // or from digitizer parameters for simulated data.
170 //Return energy with default parameters if calibration is not available
171 if (!fCalibData&&!fCaloPed) {
172 if (fIsInputCalibrated == kTRUE)
174 AliDebug(10, Form("Input already calibrated!"));
178 return -fADCpedestalECA + amp * fADCchannelECA ;
182 AliFatal("Did not get geometry from EMCALLoader") ;
191 Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
193 fGeom->PrintGeometry();
194 AliError(Form("Wrong cell id number : %i", absId));
195 //assert(0); // GCB: This aborts reconstruction of raw simulations
196 //where simulation had more SM than default geometry,
197 //change to return 0, to avoid aborting good generations.
201 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
203 // Check if channel is bad (dead or hot), in this case return 0.
204 // Gustavo: 15-12-09 In case of RAW data this selection is already done, but not in simulation.
205 // for the moment keep it here but remember to do the selection at the sdigitizer level
206 // and remove it from here
208 Int_t channelStatus = (Int_t)(fCaloPed->GetDeadMap(iSupMod))->GetBinContent(ieta,iphi);
209 if(channelStatus == AliCaloCalibPedestal::kHot || channelStatus == AliCaloCalibPedestal::kDead) {
210 AliDebug(2,Form("Tower from SM %d, ieta %d, iphi %d is BAD : status %d !!!",iSupMod,ieta,iphi, channelStatus));
214 //Check if time is too large or too small, indication of a noisy channel, remove in this case
215 if(time > fTimeMax || time < fTimeMin) return 0;
217 if (fIsInputCalibrated||!fCalibData)
219 AliDebug(10, Form("Input already calibrated!"));
223 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
224 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
226 return -fADCpedestalECA + amp * fADCchannelECA ;
229 //____________________________________________________________________________
230 void AliEMCALClusterizer::GetCalibrationParameters()
232 // Set calibration parameters:
233 // If calibration database exists, they are read from database,
234 // otherwise, they are taken from digitizer.
235 // It is a user responsilibity to open CDB before reconstruction,
237 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
239 if (fIsInputCalibrated)
242 //Check if calibration is stored in data base
245 AliCDBEntry *entry = (AliCDBEntry*)
246 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
247 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
251 AliFatal("Calibration parameters not found in CDB!");
254 //____________________________________________________________________________
255 void AliEMCALClusterizer::GetCaloCalibPedestal()
257 // Set calibration parameters:
258 // if calibration database exists, they are read from database,
259 // otherwise, they are taken from digitizer.
261 // It is a user responsilibity to open CDB before reconstruction,
263 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
265 if (fIsInputCalibrated)
268 // Check if calibration is stored in data base
271 AliCDBEntry *entry = (AliCDBEntry*)
272 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
273 if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
277 AliFatal("Pedestal info not found in CDB!");
280 //____________________________________________________________________________
281 void AliEMCALClusterizer::Init()
283 // Make all memory allocations which can not be done in default constructor.
284 // Attach the Clusterizer task to the list of EMCAL tasks
286 AliRunLoader *rl = AliRunLoader::Instance();
287 if (rl->GetAliRun()){
288 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
289 if(emcal)fGeom = emcal->GetGeometry();
293 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
296 AliDebug(1,Form("geom %p",fGeom));
299 gMinuit = new TMinuit(100) ;
302 for (i = 0; i < 8; i++)
304 for (i = 0; i < 3; i++) {
310 //____________________________________________________________________________
311 void AliEMCALClusterizer::InitParameters()
313 // Initializes the parameters for the Clusterizer from AliEMCALReconstructor::GetRecParam().
315 return InitParameters(AliEMCALReconstructor::GetRecParam());
318 //____________________________________________________________________________
319 void AliEMCALClusterizer::InitParameters(const AliEMCALRecParam* recParam)
321 // Initializes the parameters for the Clusterizer
323 fNumberOfECAClusters = 0 ;
328 AliFatal("Reconstruction parameters for EMCAL not set!");
331 fECAClusteringThreshold = recParam->GetClusteringThreshold();
332 fECAW0 = recParam->GetW0();
333 fMinECut = recParam->GetMinECut();
334 fToUnfold = recParam->GetUnfold();
335 fECALocMaxCut = recParam->GetLocMaxCut();
336 fTimeCut = recParam->GetTimeCut();
337 fTimeMin = recParam->GetTimeMin();
338 fTimeMax = recParam->GetTimeMax();
340 AliDebug(1,Form("Reconstruction parameters: fECAClusteringThreshold=%.3f GeV, fECAW=%.3f, fMinECut=%.3f GeV, "
341 "fToUnfold=%d, fECALocMaxCut=%.3f GeV, fTimeCut=%e s,fTimeMin=%e s,fTimeMax=%e s",
342 fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
346 for (i = 0; i < 8; i++) {
347 fSSPars[i] = recParam->GetSSPars(i);
348 } //end of loop over parameters
349 for (i = 0; i < 3; i++) {
350 fPar5[i] = recParam->GetPar5(i);
351 fPar6[i] = recParam->GetPar6(i);
352 } //end of loop over parameters
354 InitClusterUnfolding();
356 for (i = 0; i < 8; i++) {
357 AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
359 for (i = 0; i < 3; i++) {
360 AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
361 AliDebug(1,Form("unfolding parameter 6: fPar6=%f \n",fPar6[i]));
366 //____________________________________________________________________________
367 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
369 // Print clusterizer parameters
371 TString message("\n") ;
373 if (strcmp(GetName(),"") == 0) {
374 printf("AliEMCALClusterizer not initialized\n");
379 TString taskName(Version()) ;
381 printf("--------------- ");
382 printf("%s",taskName.Data()) ;
384 printf("Clusterizing digits: ");
385 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
386 printf("\n ECA Logarithmic weight = %f", fECAW0);
388 printf("\nUnfolding on\n");
389 printf("Unfolding parameters: fSSpars: \n");
391 for (i = 0; i < 8; i++) {
392 printf("fSSPars[%d] = %f \n", i, fSSPars[i]);
394 printf("Unfolding parameter 5 and 6: fPar5 and fPar6: \n");
395 for (i = 0; i < 3; i++) {
396 printf("fPar5[%d] = %f \n", i, fPar5[i]);
397 printf("fPar6[%d] = %f \n", i, fPar6[i]);
401 printf("\nUnfolding off\n");
403 printf("------------------------------------------------------------------");
406 //____________________________________________________________________________
407 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
409 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
411 if (strstr(option,"deb")) {
412 printf("PrintRecPoints: Clusterization result:") ;
413 printf(" Found %d ECA Rec Points\n ",
414 fRecPoints->GetEntriesFast()) ;
417 if (strstr(option,"all")) {
418 if (strstr(option,"deb")) {
419 printf("\n-----------------------------------------------------------------------\n") ;
420 printf("Clusters in ECAL section\n") ;
421 printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
424 for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
425 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
430 //rp->GetGlobalPosition(globalpos);
432 rp->GetLocalPosition(localpos);
434 rp->GetElipsAxis(lambda);
437 Int_t * primaries = rp->GetPrimaries(nprimaries);
438 if(strstr(option,"deb"))
439 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
440 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
441 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
442 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
443 if(strstr(option,"deb")){
444 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
445 printf("%d ", primaries[iprimary] ) ;
450 if(strstr(option,"deb"))
451 printf("\n-----------------------------------------------------------------------\n");
455 //___________________________________________________________________
456 void AliEMCALClusterizer::PrintRecoInfo()
458 // Print reco version
460 printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
463 //____________________________________________________________________________
464 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
466 // Read the digits from the input tree
468 TBranch *branch = digitsTree->GetBranch("EMCAL");
470 AliError("can't get the branch with the EMCAL digits !");
474 fDigitsArr = new TClonesArray("AliEMCALDigit",100);
475 branch->SetAddress(&fDigitsArr);
479 //____________________________________________________________________________
480 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
482 // Read the digits from the input tree
484 AliDebug(9, "Making array for EMCAL clusters");
485 fRecPoints = new TObjArray(1000);
487 fTreeR = clustersTree;
489 Int_t bufsize = 32000;
490 fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);
494 //___________________________________________________________________
495 void AliEMCALClusterizer::SetInputCalibrated(Bool_t val)
497 // Flag to indicate that input is calibrated - the case when we run already on ESD
499 fIsInputCalibrated = val;
502 //___________________________________________________________________
503 void AliEMCALClusterizer::SetJustClusters(Bool_t val)
505 // Flag to indicate that we are running on ESDs, when calling
506 // rp->EvalAll(fECAW0,fDigitsArr,fJustClusters); in derived classes