new file to test SDigits
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerFixedWindow.cxx
CommitLineData
89504c57 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
27bad386 16// This class derives from AliEMCALClustrerizer
89504c57 17
27bad386 18// --- Root ---
89504c57 19#include <TMath.h>
20#include <TMinuit.h>
21#include <TTree.h>
22#include <TBenchmark.h>
23#include <TBrowser.h>
24#include <TROOT.h>
25#include <TClonesArray.h>
26#include <TH1I.h>
27
27bad386 28// --- AliRoot ---
89504c57 29#include "AliLog.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"
89504c57 37
27bad386 38#include "AliEMCALClusterizerFixedWindow.h"
39
89504c57 40ClassImp(AliEMCALClusterizerFixedWindow)
41
27bad386 42//__________________________________________________________________________________________
85d48773 43AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow() :
44 AliEMCALClusterizer(),
45 fNphi(4),
46 fNeta(4),
47 fShiftPhi(2),
48 fShiftEta(2),
49 fTRUshift(0),
50 fNEtaDigitsSupMod(0),
51 fNPhiDigitsSupMod(0),
52 fNTRUPhi(0),
53 fNTRUEta(0),
54 fNEtaDigits(0),
55 fNPhiDigits(0),
56 fMaxShiftPhi(0),
57 fMaxShiftEta(0),
58 fNDigitsCluster(0),
59 fNClusEtaNoShift(0),
60 fNClusPhiNoShift(0),
61 fNClusters(0),
62 fNTotalClus(0),
63 fClustersArray(0),
64 fInitialized(0)
89504c57 65{
06ca320f 66 // Constructor
89504c57 67}
68
27bad386 69//__________________________________________________________________________________________
85d48773 70AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry) :
71 AliEMCALClusterizer(geometry),
72 fNphi(4),
73 fNeta(4),
74 fShiftPhi(2),
75 fShiftEta(2),
76 fTRUshift(0),
77 fNEtaDigitsSupMod(0),
78 fNPhiDigitsSupMod(0),
79 fNTRUPhi(0),
80 fNTRUEta(0),
81 fNEtaDigits(0),
82 fNPhiDigits(0),
83 fMaxShiftPhi(0),
84 fMaxShiftEta(0),
85 fNDigitsCluster(0),
86 fNClusEtaNoShift(0),
87 fNClusPhiNoShift(0),
88 fNClusters(0),
89 fNTotalClus(0),
90 fClustersArray(0),
91 fInitialized(0)
89504c57 92{
27bad386 93 // Constructor
89504c57 94}
95
27bad386 96//__________________________________________________________________________________________
85d48773 97AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped) :
98 AliEMCALClusterizer(geometry, calib, caloped),
99 fNphi(4),
100 fNeta(4),
101 fShiftPhi(2),
102 fShiftEta(2),
103 fTRUshift(0),
104 fNEtaDigitsSupMod(0),
105 fNPhiDigitsSupMod(0),
106 fNTRUPhi(0),
107 fNTRUEta(0),
108 fNEtaDigits(0),
109 fNPhiDigits(0),
110 fMaxShiftPhi(0),
111 fMaxShiftEta(0),
112 fNDigitsCluster(0),
113 fNClusEtaNoShift(0),
114 fNClusPhiNoShift(0),
115 fNClusters(0),
116 fNTotalClus(0),
117 fClustersArray(0),
118 fInitialized(0)
89504c57 119{
27bad386 120 // Constructor
89504c57 121}
122
27bad386 123//__________________________________________________________________________________________
89504c57 124AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
125{
06ca320f 126 // Destructor
85d48773 127
128 if (fClustersArray) {
129 for (Int_t i = 0; i < fNTotalClus; i++) {
130 if (fClustersArray[i]) {
131 delete[] fClustersArray[i];
132 fClustersArray[i] = 0;
133 }
134 delete[] fClustersArray;
135 fClustersArray = 0;
136 }
137 }
138
27bad386 139}
140
141//__________________________________________________________________________________________
142void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n)
143{
144 // Set fNphi; if clusterizer already initialized gives a warning and does nothing
145
85d48773 146 if (fInitialized != 0)
27bad386 147 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
148 else
149 fNphi = n;
150}
151
152//__________________________________________________________________________________________
153void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n)
154{
155 // Set fNeta; if clusterizer already initialized gives a warning and does nothing
156
85d48773 157 if (fInitialized != 0)
27bad386 158 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
159 else
160 fNeta = n;
161}
162
163//__________________________________________________________________________________________
164void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s)
165{
166 // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
167
85d48773 168 if (fInitialized != 0)
27bad386 169 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
170 else
171 fShiftPhi = s;
172}
173
174//__________________________________________________________________________________________
175void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s)
176{
177 // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
178
85d48773 179 if (fInitialized != 0)
27bad386 180 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
181 else
182 fShiftEta = s;
183}
184
185//__________________________________________________________________________________________
186void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b)
187{
188 // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
189
85d48773 190 if (fInitialized != 0)
27bad386 191 AliWarning("Clusterizer already initialized. Unable to change the parameters.");
192 else
193 fTRUshift = b;
89504c57 194}
195
27bad386 196
197//__________________________________________________________________________________________
89504c57 198void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
199{
06ca320f 200 // Steering method to perform clusterization for the current event
85d48773 201
202 static Float_t cputime = 0;
203 static Float_t realtime = 0;
204
06ca320f 205 if (strstr(option,"tim"))
206 gBenchmark->Start("EMCALClusterizer");
89504c57 207
06ca320f 208 if (strstr(option,"print"))
209 Print("");
89504c57 210
06ca320f 211 //Get calibration parameters from file or digitizer default values.
212 GetCalibrationParameters();
89504c57 213
06ca320f 214 //Get dead channel map from file or digitizer default values.
215 GetCaloCalibPedestal();
89504c57 216
06ca320f 217 MakeClusters(); //only the real clusters
89504c57 218
06ca320f 219 if (fToUnfold) {
220 fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
221 fClusterUnfolding->MakeUnfolding();
222 }
89504c57 223
06ca320f 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));
227 if (rp) {
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
231 if (fCaloPed)
232 rp->EvalDistanceToBadChannels(fCaloPed);
233 }
234 }
6265e793 235
236 fRecPoints->Sort();
89504c57 237
06ca320f 238 for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
239 AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
240 if (rp) {
241 rp->SetIndexInList(index);
242 }
243 else AliFatal("RecPoint NULL!!");
244 }
89504c57 245
06ca320f 246 if (fTreeR)
247 fTreeR->Fill();
89504c57 248
06ca320f 249 if (strstr(option,"deb") || strstr(option,"all"))
250 PrintRecPoints(option);
89504c57 251
06ca320f 252 AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
89504c57 253
06ca320f 254 if (strstr(option,"tim")) {
255 gBenchmark->Stop("EMCALClusterizer");
85d48773 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");
06ca320f 260 }
89504c57 261}
262
27bad386 263//__________________________________________________________________________________________
85d48773 264void AliEMCALClusterizerFixedWindow::ExecOnce()
89504c57 265{
85d48773 266 // Initialize clusterizer.
267
268 fInitialized = -1;
269
270 if (!fGeom) {
271 AliError("Did not get geometry!");
272 return;
273 }
89504c57 274
06ca320f 275 // Defining geometry and clusterization parameter
85d48773 276 fNEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
277 fNPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
89504c57 278
85d48773 279 fNTRUPhi = 1;
280 fNTRUEta = 1;
89504c57 281
85d48773 282 fNEtaDigits = fNEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
283 fNPhiDigits = fNPhiDigitsSupMod * fGeom->GetNPhiSuperModule();
89504c57 284
27bad386 285 if (fTRUshift){
85d48773 286 fNTRUPhi = fGeom->GetNPhiSuperModule() * 3;
287 fNTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
288 fNEtaDigits /= fNTRUEta;
289 fNPhiDigits /= fNTRUPhi;
89504c57 290 }
291
292 // Check if clusterizer parameter are compatible with calorimeter geometry
85d48773 293 if (fNEtaDigits < fNeta){
294 AliError(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,fNEtaDigits));
06ca320f 295 return;
296 }
85d48773 297 if (fNPhiDigits < fNphi){
298 AliError(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,fNPhiDigits));
06ca320f 299 return;
300 }
85d48773 301 if (fNEtaDigits % fShiftEta != 0){
302 AliError(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,fNEtaDigits));
06ca320f 303 return;
304 }
85d48773 305 if (fNPhiDigits % fShiftPhi != 0){
306 AliError(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,fNPhiDigits));
06ca320f 307 return;
308 }
309 if (fNeta % fShiftEta != 0){
85d48773 310 AliError(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
06ca320f 311 return;
312 }
313 if (fNphi % fShiftPhi != 0){
85d48773 314 AliError(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
06ca320f 315 return;
316 }
89504c57 317
85d48773 318 fMaxShiftPhi = fNphi / fShiftPhi;
319 fMaxShiftEta = fNeta / fShiftEta;
89504c57 320
85d48773 321 fNClusEtaNoShift = fNEtaDigits / fNeta;
322 fNClusPhiNoShift = fNPhiDigits / fNphi;
323
324 fNClusters = fNClusEtaNoShift * fNClusPhiNoShift * fNTRUEta * fNTRUPhi;
325
326 fNTotalClus = fNClusters * fMaxShiftEta * fMaxShiftPhi;
327
328 fNDigitsCluster = fNphi * fNeta;
329
330 if (fClustersArray) {
331 for (Int_t i = 0; i < fNTotalClus; i++) {
332 if (fClustersArray[i]) {
333 delete[] fClustersArray[i];
334 fClustersArray[i] = 0;
335 }
336 delete[] fClustersArray;
337 fClustersArray = 0;
338 }
339 }
340
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;
89504c57 346 }
347 }
85d48773 348
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));
358
359 fInitialized = 1;
360}
361
362//__________________________________________________________________________________________
363void AliEMCALClusterizerFixedWindow::MakeClusters()
364{
365 // Make clusters.
366
367 fNumberOfECAClusters = 0;
368 fRecPoints->Delete();
369
370 if (fInitialized == 0)
371 ExecOnce();
372
373 if (fInitialized == -1) {
374 AliError(Form("%s: error initializing the clusterizer. No clusterization will be performed.",GetName()));
375 return;
376 }
89504c57 377
e93ec1f7 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);
959d5ca4 389 if (dEnergyCalibrated < fMinECut || time > fTimeMax || time < fTimeMin) {
e93ec1f7 390 continue;
391 }
392 else if (!fGeom->CheckAbsCellId(digit->GetId())) {
393 continue;
394 }
395 else {
396 ehs += dEnergyCalibrated;
397 digitsC->AddLast(digit);
398 }
399 }
89504c57 400
e93ec1f7 401 AliDebug(1,Form("MakeClusters: Number of digits %d -> (e %f), ehs %f\n",
402 fDigitsArr->GetEntries(),fMinECut,ehs));
403
85d48773 404 Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0;
405 Int_t iphi=0, ieta=0; // cell eta-phi indexes in SM
406
407 for (Int_t ishiftPhi = 0; ishiftPhi < fMaxShiftPhi; ishiftPhi++){
408 Int_t nClusPhi = (fNPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
89504c57 409
85d48773 410 for (Int_t ishiftEta = 0; ishiftEta < fMaxShiftEta; ishiftEta++) {
89504c57 411
85d48773 412 Int_t nClusEta = (fNEtaDigits - fShiftEta * ishiftEta) / fNeta;
89504c57 413
85d48773 414 Int_t iTotalClus = fNClusters * (ishiftPhi * fMaxShiftEta + ishiftEta);
89504c57 415
e93ec1f7 416 TIter nextdigitC(digitsC);
417 while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
418
89504c57 419 fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
420 fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
421
85d48773 422 Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + fNPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
89504c57 423
85d48773 424 Int_t iTRUphi = iphi_eff / fNPhiDigits;
89504c57 425
85d48773 426 iphi_eff -= iTRUphi * fNPhiDigits;
89504c57 427
27bad386 428 Int_t iClusPhi = iphi_eff / fNphi;
89504c57 429
430 if (iphi_eff < 0 || iClusPhi >= nClusPhi)
431 continue;
432
85d48773 433 Int_t ieta_eff = ieta - fShiftEta * ishiftEta + fNEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
89504c57 434
85d48773 435 Int_t iTRUeta = ieta_eff / fNEtaDigits;
89504c57 436
85d48773 437 ieta_eff -= iTRUeta * fNEtaDigits;
89504c57 438
27bad386 439 Int_t iClusEta = ieta_eff / fNeta;
89504c57 440
441 if (ieta_eff < 0 || iClusEta >= nClusEta)
442 continue;
443
85d48773 444 iphi_eff += iTRUphi * fNPhiDigits;
27bad386 445 iClusPhi = iphi_eff / fNphi;
89504c57 446
85d48773 447 ieta_eff += iTRUeta * fNEtaDigits;
27bad386 448 iClusEta = ieta_eff / fNeta;
89504c57 449
85d48773 450 Int_t iCluster = iClusPhi + iClusEta * fNClusPhiNoShift * fNTRUPhi;
27bad386 451 Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
89504c57 452
85d48773 453 if (iCluster >= fNClusters){
454 AliError(Form("iCluster out of range! iCluster = %d, fNClusters = %d (should never happen...)", iCluster, fNClusters));
89504c57 455 return;
456 }
457
458 iCluster += iTotalClus;
85d48773 459
460 if (iCluster >= fNTotalClus){
461 AliError(Form("iCluster out of range! iCluster = %d, fNTotalClus = %d (should never happen...)", iCluster, fNTotalClus));
462 return;
89504c57 463 }
85d48773 464
465 if (iDigit >= fNDigitsCluster){
466 AliError(Form("iDigit out of range! iDigit = %d, fNDigitsCluster = %d (should never happen...)", iDigit, fNDigitsCluster));
467 return;
468 }
469
470 if (fClustersArray[iCluster][iDigit] != 0){
471 AliError("Digit already added! (should never happen...)");
89504c57 472 return;
473 }
474
27bad386 475 fClustersArray[iCluster][iDigit] = digit;
89504c57 476
477 } // loop on digit
478
479 } // loop on eta shift
480
06ca320f 481 } // loop on phi shift
641122eb 482
85d48773 483 AliEMCALRecPoint *recPoint = 0;
484 Bool_t recPointOk = kFALSE;
485 for (Int_t iCluster = 0; iCluster < fNTotalClus; iCluster++) {
486
487 if (!recPoint) {
488 if(fNumberOfECAClusters >= fRecPoints->GetSize()) fRecPoints->Expand(fNumberOfECAClusters*2+1);
641122eb 489
85d48773 490 (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("");
491 recPoint = static_cast<AliEMCALRecPoint*>(fRecPoints->At(fNumberOfECAClusters));
492 }
89504c57 493
06ca320f 494 if (recPoint) {
06ca320f 495 recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
496 recPoint->SetUniqueID(iCluster);
85d48773 497 fNumberOfECAClusters++;
27bad386 498
85d48773 499 for (Int_t iDigit = 0; iDigit < fNDigitsCluster; iDigit++) {
500 if (fClustersArray[iCluster][iDigit] == 0) continue;
06ca320f 501 digit = fClustersArray[iCluster][iDigit];
e93ec1f7 502 recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
85d48773 503 fClustersArray[iCluster][iDigit] = 0;
504 recPointOk = kTRUE;
505 }
506
507 if (recPointOk) { // unset the pointer so that a new rec point will be allocated in the next iteration
508 recPoint = 0;
509 recPointOk = kFALSE;
06ca320f 510 }
511 }
85d48773 512 else {
513 AliError("Error allocating rec points!");
514 break;
515 }
516 }
517
518 if (!recPointOk) {
519 fRecPoints->RemoveLast();
520 fNumberOfECAClusters--;
06ca320f 521 }
39779e58 522
523 delete digitsC;
06ca320f 524 AliDebug(1, Form("MakeClusters: Number of digits %d -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
06ca320f 525 AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast()));
89504c57 526}