b5f5d471ce3066d51e20e3ab8ba3b486012be45f
[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         fBuffer(0)\r
64 {\r
65     //See header file for documentation\r
66     //fEmcClusteringThreshold = 0.2;\r
67     //fEmcMinEnergyThreshold = 0.03;\r
68 \r
69     fEmcClusteringThreshold = 0.1;\r
70     fEmcMinEnergyThreshold = 0.01;\r
71     fEmcTimeGate = 1.e-6 ;\r
72 \r
73     fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
74 \r
75 \r
76     fArraySize = 10;\r
77     fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
78 \r
79     fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;\r
80     fBuffer = new UChar_t[fAvailableSize]; //FR\r
81     \r
82     fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);\r
83     fRecPointDataPtr = fFirstRecPointPtr;\r
84 \r
85 }//end\r
86 \r
87 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()\r
88 {\r
89     //See header file for documentation\r
90   delete [] fBuffer; //FR\r
91   fBuffer = NULL;    //FR\r
92   fAvailableSize = 0;  //FR\r
93 \r
94   delete [] fRecPointArray;\r
95 }\r
96 \r
97 void\r
98 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
99 {\r
100     // See header file for documentation\r
101     fRecPointDataPtr = recPointDataPtr;\r
102 }\r
103 \r
104 Int_t\r
105 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
106 {\r
107     //see header file for documentation\r
108     Int_t nRecPoints = 0;\r
109     fNRecPoints = 0;\r
110     fUsedSize = 0;\r
111     fNDigits = nDigits;\r
112     fRecPointDataPtr = fFirstRecPointPtr;\r
113 \r
114     // Sort our digits\r
115     SortDigits();\r
116 \r
117     //Clusterization starts\r
118     for (Int_t i = 0; i < nDigits; i++)\r
119     {\r
120         fDigitsInCluster = 0;\r
121 \r
122          HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);\r
123         \r
124         if (fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold && fSortedByEnergy)\r
125         {\r
126            // Since we have sorted by energy the next digit will have even lower energy, so we return \r
127            return fNRecPoints;\r
128         }\r
129 \r
130         if(fDigitsPointerArray[i]->fAssociatedCluster != -1)\r
131         {\r
132            // The digit is added to a previous cluster, continue\r
133            continue;\r
134         }\r
135 \r
136         CheckArray();\r
137         CheckBuffer();\r
138 \r
139         // First digit is placed at the fDigits member variable in the recpoint\r
140         fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
141 \r
142         fRecPointDataPtr->fAmp = 0;\r
143         fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
144 \r
145         // Assigning the digit to this rec point\r
146         fRecPointDataPtr->fDigits = i;\r
147         fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
148 \r
149         // Incrementing the pointer to be ready for new entry\r
150         fDigitIndexPtr++;\r
151 \r
152         fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
153     \r
154         \r
155         //fDigitsPointerArray[i]->fEnergy = 0;\r
156         fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;\r
157         \r
158         \r
159         fDigitsInCluster++;\r
160         nRecPoints++;\r
161 \r
162         // Scanning for the neighbours\r
163         if (ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
164         {\r
165             return -1;\r
166         }\r
167 \r
168         //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);\r
169 \r
170         fRecPointDataPtr->fMultiplicity = fDigitsInCluster;\r
171         fRecPointArray[fNRecPoints] = fRecPointDataPtr;\r
172 \r
173         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
174 \r
175         fNRecPoints++;\r
176 \r
177     }//end of clusterization\r
178 \r
179     return nRecPoints;\r
180 }\r
181 \r
182 Int_t\r
183 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
184 {\r
185     //see header file for documentation\r
186 \r
187     // The following cuts can be used if we sort by posisiton. Not tested, but it should be fine...\r
188     Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
189     Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
190 \r
191     // All digits for now\r
192     max = fNDigits;\r
193     min = 0;\r
194 \r
195     for (Int_t j = min; j < max; j++)\r
196     {\r
197         if (fDigitsPointerArray[j]->fAssociatedCluster == -1 &&  fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
198         {\r
199             if (j != index)\r
200             {\r
201                 if (AreNeighbours(fDigitsPointerArray[index],\r
202                                   fDigitsPointerArray[j]))\r
203                 {\r
204                     // Check that the buffer is large enough for adding a digit (can be heavily improved wrt performance)\r
205                     CheckBuffer();\r
206 \r
207                     // Assigning index to digit\r
208                     *fDigitIndexPtr = j;\r
209                     fUsedSize += sizeof(Int_t);\r
210 \r
211                     // Incrementing digit pointer to be ready for new entry\r
212                     fDigitIndexPtr++;\r
213 \r
214                     // Adding the digit energy to the rec point\r
215                     fRecPointDataPtr->fAmp += fDigitsPointerArray[j]->fEnergy;\r
216 \r
217                     // Setting energy to 0\r
218                     //fDigitsPointerArray[j]->fEnergy = 0;\r
219                     \r
220                     // Setting the associated cluster \r
221                     fDigitsPointerArray[j]->fAssociatedCluster = fNRecPoints;\r
222                     \r
223                     HLTDebug("Added digit with index: %d, energy: %f, to associated cluster: %d", fDigitsPointerArray[j]->fID, fDigitsPointerArray[j]->fEnergy, fDigitsPointerArray[j]->fAssociatedCluster);\r
224                     \r
225                     fDigitsInCluster++;\r
226 \r
227                     // Scan for neighbours of this digit\r
228                     ScanForNeighbourDigits(j, recPoint);\r
229                 }\r
230             }\r
231         }\r
232     }\r
233     return 0;\r
234 }\r
235 \r
236 Int_t\r
237 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,\r
238                                      AliHLTCaloDigitDataStruct* digit2)\r
239 {\r
240     //see header file for documentation\r
241     if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
242     {\r
243         Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );\r
244         Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );\r
245 \r
246         // As in the offline code we define neighbours as cells that share an edge, a corner is not  enough\r
247         //        if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
248         if (( coldiff <= 1) || ( rowdiff <= 1 ))\r
249         {\r
250             // Check also for time\r
251             if (TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
252             {\r
253                 return 1;\r
254             }\r
255         }\r
256     }\r
257     return 0;\r
258 }\r
259 \r
260 \r
261 \r
262 Int_t AliHLTCaloClusterizer::CheckArray()\r
263 {\r
264     // See header file for class documentation\r
265     if (fArraySize == fNRecPoints)\r
266     {\r
267         fArraySize *= 2;\r
268         AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
269         memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));\r
270         delete [] fRecPointArray;\r
271         fRecPointArray = tmp;\r
272     }\r
273     return 0;\r
274 }\r
275 \r
276 Int_t AliHLTCaloClusterizer::CheckBuffer()\r
277 {\r
278     // See header file for class documentation\r
279     if ((fAvailableSize - fUsedSize) < (Int_t)sizeof(AliHLTCaloRecPointDataStruct))\r
280     {\r
281         Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);\r
282         Int_t digitIndexOffset = reinterpret_cast<UChar_t*>(fDigitIndexPtr) - reinterpret_cast<UChar_t*>(fRecPointDataPtr);\r
283         UChar_t *tmp = new UChar_t[fAvailableSize*2];\r
284         \r
285         if (tmp == NULL)\r
286           {\r
287             HLTError("Pointer error");\r
288             return(-1);\r
289           }\r
290         \r
291         memcpy(tmp, fFirstRecPointPtr, fUsedSize);\r
292         fAvailableSize *= 2;\r
293         for (Int_t n = 0; n < fNRecPoints; n++)\r
294         {\r
295             fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
296         }\r
297         delete [] fBuffer; //FR\r
298         fBuffer = tmp; //FR\r
299         \r
300         fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
301         fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);\r
302         fDigitIndexPtr = reinterpret_cast<Int_t*>(reinterpret_cast<UChar_t*>(fRecPointDataPtr) + digitIndexOffset);\r
303         //fUsedSize = 0;\r
304     }\r
305     return 0;\r
306 }\r
307 \r
308 void AliHLTCaloClusterizer::SetSortDigitsByPosition()\r
309 {\r
310     // Sort the digit pointers by position\r
311     fCompareFunction = &CompareDigitsByPosition;\r
312     fSortDigits = true;\r
313     fSortedByPosition = true;\r
314 }\r
315 \r
316 void AliHLTCaloClusterizer::SetSortDigitsByEnergy()\r
317 {\r
318     // See header file for class documentation\r
319     fCompareFunction = &CompareDigitsByEnergy;\r
320     fSortDigits = true;\r
321     fSortedByEnergy = true;\r
322 }\r
323 \r
324 void AliHLTCaloClusterizer::SortDigits()\r
325 {\r
326     // See header file for class documentation\r
327     if (fSortDigits) qsort(fDigitsPointerArray, fNDigits, sizeof(AliHLTCaloDigitDataStruct*), fCompareFunction);\r
328 }\r
329 \r
330 Int_t\r
331 AliHLTCaloClusterizer::CompareDigitsByPosition(const void *dig0, const void *dig1)\r
332 {\r
333     // See header file for documentation\r
334     return (*((AliHLTCaloDigitDataStruct**)(dig0)))->fID - (*((AliHLTCaloDigitDataStruct**)(dig1)))->fID;\r
335 }\r
336 \r
337 Int_t\r
338 AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)\r
339 {\r
340     // See header file for documentation\r
341   if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;\r
342   return 1;\r
343 }\r