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