put back previous default value, although decission needs to be made to which one...
[u/mrichter/AliRoot.git] / EMCAL / AliEMCALClusterizerFixedWindow.cxx
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 "AliESDCaloCluster.h"
36 #include "AliEMCALUnfolding.h"
37
38 #include "AliEMCALClusterizerFixedWindow.h"
39
40 ClassImp(AliEMCALClusterizerFixedWindow)
41
42 //__________________________________________________________________________________________
43 AliEMCALClusterizerFixedWindow::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)
65 {
66   // Constructor
67 }
68
69 //__________________________________________________________________________________________
70 AliEMCALClusterizerFixedWindow::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)
92 {
93   // Constructor
94 }
95
96 //__________________________________________________________________________________________
97 AliEMCALClusterizerFixedWindow::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)
119 {
120   // Constructor
121 }
122
123 //__________________________________________________________________________________________
124 AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
125 {
126   // Destructor
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     }
135     delete[] fClustersArray;
136     fClustersArray = 0;
137   }
138
139 }
140
141 //__________________________________________________________________________________________
142 void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n) 
143 {
144   // Set fNphi; if clusterizer already initialized gives a warning and does nothing
145   
146   if (fInitialized != 0)
147     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
148   else
149     fNphi = n;
150 }
151
152 //__________________________________________________________________________________________
153 void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n) 
154 {
155   // Set fNeta; if clusterizer already initialized gives a warning and does nothing
156   
157   if (fInitialized != 0)
158     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
159   else
160     fNeta = n;
161 }
162
163 //__________________________________________________________________________________________
164 void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s) 
165 {
166   // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
167   
168   if (fInitialized != 0)
169     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
170   else
171     fShiftPhi = s;
172 }
173
174 //__________________________________________________________________________________________
175 void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s) 
176 {
177   // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
178   
179   if (fInitialized != 0)
180     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
181   else
182     fShiftEta = s;
183 }
184
185 //__________________________________________________________________________________________
186 void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b) 
187 {
188   // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
189   
190   if (fInitialized != 0)
191     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
192   else
193     fTRUshift = b;
194 }
195
196
197 //__________________________________________________________________________________________
198 void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
199 {
200   // Steering method to perform clusterization for the current event 
201  
202   static Float_t cputime = 0;
203   static Float_t realtime = 0;
204
205   if (strstr(option,"tim"))
206     gBenchmark->Start("EMCALClusterizer"); 
207         
208   if (strstr(option,"print"))
209     Print(""); 
210         
211   //Get calibration parameters from file or digitizer default values.
212   GetCalibrationParameters();
213         
214   //Get dead channel map from file or digitizer default values.
215   GetCaloCalibPedestal();
216         
217   MakeClusters();  //only the real clusters
218         
219   if (fToUnfold) {
220     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
221     fClusterUnfolding->MakeUnfolding();
222   }
223         
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   }
235   
236   fRecPoints->Sort();
237         
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   }
245         
246   if (fTreeR)
247     fTreeR->Fill();
248         
249   if (strstr(option,"deb") || strstr(option,"all"))  
250     PrintRecPoints(option);
251         
252   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
253         
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");
260   }    
261 }
262
263 //__________________________________________________________________________________________
264 void AliEMCALClusterizerFixedWindow::ExecOnce()
265 {
266   // Initialize clusterizer.
267   
268   fInitialized = -1;
269
270   if (!fGeom) {
271     AliError("Did not get geometry!");
272     return;
273   }
274         
275   // Defining geometry and clusterization parameter
276   fNEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
277   fNPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
278   
279   fNTRUPhi = 1;
280   fNTRUEta = 1;
281   
282   fNEtaDigits = fNEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
283   fNPhiDigits = fNPhiDigitsSupMod * fGeom->GetNPhiSuperModule();    
284   
285   if (fTRUshift){
286     fNTRUPhi = fGeom->GetNPhiSuperModule() * 3;
287     fNTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
288     fNEtaDigits /= fNTRUEta;
289     fNPhiDigits /= fNTRUPhi;
290   }
291
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));
295     return;
296   }
297   if (fNPhiDigits < fNphi){
298     AliError(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,fNPhiDigits));
299     return;
300   }
301   if (fNEtaDigits % fShiftEta != 0){
302     AliError(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,fNEtaDigits));
303     return;
304   }
305   if (fNPhiDigits % fShiftPhi != 0){
306     AliError(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,fNPhiDigits));
307     return;
308   }
309   if (fNeta % fShiftEta != 0){
310     AliError(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
311     return;
312   }
313   if (fNphi % fShiftPhi != 0){
314     AliError(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
315     return;
316   }
317   
318   fMaxShiftPhi = fNphi / fShiftPhi;
319   fMaxShiftEta = fNeta / fShiftEta;
320   
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     }
337     delete[] fClustersArray;
338     fClustersArray = 0;
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;
346     }
347   }
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 //__________________________________________________________________________________________
363 void 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   }
377   
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) {
390       continue;
391     }
392     else if (!fGeom->CheckAbsCellId(digit->GetId())) {
393       continue;
394     }
395     else {
396       ehs += dEnergyCalibrated;
397       digitsC->AddLast(digit);
398     }
399   } 
400   
401   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
402                   fDigitsArr->GetEntries(),fMinECut,ehs));
403    
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;
409     
410     for (Int_t ishiftEta = 0; ishiftEta < fMaxShiftEta; ishiftEta++) {
411       
412       Int_t nClusEta = (fNEtaDigits - fShiftEta * ishiftEta) / fNeta; 
413       
414       Int_t iTotalClus = fNClusters * (ishiftPhi * fMaxShiftEta + ishiftEta);
415       
416       TIter nextdigitC(digitsC);
417       while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
418         
419         fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
420         fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
421         
422         Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + fNPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
423         
424         Int_t iTRUphi = iphi_eff / fNPhiDigits;
425         
426         iphi_eff -= iTRUphi * fNPhiDigits;
427         
428         Int_t iClusPhi = iphi_eff / fNphi; 
429         
430         if (iphi_eff < 0 || iClusPhi >= nClusPhi) 
431           continue;
432         
433         Int_t ieta_eff = ieta - fShiftEta * ishiftEta + fNEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
434         
435         Int_t iTRUeta = ieta_eff / fNEtaDigits;
436         
437         ieta_eff -= iTRUeta * fNEtaDigits;
438         
439         Int_t iClusEta = ieta_eff / fNeta; 
440         
441         if (ieta_eff < 0 || iClusEta >= nClusEta) 
442           continue;
443         
444         iphi_eff += iTRUphi * fNPhiDigits;
445         iClusPhi = iphi_eff / fNphi; 
446         
447         ieta_eff += iTRUeta * fNEtaDigits;
448         iClusEta = ieta_eff / fNeta; 
449         
450         Int_t iCluster = iClusPhi + iClusEta * fNClusPhiNoShift * fNTRUPhi; 
451         Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
452
453         if (iCluster >= fNClusters){
454           AliError(Form("iCluster out of range! iCluster = %d, fNClusters = %d (should never happen...)", iCluster, fNClusters));
455           return;
456         }
457         
458         iCluster += iTotalClus;
459
460         if (iCluster >= fNTotalClus){
461           AliError(Form("iCluster out of range! iCluster = %d, fNTotalClus = %d (should never happen...)", iCluster, fNTotalClus));
462           return;
463         }
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...)");
472           return;
473         }
474         
475         fClustersArray[iCluster][iDigit] = digit;
476         
477       } // loop on digit
478       
479     } // loop on eta shift
480     
481   } // loop on phi shift
482   
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);
489
490       (*fRecPoints)[fNumberOfECAClusters] = new AliEMCALRecPoint("");
491       recPoint = static_cast<AliEMCALRecPoint*>(fRecPoints->At(fNumberOfECAClusters));
492     }
493                 
494     if (recPoint) {
495       recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
496       recPoint->SetUniqueID(iCluster);
497       fNumberOfECAClusters++;
498
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;
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;
510       }
511     }
512     else {
513       AliError("Error allocating rec points!");
514       break;
515     }
516   }
517
518   if (!recPointOk) {
519     fRecPoints->RemoveLast();
520     fNumberOfECAClusters--;
521   }
522
523   delete digitsC;
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())); 
526 }