- fixing major bug in clusterizer causing almost all events to contain only one cluster
[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterizer.cxx
1 // $Id$\r
2 \r
3 /**************************************************************************\r
4  * This file is property of and copyright by the ALICE HLT Project        *\r
5  * All rights reserved.                                                   *\r
6  *                                                                        *\r
7  * Primary Authors: Oystein Djuvsland                                     *\r
8  *                                                                        *\r
9  * Permission to use, copy, modify and distribute this software and its   *\r
10  * documentation strictly for non-commercial purposes is hereby granted   *\r
11  * without fee, provided that the above copyright notice appears in all   *\r
12  * copies and that both the copyright notice and this permission notice   *\r
13  * appear in the supporting documentation. The authors make no claims     *\r
14  * about the suitability of this software for any purpose. It is          *\r
15  * provided "as is" without express or implied warranty.                  *\r
16  **************************************************************************/\r
17 \r
18 /**\r
19  * @file   AliHLTCaloClusterizer.cxx\r
20  * @author Oystein Djuvsland\r
21  * @date\r
22  * @brief  Clusterizer for PHOS HLT\r
23  */\r
24 \r
25 // see header file for class documentation\r
26 // or\r
27 // refer to README to build package\r
28 // or\r
29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt\r
30 \r
31 #include "AliHLTCaloClusterizer.h"\r
32 #include "AliHLTLogging.h"\r
33 #include "TMath.h"\r
34 #include "AliHLTCaloRecPointDataStruct.h"\r
35 #include "AliHLTCaloDigitDataStruct.h"\r
36 #include "AliHLTCaloDigitContainerDataStruct.h"\r
37 #include "AliHLTCaloConstantsHandler.h"\r
38 \r
39 ClassImp(AliHLTCaloClusterizer);\r
40 \r
41 AliHLTCaloClusterizer::AliHLTCaloClusterizer(TString det):\r
42         AliHLTCaloConstantsHandler(det),\r
43         fCompareFunction(CompareDigitsByPosition),\r
44         fRecPointArray(0),\r
45         fRecPointDataPtr(0),\r
46         fFirstRecPointPtr(0),\r
47         fArraySize(0),\r
48         fAvailableSize(0),\r
49         fUsedSize(0),\r
50         fNRecPoints(0),\r
51         fDigitIndexPtr(0),\r
52         fEmcClusteringThreshold(0),\r
53         fEmcMinEnergyThreshold(0),\r
54         fEmcTimeGate(0),\r
55         fDigitsInCluster(0),\r
56         fDigitsPointerArray(0),\r
57         fDigitContainerPtr(0),\r
58         fMaxDigitIndexDiff(0),\r
59         fNDigits(0),\r
60         fSortedByPosition(false),\r
61         fSortedByEnergy(false),\r
62         fSortDigits(false),\r
63         fIsEMCAL(false),\r
64         fBuffer(0)\r
65         \r
66 {\r
67     //See header file for documentation\r
68     //fEmcClusteringThreshold = 0.2;\r
69     //fEmcMinEnergyThreshold = 0.03;\r
70 \r
71     fEmcClusteringThreshold = 0.1;\r
72     fEmcMinEnergyThreshold = 0.01;\r
73     fEmcTimeGate = 1.e-6 ;\r
74 \r
75     fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
76 \r
77 \r
78     fArraySize = 10;\r
79     fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
80 \r
81     fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;\r
82     fBuffer = new UChar_t[fAvailableSize]; //FR\r
83     \r
84     fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);\r
85     fRecPointDataPtr = fFirstRecPointPtr;\r
86 \r
87 }//end\r
88 \r
89 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()\r
90 {\r
91     //See header file for documentation\r
92   delete [] fBuffer; //FR\r
93   fBuffer = NULL;    //FR\r
94   fAvailableSize = 0;  //FR\r
95 \r
96   delete [] fRecPointArray;\r
97 }\r
98 \r
99 void\r
100 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
101 {\r
102     // See header file for documentation\r
103     fRecPointDataPtr = recPointDataPtr;\r
104 }\r
105 \r
106 Int_t\r
107 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
108 {\r
109     //see header file for documentation\r
110     Int_t nRecPoints = 0;\r
111     fNRecPoints = 0;\r
112     fUsedSize = 0;\r
113     fNDigits = nDigits;\r
114     fRecPointDataPtr = fFirstRecPointPtr;\r
115 \r
116     // Sort our digits\r
117     SortDigits();\r
118 \r
119     //Clusterization starts\r
120     for (Int_t i = 0; i < nDigits; i++)\r
121     {\r
122         fDigitsInCluster = 0;\r
123 \r
124          HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);\r
125         \r
126         if (fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold && fSortedByEnergy)\r
127         {\r
128            // Since we have sorted by energy the next digit will have even lower energy, so we return \r
129            return fNRecPoints;\r
130         }\r
131 \r
132         if(fDigitsPointerArray[i]->fAssociatedCluster != -1)\r
133         {\r
134            // The digit is added to a previous cluster, continue\r
135            continue;\r
136         }\r
137 \r
138         CheckArray();\r
139         CheckBuffer();\r
140 \r
141         // First digit is placed at the fDigits member variable in the recpoint\r
142         fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
143 \r
144         fRecPointDataPtr->fAmp = 0;\r
145         fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
146 \r
147         // Assigning the digit to this rec point\r
148         fRecPointDataPtr->fDigits = i;\r
149         fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
150 \r
151         // Incrementing the pointer to be ready for new entry\r
152         fDigitIndexPtr++;\r
153 \r
154         fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
155     \r
156         \r
157         //fDigitsPointerArray[i]->fEnergy = 0;\r
158         fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;\r
159         \r
160         \r
161         fDigitsInCluster++;\r
162         nRecPoints++;\r
163 \r
164         // Scanning for the neighbours\r
165         if (ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
166         {\r
167             return -1;\r
168         }\r
169 \r
170         //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);\r
171 \r
172         fRecPointDataPtr->fMultiplicity = fDigitsInCluster;\r
173         fRecPointArray[fNRecPoints] = fRecPointDataPtr;\r
174 \r
175         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
176 \r
177         fNRecPoints++;\r
178 \r
179     }//end of clusterization\r
180 \r
181     return nRecPoints;\r
182 }\r
183 \r
184 Int_t\r
185 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
186 {\r
187     //see header file for documentation\r
188 \r
189     // The following cuts can be used if we sort by posisiton. Not tested, but it should be fine...\r
190     Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
191     Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
192 \r
193     // All digits for now\r
194     max = fNDigits;\r
195     min = 0;\r
196 \r
197     for (Int_t j = min; j < max; j++)\r
198     {\r
199         if (fDigitsPointerArray[j]->fAssociatedCluster == -1 &&  fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
200         {\r
201             if (j != index)\r
202             {\r
203                 if (AreNeighbours(fDigitsPointerArray[index],\r
204                                   fDigitsPointerArray[j]))\r
205                 {\r
206                     // Check that the buffer is large enough for adding a digit (can be heavily improved wrt performance)\r
207                     CheckBuffer();\r
208 \r
209                     // Assigning index to digit\r
210                     *fDigitIndexPtr = j;\r
211                     fUsedSize += sizeof(Int_t);\r
212 \r
213                     // Incrementing digit pointer to be ready for new entry\r
214                     fDigitIndexPtr++;\r
215 \r
216                     // Adding the digit energy to the rec point\r
217                     fRecPointDataPtr->fAmp += fDigitsPointerArray[j]->fEnergy;\r
218 \r
219                     // Setting energy to 0\r
220                     //fDigitsPointerArray[j]->fEnergy = 0;\r
221                     \r
222                     // Setting the associated cluster \r
223                     fDigitsPointerArray[j]->fAssociatedCluster = fNRecPoints;\r
224                     \r
225                     HLTDebug("Added digit with index: %d, energy: %f, to associated cluster: %d", fDigitsPointerArray[j]->fID, fDigitsPointerArray[j]->fEnergy, fDigitsPointerArray[j]->fAssociatedCluster);\r
226                     \r
227                     fDigitsInCluster++;\r
228 \r
229                     // Scan for neighbours of this digit\r
230                     ScanForNeighbourDigits(j, recPoint);\r
231                 }\r
232             }\r
233         }\r
234     }\r
235     return 0;\r
236 }\r
237 \r
238 Int_t\r
239 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,\r
240                                      AliHLTCaloDigitDataStruct* digit2)\r
241 {\r
242     //see header file for documentation\r
243     if ( (digit1->fModule == digit2->fModule) || AreEdgeCells(digit1, digit2))\r
244     {\r
245         Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );\r
246         Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );\r
247 \r
248         // Common edge defines neighbour\r
249         //if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
250           // Common edge and corner defines neighbour\r
251         if (( coldiff <= 1   &&  rowdiff <= 1 ))\r
252         {\r
253             // Check also for time\r
254             if (TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
255             {\r
256                 return 1;\r
257             }\r
258         }\r
259     }\r
260     return 0;\r
261 }\r
262 \r
263 \r
264 \r
265 Int_t AliHLTCaloClusterizer::CheckArray()\r
266 {\r
267     // See header file for class documentation\r
268     if (fArraySize == fNRecPoints)\r
269     {\r
270         fArraySize *= 2;\r
271         AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
272         memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));\r
273         delete [] fRecPointArray;\r
274         fRecPointArray = tmp;\r
275     }\r
276     return 0;\r
277 }\r
278 \r
279 Int_t AliHLTCaloClusterizer::CheckBuffer()\r
280 {\r
281     // See header file for class documentation\r
282     if ((fAvailableSize - fUsedSize) < (Int_t)sizeof(AliHLTCaloRecPointDataStruct))\r
283     {\r
284         Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);\r
285         Int_t digitIndexOffset = reinterpret_cast<UChar_t*>(fDigitIndexPtr) - reinterpret_cast<UChar_t*>(fRecPointDataPtr);\r
286         UChar_t *tmp = new UChar_t[fAvailableSize*2];\r
287         \r
288         if (tmp == NULL)\r
289           {\r
290             HLTError("Pointer error");\r
291             return(-1);\r
292           }\r
293         \r
294         memcpy(tmp, fFirstRecPointPtr, fUsedSize);\r
295         fAvailableSize *= 2;\r
296         for (Int_t n = 0; n < fNRecPoints; n++)\r
297         {\r
298             fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
299         }\r
300         delete [] fBuffer; //FR\r
301         fBuffer = tmp; //FR\r
302         \r
303         fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
304         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);\r
305         fDigitIndexPtr = reinterpret_cast<Int_t*>(reinterpret_cast<UChar_t*>(fRecPointDataPtr) + digitIndexOffset);\r
306         //fUsedSize = 0;\r
307     }\r
308     return 0;\r
309 }\r
310 \r
311 void AliHLTCaloClusterizer::SetSortDigitsByPosition()\r
312 {\r
313     // Sort the digit pointers by position\r
314     fCompareFunction = &CompareDigitsByPosition;\r
315     fSortDigits = true;\r
316     fSortedByPosition = true;\r
317 }\r
318 \r
319 void AliHLTCaloClusterizer::SetSortDigitsByEnergy()\r
320 {\r
321     // See header file for class documentation\r
322     fCompareFunction = &CompareDigitsByEnergy;\r
323     fSortDigits = true;\r
324     fSortedByEnergy = true;\r
325 }\r
326 \r
327 void AliHLTCaloClusterizer::SortDigits()\r
328 {\r
329     // See header file for class documentation\r
330     if (fSortDigits) qsort(fDigitsPointerArray, fNDigits, sizeof(AliHLTCaloDigitDataStruct*), fCompareFunction);\r
331 }\r
332 \r
333 Int_t\r
334 AliHLTCaloClusterizer::CompareDigitsByPosition(const void *dig0, const void *dig1)\r
335 {\r
336     // See header file for documentation\r
337     return (*((AliHLTCaloDigitDataStruct**)(dig0)))->fID - (*((AliHLTCaloDigitDataStruct**)(dig1)))->fID;\r
338 }\r
339 \r
340 Int_t\r
341 AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)\r
342 {\r
343     // See header file for documentation\r
344   if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;\r
345   return 1;\r
346 }\r
347 \r
348 void AliHLTCaloClusterizer::SetDetector(TString det)\r
349 {\r
350   if(det.CompareTo("EMCAL"))\r
351   {\r
352     fIsEMCAL = true;\r
353   }\r
354   else\r
355   {\r
356     fIsEMCAL = false;\r
357   }\r
358 }\r
359 \r
360 Bool_t AliHLTCaloClusterizer::AreEdgeCells(AliHLTCaloDigitDataStruct *digit0, AliHLTCaloDigitDataStruct *digit1)\r
361 {\r
362   if(fIsEMCAL)\r
363   {\r
364     Int_t modDiff = digit0->fModule - digit1->fModule;\r
365     if(TMath::Abs(modDiff) > 1) return kFALSE;\r
366     if(digit0->fModule > digit1->fModule && digit1->fModule%2 == 0) \r
367     {\r
368       if(digit0->fZ == 0 && digit1->fZ == (fCaloConstants->GetNZROWSMOD()-1))\r
369       return kTRUE;\r
370     }\r
371     if(digit1->fModule > digit0->fModule && digit0->fModule%2 == 0) \r
372     {\r
373       if(digit1->fZ == 0 && digit0->fZ == (fCaloConstants->GetNZROWSMOD()-1))\r
374       return kTRUE;\r
375     }\r
376   }\r
377   \r
378   return false;\r
379 \r
380 }\r