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"
59 #include "AliEMCALUnfolding.h"
61 ClassImp(AliEMCALClusterizer)
63 //____________________________________________________________________________
64 AliEMCALClusterizer::AliEMCALClusterizer():
71 fADCchannelECA(0.),fADCpedestalECA(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)
84 //____________________________________________________________________________
85 AliEMCALClusterizer::AliEMCALClusterizer(AliEMCALGeometry* geometry):
92 fADCchannelECA(0.),fADCpedestalECA(0.),
93 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
94 fDefaultInit(kFALSE),fToUnfold(kFALSE),
95 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
96 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
97 fClusterUnfolding(NULL)
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
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
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, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped):
129 fADCchannelECA(0.),fADCpedestalECA(0.),
130 fTimeMin(-1.),fTimeMax(1.),fTimeCut(1.),
131 fDefaultInit(kFALSE),fToUnfold(kFALSE),
132 fNumberOfECAClusters(0), fECAClusteringThreshold(0.),
133 fECALocMaxCut(0.),fECAW0(0.),fMinECut(0.),
134 fClusterUnfolding(NULL)
136 // ctor, geometry and calibration are initialized elsewhere.
139 AliFatal("Geometry not initialized.");
142 for (i = 0; i < 8; i++)
144 for (i = 0; i < 3; i++) {
152 //____________________________________________________________________________
153 AliEMCALClusterizer::~AliEMCALClusterizer()
156 //Already deleted in AliEMCALReconstructor.
158 // if(fGeom) delete fGeom;
159 // if(fCalibData) delete fCalibData;
160 // if(fCaloPed) delete fCaloPed;
163 // fDigitsArr->Clear("C");
164 // delete fDigitsArr;
167 // fRecPoints->Delete();
168 // delete fRecPoints;
171 if(fClusterUnfolding) delete fClusterUnfolding;
174 //____________________________________________________________________________
175 Float_t AliEMCALClusterizer::Calibrate(const Float_t amp, const Float_t time, const Int_t absId)
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.
184 AliFatal("Did not get geometry from EMCALLoader") ;
193 Bool_t bCell = fGeom->GetCellIndex(absId, iSupMod, nModule, nIphi, nIeta) ;
195 fGeom->PrintGeometry();
196 Error("Calibrate()"," Wrong cell id number : %i", absId);
200 fGeom->GetCellPhiEtaIndexInSModule(iSupMod,nModule,nIphi, nIeta,iphi,ieta);
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));
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;
214 fADCchannelECA = fCalibData->GetADCchannel (iSupMod,ieta,iphi);
215 fADCpedestalECA = fCalibData->GetADCpedestal(iSupMod,ieta,iphi);
217 return -fADCpedestalECA + amp * fADCchannelECA ;
220 else //Return energy with default parameters if calibration is not available
221 return -fADCpedestalECA + amp * fADCchannelECA ;
225 //____________________________________________________________________________
226 void AliEMCALClusterizer::GetCalibrationParameters()
228 // Set calibration parameters:
229 // if calibration database exists, they are read from database,
230 // otherwise, they are taken from digitizer.
232 // It is a user responsilibity to open CDB before reconstruction,
234 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
236 //Check if calibration is stored in data base
240 AliCDBEntry *entry = (AliCDBEntry*)
241 AliCDBManager::Instance()->Get("EMCAL/Calib/Data");
242 if (entry) fCalibData = (AliEMCALCalibData*) entry->GetObject();
246 AliFatal("Calibration parameters not found in CDB!");
250 //____________________________________________________________________________
251 void AliEMCALClusterizer::GetCaloCalibPedestal()
253 // Set calibration parameters:
254 // if calibration database exists, they are read from database,
255 // otherwise, they are taken from digitizer.
257 // It is a user responsilibity to open CDB before reconstruction,
259 // AliCDBStorage* storage = AliCDBManager::Instance()->GetStorage("local://CalibDB");
261 //Check if calibration is stored in data base
265 AliCDBEntry *entry = (AliCDBEntry*)
266 AliCDBManager::Instance()->Get("EMCAL/Calib/Pedestals");
267 if (entry) fCaloPed = (AliCaloCalibPedestal*) entry->GetObject();
271 AliFatal("Pedestal info not found in CDB!");
275 //____________________________________________________________________________
276 void AliEMCALClusterizer::Init()
278 // Make all memory allocations which can not be done in default constructor.
279 // Attach the Clusterizer task to the list of EMCAL tasks
281 AliRunLoader *rl = AliRunLoader::Instance();
282 if (rl->GetAliRun()){
283 AliEMCAL* emcal = dynamic_cast<AliEMCAL*>(rl->GetAliRun()->GetDetector("EMCAL"));
284 if(emcal)fGeom = emcal->GetGeometry();
288 fGeom = AliEMCALGeometry::GetInstance(AliEMCALGeometry::GetDefaultGeometryName());
291 AliDebug(1,Form("geom %p",fGeom));
294 gMinuit = new TMinuit(100) ;
297 for (i = 0; i < 8; i++)
299 for (i = 0; i < 3; i++) {
306 //____________________________________________________________________________
307 void AliEMCALClusterizer::InitParameters()
309 // Initializes the parameters for the Clusterizer
310 fNumberOfECAClusters = 0 ;
314 const AliEMCALRecParam* recParam = AliEMCALReconstructor::GetRecParam();
316 AliFatal("Reconstruction parameters for EMCAL not set!");
318 fECAClusteringThreshold = recParam->GetClusteringThreshold();
319 fECAW0 = recParam->GetW0();
320 fMinECut = recParam->GetMinECut();
321 fToUnfold = recParam->GetUnfold();
322 if(fToUnfold) AliWarning("Cluster Unfolding ON. Implementing only for eta=0 case!!!");
323 fECALocMaxCut = recParam->GetLocMaxCut();
324 fTimeCut = recParam->GetTimeCut();
325 fTimeMin = recParam->GetTimeMin();
326 fTimeMax = recParam->GetTimeMax();
328 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",
329 fECAClusteringThreshold,fECAW0,fMinECut,fToUnfold,fECALocMaxCut,fTimeCut, fTimeMin, fTimeMax));
334 for (i = 0; i < 8; i++) {
335 fSSPars[i] = recParam->GetSSPars(i);
336 }//end of loop over parameters
337 for (i = 0; i < 3; i++) {
338 fPar5[i] = recParam->GetPar5(i);
339 fPar6[i] = recParam->GetPar6(i);
340 }//end of loop over parameters
342 fClusterUnfolding=new AliEMCALUnfolding(fGeom,fECALocMaxCut,fSSPars,fPar5,fPar6);
344 for (i = 0; i < 8; i++) {
345 AliDebug(1,Form("unfolding shower shape parameters: fSSPars=%f \n",fSSPars[i]));
347 for (i = 0; i < 3; i++) {
348 AliDebug(1,Form("unfolding parameter 5: fPar5=%f \n",fPar5[i]));
349 AliDebug(1,Form("unfolding parameter 6: fPar5=%f \n",fPar6[i]));
357 //____________________________________________________________________________
358 void AliEMCALClusterizer::Print(Option_t * /*option*/)const
360 // Print clusterizer parameters
362 TString message("\n") ;
364 if( strcmp(GetName(), "") !=0 ){
368 TString taskName(Version()) ;
370 printf("--------------- ");
371 printf("%s",taskName.Data()) ;
373 printf("Clusterizing digits: ");
374 printf("\n ECA Local Maximum cut = %f", fECALocMaxCut);
375 printf("\n ECA Logarithmic weight = %f", fECAW0);
377 printf("\nUnfolding on\n");
379 printf("\nUnfolding off\n");
381 printf("------------------------------------------------------------------");
384 printf("AliEMCALClusterizer not initialized ") ;
387 //____________________________________________________________________________
388 void AliEMCALClusterizer::PrintRecPoints(Option_t * option)
390 // Prints list of RecPoints produced at the current pass of AliEMCALClusterizer
391 if(strstr(option,"deb")) {
392 printf("PrintRecPoints: Clusterization result:") ;
394 printf(" Found %d ECA Rec Points\n ",
395 fRecPoints->GetEntriesFast()) ;
398 if(strstr(option,"all")) {
399 if(strstr(option,"deb")) {
400 printf("\n-----------------------------------------------------------------------\n") ;
401 printf("Clusters in ECAL section\n") ;
402 printf("Index Ene(GeV) Multi Module GX GY GZ lX lY lZ Dispersion Lambda 1 Lambda 2 # of prim Primaries list\n") ;
405 for (index = 0 ; index < fRecPoints->GetEntries() ; index++) {
406 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint * >(fRecPoints->At(index)) ;
409 //rp->GetGlobalPosition(globalpos);
411 rp->GetLocalPosition(localpos);
413 rp->GetElipsAxis(lambda);
416 Int_t * primaries = rp->GetPrimaries(nprimaries);
418 if(strstr(option,"deb"))
419 printf("\n%6d %8.4f %3d %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4.1f %4f %4f %2d : ",
420 rp->GetIndexInList(), rp->GetEnergy(), rp->GetMultiplicity(),
421 globalpos.X(), globalpos.Y(), globalpos.Z(), localpos.X(), localpos.Y(), localpos.Z(),
422 rp->GetDispersion(), lambda[0], lambda[1], nprimaries) ;
423 if(strstr(option,"deb")){
424 for (Int_t iprimary=0; iprimary<nprimaries; iprimary++) {
425 printf("%d ", primaries[iprimary] ) ;
431 if(strstr(option,"deb"))
432 printf("\n-----------------------------------------------------------------------\n");
436 //___________________________________________________________________
437 void AliEMCALClusterizer::PrintRecoInfo()
440 printf(" AliEMCALClusterizer::PrintRecoInfo() : version %s \n", Version() );
444 //____________________________________________________________________________
445 void AliEMCALClusterizer::SetInput(TTree *digitsTree)
447 // Read the digits from the input tree
448 TBranch *branch = digitsTree->GetBranch("EMCAL");
450 AliError("can't get the branch with the EMCAL digits !");
454 fDigitsArr = new TClonesArray("AliEMCALDigit",100);
455 branch->SetAddress(&fDigitsArr);
459 //____________________________________________________________________________
460 void AliEMCALClusterizer::SetOutput(TTree *clustersTree)
462 // Read the digits from the input tree
463 fTreeR = clustersTree;
465 AliDebug(9, "Making array for EMCAL clusters");
466 fRecPoints = new TObjArray(100) ;
468 Int_t bufsize = 32000;
469 fTreeR->Branch("EMCALECARP", "TObjArray", &fRecPoints, bufsize, split);