Adding a reminder for coders
[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       fClustersArray(0)
51 {
52   // Constructor
53 }
54
55 //__________________________________________________________________________________________
56 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry)
57   : AliEMCALClusterizer(geometry), 
58     fNphi(4), 
59     fNeta(4), 
60     fShiftPhi(2), 
61     fShiftEta(2),
62     fTRUshift(0),
63     fClustersArray(0)
64 {
65   // Constructor
66 }
67
68 //__________________________________________________________________________________________
69 AliEMCALClusterizerFixedWindow::AliEMCALClusterizerFixedWindow(AliEMCALGeometry* geometry, AliEMCALCalibData * calib, AliCaloCalibPedestal * caloped)
70   : AliEMCALClusterizer(geometry, calib, caloped),
71     fNphi(4), 
72     fNeta(4), 
73     fShiftPhi(2), 
74     fShiftEta(2),
75     fTRUshift(0),
76     fClustersArray(0)
77 {
78   // Constructor
79 }
80
81 //__________________________________________________________________________________________
82 AliEMCALClusterizerFixedWindow::~AliEMCALClusterizerFixedWindow()
83 {
84   // Destructor
85   
86   delete fClustersArray;
87 }
88
89 //__________________________________________________________________________________________
90 void AliEMCALClusterizerFixedWindow::SetNphi (Int_t n) 
91 {
92   // Set fNphi; if clusterizer already initialized gives a warning and does nothing
93   
94   if (fClustersArray)
95     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
96   else
97     fNphi = n;
98 }
99
100 //__________________________________________________________________________________________
101 void AliEMCALClusterizerFixedWindow::SetNeta (Int_t n) 
102 {
103   // Set fNeta; if clusterizer already initialized gives a warning and does nothing
104   
105   if (fClustersArray)
106     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
107   else
108     fNeta = n;
109 }
110
111 //__________________________________________________________________________________________
112 void AliEMCALClusterizerFixedWindow::SetShiftPhi (Int_t s) 
113 {
114   // Set fShiftPhi; if clusterizer already initialized gives a warning and does nothing
115   
116   if (fClustersArray)
117     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
118   else
119     fShiftPhi = s;
120 }
121
122 //__________________________________________________________________________________________
123 void AliEMCALClusterizerFixedWindow::SetShiftEta (Int_t s) 
124 {
125   // Set fShiftEta; if clusterizer already initialized gives a warning and does nothing
126   
127   if (fClustersArray)
128     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
129   else
130     fShiftEta = s;
131 }
132
133 //__________________________________________________________________________________________
134 void AliEMCALClusterizerFixedWindow::SetTRUshift(Bool_t b) 
135 {
136   // Set fTRUshift; if clusterizer already initialized gives a warning and does nothing
137   
138   if (fClustersArray)
139     AliWarning("Clusterizer already initialized. Unable to change the parameters.");
140   else
141     fTRUshift = b;
142 }
143
144
145 //__________________________________________________________________________________________
146 void AliEMCALClusterizerFixedWindow::Digits2Clusters(Option_t * option)
147 {
148   // Steering method to perform clusterization for the current event 
149         
150   if (strstr(option,"tim"))
151     gBenchmark->Start("EMCALClusterizer"); 
152         
153   if (strstr(option,"print"))
154     Print(""); 
155         
156   //Get calibration parameters from file or digitizer default values.
157   GetCalibrationParameters();
158         
159   //Get dead channel map from file or digitizer default values.
160   GetCaloCalibPedestal();
161         
162   MakeClusters();  //only the real clusters
163         
164   if (fToUnfold) {
165     fClusterUnfolding->SetInput(fNumberOfECAClusters,fRecPoints,fDigitsArr);
166     fClusterUnfolding->MakeUnfolding();
167   }
168         
169   //Evaluate position, dispersion and other RecPoint properties for EC section 
170   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) { 
171     AliEMCALRecPoint * rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
172     if (rp) {
173       rp->EvalAll(fECAW0,fDigitsArr,fJustClusters);
174       AliDebug(5, Form("MAX INDEX %d ", rp->GetMaximalEnergyIndex()));
175       //For each rec.point set the distance to the nearest bad crystal
176       if (fCaloPed)
177         rp->EvalDistanceToBadChannels(fCaloPed);
178     }
179   }
180   
181   fRecPoints->Sort();
182         
183   for (Int_t index = 0; index < fRecPoints->GetEntries(); index++) {
184     AliEMCALRecPoint *rp = dynamic_cast<AliEMCALRecPoint *>(fRecPoints->At(index));
185     if (rp) {
186       rp->SetIndexInList(index);
187     }
188     else AliFatal("RecPoint NULL!!");
189   }
190         
191   if (fTreeR)
192     fTreeR->Fill();
193         
194   if (strstr(option,"deb") || strstr(option,"all"))  
195     PrintRecPoints(option);
196         
197   AliDebug(1,Form("EMCAL Clusterizer found %d Rec Points",fRecPoints->GetEntriesFast()));
198         
199   if (strstr(option,"tim")) {
200     gBenchmark->Stop("EMCALClusterizer");
201     printf("Exec took %f seconds for Clusterizing", 
202            gBenchmark->GetCpuTime("EMCALClusterizer"));
203   }    
204 }
205
206 //__________________________________________________________________________________________
207 void AliEMCALClusterizerFixedWindow::MakeClusters()
208 {
209   // Make clusters
210         
211   if (fGeom == 0) 
212     AliFatal("Did not get geometry from EMCALLoader");
213         
214   fNumberOfECAClusters = 0;
215   fRecPoints->Delete();
216         
217   Int_t nSupMod=0, nModule=0, nIphi=0, nIeta=0, iphi=0, ieta=0;
218         
219   // Defining geometry and clusterization parameter
220   Int_t nEtaDigitsSupMod = fGeom->GetNEta() * fGeom->GetNETAdiv(); // always 48?;
221   Int_t nPhiDigitsSupMod = fGeom->GetNPhi() * fGeom->GetNPHIdiv(); // always 24?;
222   
223   Int_t nTRUPhi = 1;
224   Int_t nTRUEta = 1;
225   
226   Int_t nEtaDigits = nEtaDigitsSupMod * fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
227   Int_t nPhiDigits = nPhiDigitsSupMod * fGeom->GetNPhiSuperModule();    
228   
229   if (fTRUshift){
230     nTRUPhi = fGeom->GetNPhiSuperModule() * 3;
231     nTRUEta = fGeom->GetNumberOfSuperModules() / fGeom->GetNPhiSuperModule();
232     nEtaDigits /= nTRUEta;
233     nPhiDigits /= nTRUPhi;
234   }
235
236   // Check if clusterizer parameter are compatible with calorimeter geometry
237   if (nEtaDigits < fNeta){
238     AliFatal(Form("Error: fNeta = %d is greater than nEtaDigits = %d.",fNeta,nEtaDigits));
239     return;
240   }
241   if (nPhiDigits < fNphi){
242     AliFatal(Form("Error: fNphi = %d is greater than nPhiDigits = %d.",fNphi,nPhiDigits));
243     return;
244   }
245   if (nEtaDigits % fShiftEta != 0){
246     AliFatal(Form("Error: fShiftEta = %d is such that clusters cannot slide the whole calorimeter (nEtaDigits = %d).",fShiftEta,nEtaDigits));
247     return;
248   }
249   if (nPhiDigits % fShiftPhi != 0){
250     AliFatal(Form("Error: fShiftPhi = %d is such that clusters cannot slide the whole calorimeter (nPhiDigits = %d).",fShiftPhi,nPhiDigits));
251     return;
252   }
253   if (fNeta % fShiftEta != 0){
254     AliFatal(Form("Error: fShiftEta = %d is not divisor of fNeta = %d.",fShiftEta,fNeta));
255     return;
256   }
257   if (fNphi % fShiftPhi != 0){
258     AliFatal(Form("Error: fShiftPhi = %d is not divisor of fNphi = %d).",fShiftPhi,fNphi));
259     return;
260   }
261   
262   Int_t maxiShiftPhi = fNphi / fShiftPhi;
263   Int_t maxiShiftEta = fNeta / fShiftEta;
264         
265   Int_t nDigitsCluster = fNphi * fNeta;
266   
267   Int_t nClusEtaNoShift = nEtaDigits / fNeta;
268   Int_t nClusPhiNoShift = nPhiDigits / fNphi;
269   
270   Int_t nClusters =  nClusEtaNoShift * nClusPhiNoShift * nTRUEta * nTRUPhi;
271   
272   Int_t nTotalClus = nClusters * maxiShiftEta * maxiShiftPhi;
273   
274   if (!fClustersArray){
275     fClustersArray = new AliEMCALDigit**[nTotalClus];
276     for (Int_t i = 0; i < nTotalClus; i++)
277     {
278       fClustersArray[i] = NULL;
279     }
280   }
281   
282   AliEMCALDigit *digit = 0;
283   
284   for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++){
285     Int_t nClusPhi = (nPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
286     
287     for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++){
288       
289       Int_t nClusEta = (nEtaDigits - fShiftEta * ishiftEta) / fNeta; 
290       
291       Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
292       
293       TIter nextdigit(fDigitsArr);
294       
295       nextdigit.Reset();
296       
297       while ((digit = static_cast<AliEMCALDigit*>(nextdigit()))){
298         fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
299         fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
300         
301         Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
302         
303         Int_t iTRUphi = iphi_eff / nPhiDigits;
304         
305         iphi_eff -= iTRUphi * nPhiDigits;
306         
307         Int_t iClusPhi = iphi_eff / fNphi; 
308         
309         if (iphi_eff < 0 || iClusPhi >= nClusPhi) 
310           continue;
311         
312         Int_t ieta_eff = ieta - fShiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
313         
314         Int_t iTRUeta = ieta_eff / nEtaDigits;
315         
316         ieta_eff -= iTRUeta * nEtaDigits;
317         
318         Int_t iClusEta = ieta_eff / fNeta; 
319         
320         if (ieta_eff < 0 || iClusEta >= nClusEta) 
321           continue;
322         
323         iphi_eff += iTRUphi * nPhiDigits;
324         iClusPhi = iphi_eff / fNphi; 
325         
326         ieta_eff += iTRUeta * nEtaDigits;
327         iClusEta = ieta_eff / fNeta; 
328         
329         Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi; 
330         Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
331
332         if (iCluster >= nClusters){
333           AliWarning(Form("iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
334           return;
335         }
336         
337         iCluster += iTotalClus;
338         
339         if (fClustersArray[iCluster] == NULL){
340           fNumberOfECAClusters++;
341           fClustersArray[iCluster] = new AliEMCALDigit*[nDigitsCluster];
342           for (Int_t i = 0; i < nDigitsCluster; i++){
343             fClustersArray[iCluster][i] = NULL;
344           }
345         }
346         
347         if (fClustersArray[iCluster][iDigit] != NULL){
348           AliWarning("Digit already added!");
349           return;
350         }
351         
352         fClustersArray[iCluster][iDigit] = digit;
353         
354       } // loop on digit
355       
356     } // loop on eta shift
357     
358   } // loop on phi shift
359         
360   Int_t iRecPoint = 0;
361   for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++){
362     
363     if (fClustersArray[iCluster] == NULL) continue;
364                 
365     (*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
366     AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
367                 
368     if (recPoint) {
369       iRecPoint++;       
370       recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
371       recPoint->SetUniqueID(iCluster);
372
373       for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
374         if (fClustersArray[iCluster][iDigit] == NULL) continue;
375         digit = fClustersArray[iCluster][iDigit];
376         Float_t dEnergyCalibrated = digit->GetAmplitude();
377         Float_t time              = digit->GetTime();
378         Calibrate(dEnergyCalibrated,time,digit->GetId());
379         digit->SetCalibAmp(dEnergyCalibrated);
380         recPoint->AddDigit(*digit, dEnergyCalibrated, kFALSE); //Time or TimeR?
381         fClustersArray[iCluster][iDigit] = NULL;
382       }
383     }
384     
385     delete[] fClustersArray[iCluster];
386     fClustersArray[iCluster] = NULL;
387   }
388   
389   AliDebug(1, Form("MakeClusters: Number of digits %d  -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
390         
391   AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast())); 
392 }