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 **************************************************************************/
16 // This class derives from AliEMCALClustrerizer
22 #include <TBenchmark.h>
25 #include <TClonesArray.h>
30 #include "AliEMCALRecPoint.h"
31 #include "AliEMCALDigit.h"
32 #include "AliEMCALGeometry.h"
33 #include "AliCaloCalibPedestal.h"
34 #include "AliEMCALCalibData.h"
35 #include "AliESDCaloCluster.h"
36 #include "AliEMCALUnfolding.h"
38 #include "AliEMCALClusterizerFixedWindow.h"
40 ClassImp(AliEMCALClusterizerFixedWindow)
42 //__________________________________________________________________________________________
43 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow() :
44 AliEMCALClusterizer(),
69 //__________________________________________________________________________________________
70 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry) :
71 AliEMCALClusterizer(geometry),
96 //__________________________________________________________________________________________
97 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped) :
98 AliEMCALClusterizer(geometry, calib, caloped),
104 fNEtaDigitsSupMod(0),
105 fNPhiDigitsSupMod(0),
123 //__________________________________________________________________________________________
124 AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
128 if (fClustersArray) {
129 for (Int_t i = 0; i < fNTotalClus; i++) {
130 if (fClustersArray[i]) {
131 delete[] fClustersArray[i];
132 fClustersArray[i] = 0;
134 delete[] fClustersArray;
141 //__________________________________________________________________________________________
142 void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n)
144 // Set fNphi; if clusterizer already initialized gives a warning and does nothing
146 if (fInitialized != 0)
147 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
152 //__________________________________________________________________________________________
153 void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n)
155 // Set fNeta; if clusterizer already initialized gives a warning and does nothing
157 if (fInitialized != 0)
158 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
163 //__________________________________________________________________________________________
164 void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s)
166 // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
168 if (fInitialized != 0)
169 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
174 //__________________________________________________________________________________________
175 void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s)
177 // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
179 if (fInitialized != 0)
180 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
185 //__________________________________________________________________________________________
186 void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b)
188 // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
190 if (fInitialized != 0)
191 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
197 //__________________________________________________________________________________________
198 void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
200 // Steering method to perform clusterization for the current event
202 static Float_t cputime = 0;
203 static Float_t realtime = 0;
205 if (strstr(option,"tim"))
206 gBenchmark->Start("EMCALClusterizer");
208 if (strstr(option,"print"))
211 //Get calibration parameters from file or digitizer default values.
212 GetCalibrationParameters();
214 //Get dead channel map from file or digitizer default values.
215 GetCaloCalibPedestal();
217 MakeClusters(); //only the real clusters
220 fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
221 fClusterUnfolding->MakeUnfolding();
224 //Evaluate position, dispersion and other RecPoint properties for EC section
225 for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
226 AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
228 rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
229 AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
230 //For each rec.point set the distance to the nearest bad crystal
232 rp->EvalDistanceToBadChannels(fCaloPed);
238 for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
239 AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
241 rp->SetIndexInList(index);
243 else AliFatal("RecPoint NULL!!");
249 if (strstr(option,"deb") || strstr(option,"all"))
250 PrintRecPoints(option);
252 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
254 if (strstr(option,"tim")) {
255 gBenchmark->Stop("EMCALClusterizer");
256 Printf("Exec took %f CPU time (%f real time) for clusterizing",
257 gBenchmark->GetCpuTime("EMCALClusterizer")-cputime,gBenchmark->GetRealTime("EMCALClusterizer")-realtime);
258 cputime = gBenchmark->GetCpuTime("EMCALClusterizer");
259 realtime = gBenchmark->GetRealTime("EMCALClusterizer");
263 //__________________________________________________________________________________________
264 void AliEMCALClusterizerFixedWindow::ExecOnce()
266 // Initialize clusterizer.
271 AliError("Did not get geometry!");
275 // Defining geometry and clusterization parameter
276 fNEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
277 fNPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
282 fNEtaDigits = fNEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
283 fNPhiDigits = fNPhiDigitsSupMod * fGeom->GetNPhiSuperModule();
286 fNTRUPhi = fGeom->GetNPhiSuperModule() * 3;
287 fNTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
288 fNEtaDigits /= fNTRUEta;
289 fNPhiDigits /= fNTRUPhi;
292 // Check if clusterizer parameter are compatible with calorimeter geometry
293 if (fNEtaDigits < fNeta){
294 AliError(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,fNEtaDigits));
297 if (fNPhiDigits < fNphi){
298 AliError(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,fNPhiDigits));
301 if (fNEtaDigits % fShiftEta != 0){
302 AliError(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,fNEtaDigits));
305 if (fNPhiDigits % fShiftPhi != 0){
306 AliError(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,fNPhiDigits));
309 if (fNeta % fShiftEta != 0){
310 AliError(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
313 if (fNphi % fShiftPhi != 0){
314 AliError(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
318 fMaxShiftPhi = fNphi / fShiftPhi;
319 fMaxShiftEta = fNeta / fShiftEta;
321 fNClusEtaNoShift = fNEtaDigits / fNeta;
322 fNClusPhiNoShift = fNPhiDigits / fNphi;
324 fNClusters = fNClusEtaNoShift * fNClusPhiNoShift * fNTRUEta * fNTRUPhi;
326 fNTotalClus = fNClusters * fMaxShiftEta * fMaxShiftPhi;
328 fNDigitsCluster = fNphi * fNeta;
330 if (fClustersArray) {
331 for (Int_t i = 0; i < fNTotalClus; i++) {
332 if (fClustersArray[i]) {
333 delete[] fClustersArray[i];
334 fClustersArray[i] = 0;
336 delete[] fClustersArray;
341 fClustersArray = new AliEMCALDigit**[fNTotalClus];
342 for (Int_t i = 0; i < fNTotalClus; i++) {
343 fClustersArray[i] = new AliEMCALDigit*[fNDigitsCluster];
344 for (Int_t j = 0; j < fNDigitsCluster; j++) {
345 fClustersArray[i][j] = 0;
349 AliDebug(1,Form("****ExecOnce*****\n"
350 "fNphi = %d, fNeta = %d, fShiftPhi = %d, fShiftEta = %d, fTRUshift = %d\n"
351 "fNEtaDigitsSupMod = %d, fNPhiDigitsSupMod = %d, fNTRUPhi = %d, fNTRUEta = %d, fNEtaDigits = %d, fNPhiDigits = %d\n"
352 "fMaxShiftPhi = %d, fMaxShiftEta = %d, fNDigitsCluster = %d, fNClusEtaNoShift = %d, fNClusPhiNoShift = %d\n"
353 "fNClusters = %d, fNTotalClus = %d\n",
354 fNphi,fNeta,fShiftPhi,fShiftEta,fTRUshift,
355 fNEtaDigitsSupMod,fNPhiDigitsSupMod,fNTRUPhi,fNTRUEta,fNEtaDigits,fNPhiDigits,
356 fMaxShiftPhi,fMaxShiftEta,fNDigitsCluster,fNClusEtaNoShift,fNClusPhiNoShift,
357 fNClusters,fNTotalClus));
362 //__________________________________________________________________________________________
363 void AliEMCALClusterizerFixedWindow::MakeClusters()
367 fNumberOfECAClusters = 0;
368 fRecPoints->Delete();
370 if (fInitialized == 0)
373 if (fInitialized == -1) {
374 AliError(Form("%s: error initializing the clusterizer. No clusterization will be performed.",GetName()));
378 // Set up TObjArray with pointers to digits to work on calibrated digits
379 TObjArray *digitsC = new TObjArray();
380 AliEMCALDigit *digit;
381 Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
382 TIter nextdigit(fDigitsArr);
383 while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigit()))) { // calibrate and clean up digits
384 dEnergyCalibrated = digit->GetAmplitude();
385 time = digit->GetTime();
386 Calibrate(dEnergyCalibrated, time, digit->GetId());
387 digit->SetCalibAmp(dEnergyCalibrated);
388 digit->SetTime(time);
389 if (dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin) {
392 else if (!fGeom->CheckAbsCellId(digit->GetId())) {
396 ehs += dEnergyCalibrated;
397 digitsC->AddLast(digit);
401 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
402 fDigitsArr->GetEntries(),fMinECut,ehs));
404 Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
405 Int_t iphi=0, ieta=0; // cell eta-phi indexes in SM
407 for (Int_t ishiftPhi = 0; ishiftPhi < fMaxShiftPhi; ishiftPhi++){
408 Int_t nClusPhi = (fNPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
410 for (Int_t ishiftEta = 0; ishiftEta < fMaxShiftEta; ishiftEta++) {
412 Int_t nClusEta = (fNEtaDigits - fShiftEta * ishiftEta) / fNeta;
414 Int_t iTotalClus = fNClusters * (ishiftPhi * fMaxShiftEta + ishiftEta);
416 TIter nextdigitC(digitsC);
417 while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
419 fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
420 fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
422 Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + fNPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
424 Int_t iTRUphi = iphi_eff / fNPhiDigits;
426 iphi_eff -= iTRUphi * fNPhiDigits;
428 Int_t iClusPhi = iphi_eff / fNphi;
430 if (iphi_eff < 0 || iClusPhi >= nClusPhi)
433 Int_t ieta_eff = ieta - fShiftEta * ishiftEta + fNEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
435 Int_t iTRUeta = ieta_eff / fNEtaDigits;
437 ieta_eff -= iTRUeta * fNEtaDigits;
439 Int_t iClusEta = ieta_eff / fNeta;
441 if (ieta_eff < 0 || iClusEta >= nClusEta)
444 iphi_eff += iTRUphi * fNPhiDigits;
445 iClusPhi = iphi_eff / fNphi;
447 ieta_eff += iTRUeta * fNEtaDigits;
448 iClusEta = ieta_eff / fNeta;
450 Int_t iCluster = iClusPhi + iClusEta * fNClusPhiNoShift * fNTRUPhi;
451 Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
453 if (iCluster >= fNClusters){
454 AliError(Form("iCluster out of range! iCluster = %d, fNClusters = %d (should never happen...)", iCluster, fNClusters));
458 iCluster += iTotalClus;
460 if (iCluster >= fNTotalClus){
461 AliError(Form("iCluster out of range! iCluster = %d, fNTotalClus = %d (should never happen...)", iCluster, fNTotalClus));
465 if (iDigit >= fNDigitsCluster){
466 AliError(Form("iDigit out of range! iDigit = %d, fNDigitsCluster = %d (should never happen...)", iDigit, fNDigitsCluster));
470 if (fClustersArray[iCluster][iDigit] != 0){
471 AliError("Digit already added! (should never happen...)");
475 fClustersArray[iCluster][iDigit] = digit;
479 } // loop on eta shift
481 } // loop on phi shift
483 AliEMCALRecPoint *recPoint = 0;
484 Bool_t recPointOk = kFALSE;
485 for (Int_t iCluster = 0; iCluster < fNTotalClus; iCluster++) {
488 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(fNumberOfECAClusters*2+1);
490 (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("");
491 recPoint = static_cast<AliEMCALRecPoint*>(fRecPoints->At(fNumberOfECAClusters));
495 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
496 recPoint->SetUniqueID(iCluster);
497 fNumberOfECAClusters++;
499 for (Int_t iDigit = 0; iDigit < fNDigitsCluster; iDigit++) {
500 if (fClustersArray[iCluster][iDigit] == 0) continue;
501 digit = fClustersArray[iCluster][iDigit];
502 recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
503 fClustersArray[iCluster][iDigit] = 0;
507 if (recPointOk) { // unset the pointer so that a new rec point will be allocated in the next iteration
513 AliError("Error allocating rec points!");
519 fRecPoints->RemoveLast();
520 fNumberOfECAClusters--;
524 AliDebug(1, Form("MakeClusters: Number of digits %d -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
525 AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast()));