#95494: Change AliLog message in AliAltroRawStreamV3.cxx (additional reduction of...
[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   // Set up TObjArray with pointers to digits to work on calibrated digits 
283   TObjArray *digitsC = new TObjArray();
284   AliEMCALDigit *digit;
285   Float_t dEnergyCalibrated = 0.0, ehs = 0.0, time = 0.0;
286   TIter nextdigit(fDigitsArr);
287   while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigit()))) { // calibrate and clean up digits
288     dEnergyCalibrated =  digit->GetAmplitude();
289     time              =  digit->GetTime();
290     Calibrate(dEnergyCalibrated, time, digit->GetId());
291     digit->SetCalibAmp(dEnergyCalibrated);
292     digit->SetTime(time);
293     if (dEnergyCalibrated < fMinECut) {
294       continue;
295     }
296     else if (!fGeom->CheckAbsCellId(digit->GetId())) {
297       continue;
298     }
299     else {
300       ehs += dEnergyCalibrated;
301       digitsC->AddLast(digit);
302     }
303   } 
304   
305   AliDebug(1,Form("MakeClusters: Number of digits %d  -> (e %f), ehs %f\n",
306                   fDigitsArr->GetEntries(),fMinECut,ehs));
307    
308   for (Int_t ishiftPhi = 0; ishiftPhi < maxiShiftPhi; ishiftPhi++){
309     Int_t nClusPhi = (nPhiDigits - fShiftPhi * ishiftPhi) / fNphi;
310     
311     for (Int_t ishiftEta = 0; ishiftEta < maxiShiftEta; ishiftEta++) {
312       
313       Int_t nClusEta = (nEtaDigits - fShiftEta * ishiftEta) / fNeta; 
314       
315       Int_t iTotalClus = nClusters * (ishiftPhi * maxiShiftEta + ishiftEta);
316       
317       TIter nextdigitC(digitsC);
318       while ((digit = dynamic_cast<AliEMCALDigit*>(nextdigitC()))) { // scan over the list of digitsC
319         
320         fGeom->GetCellIndex (digit->GetId(), nSupMod, nModule, nIphi, nIeta);
321         fGeom->GetCellPhiEtaIndexInSModule (nSupMod, nModule, nIphi, nIeta, iphi, ieta);
322         
323         Int_t iphi_eff = iphi - fShiftPhi * ishiftPhi + nPhiDigitsSupMod * (nSupMod / 2); // N supermodules along phi
324         
325         Int_t iTRUphi = iphi_eff / nPhiDigits;
326         
327         iphi_eff -= iTRUphi * nPhiDigits;
328         
329         Int_t iClusPhi = iphi_eff / fNphi; 
330         
331         if (iphi_eff < 0 || iClusPhi >= nClusPhi) 
332           continue;
333         
334         Int_t ieta_eff = ieta - fShiftEta * ishiftEta + nEtaDigitsSupMod * (nSupMod % 2); // 2 supermodules along eta
335         
336         Int_t iTRUeta = ieta_eff / nEtaDigits;
337         
338         ieta_eff -= iTRUeta * nEtaDigits;
339         
340         Int_t iClusEta = ieta_eff / fNeta; 
341         
342         if (ieta_eff < 0 || iClusEta >= nClusEta) 
343           continue;
344         
345         iphi_eff += iTRUphi * nPhiDigits;
346         iClusPhi = iphi_eff / fNphi; 
347         
348         ieta_eff += iTRUeta * nEtaDigits;
349         iClusEta = ieta_eff / fNeta; 
350         
351         Int_t iCluster = iClusPhi + iClusEta * nClusPhiNoShift * nTRUPhi; 
352         Int_t iDigit = iphi_eff % fNphi + (ieta_eff % fNeta) * fNphi;
353
354         if (iCluster >= nClusters){
355           AliWarning(Form("iCluster out of range! iCluster = %d, nClusters = %d", iCluster, nClusters));
356           return;
357         }
358         
359         iCluster += iTotalClus;
360         
361         if (fClustersArray[iCluster] == NULL){
362           fNumberOfECAClusters++;
363           fClustersArray[iCluster] = new AliEMCALDigit*[nDigitsCluster];
364           for (Int_t i = 0; i < nDigitsCluster; i++){
365             fClustersArray[iCluster][i] = NULL;
366           }
367         }
368         
369         if (fClustersArray[iCluster][iDigit] != NULL){
370           AliWarning("Digit already added!");
371           return;
372         }
373         
374         fClustersArray[iCluster][iDigit] = digit;
375         
376       } // loop on digit
377       
378     } // loop on eta shift
379     
380   } // loop on phi shift
381         
382   if(fNumberOfECAClusters >= fRecPoints->GetSize()) 
383     fRecPoints->Expand(fNumberOfECAClusters+1);
384   
385   Int_t iRecPoint = 0;
386   for (Int_t iCluster = 0; iCluster < nTotalClus; iCluster++) {
387     
388     if (fClustersArray[iCluster] == NULL) 
389       continue;
390
391     (*fRecPoints)[iRecPoint] = new AliEMCALRecPoint("");
392     AliEMCALRecPoint *recPoint = dynamic_cast<AliEMCALRecPoint*> (fRecPoints->At(iRecPoint));
393                 
394     if (recPoint) {
395       iRecPoint++;       
396       recPoint->SetClusterType(AliVCluster::kEMCALClusterv1);
397       recPoint->SetUniqueID(iCluster);
398
399       for (Int_t iDigit = 0; iDigit < nDigitsCluster; iDigit++){
400         if (fClustersArray[iCluster][iDigit] == NULL) continue;
401         digit = fClustersArray[iCluster][iDigit];
402         recPoint->AddDigit(*digit, digit->GetCalibAmp(), kFALSE); //Time or TimeR?
403         fClustersArray[iCluster][iDigit] = NULL;
404       }
405     }
406     
407     delete[] fClustersArray[iCluster];
408     fClustersArray[iCluster] = NULL;
409   }
410
411   delete digitsC;
412   AliDebug(1, Form("MakeClusters: Number of digits %d  -> (e %f)\n", fDigitsArr->GetEntries(),fMinECut));
413   AliDebug(1, Form("total no of clusters %d from %d digits", fNumberOfECAClusters, fDigitsArr->GetEntriesFast())); 
414 }