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