- fixing major bug in clusterizer causing almost all events to contain only one cluster
[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterizer.cxx
CommitLineData
ef44ec64 1// $Id$\r
2\r
3/**************************************************************************\r
31b89da4 4 * This file is property of and copyright by the ALICE HLT Project *\r
ef44ec64 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
31b89da4 14 * about the suitability of this software for any purpose. It is *\r
ef44ec64 15 * provided "as is" without express or implied warranty. *\r
16 **************************************************************************/\r
17\r
31b89da4 18/**\r
ef44ec64 19 * @file AliHLTCaloClusterizer.cxx\r
20 * @author Oystein Djuvsland\r
31b89da4 21 * @date\r
22 * @brief Clusterizer for PHOS HLT\r
ef44ec64 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
ef44ec64 32#include "AliHLTLogging.h"\r
33#include "TMath.h"\r
ef44ec64 34#include "AliHLTCaloRecPointDataStruct.h"\r
35#include "AliHLTCaloDigitDataStruct.h"\r
36#include "AliHLTCaloDigitContainerDataStruct.h"\r
4f4b7ba4 37#include "AliHLTCaloConstantsHandler.h"\r
ef44ec64 38\r
39ClassImp(AliHLTCaloClusterizer);\r
40\r
41AliHLTCaloClusterizer::AliHLTCaloClusterizer(TString det):\r
31b89da4 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
7cc065a2 62 fSortDigits(false),\r
36b00077 63 fIsEMCAL(false),\r
7cc065a2 64 fBuffer(0)\r
36b00077 65 \r
ef44ec64 66{\r
31b89da4 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
7cc065a2 82 fBuffer = new UChar_t[fAvailableSize]; //FR\r
dd71d443 83 \r
7cc065a2 84 fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);\r
31b89da4 85 fRecPointDataPtr = fFirstRecPointPtr;\r
c375e15d 86\r
c375e15d 87}//end\r
88\r
31b89da4 89AliHLTCaloClusterizer::~AliHLTCaloClusterizer()\r
ef44ec64 90{\r
31b89da4 91 //See header file for documentation\r
dd71d443 92 delete [] fBuffer; //FR\r
93 fBuffer = NULL; //FR\r
94 fAvailableSize = 0; //FR\r
95\r
9991396c 96 delete [] fRecPointArray;\r
ef44ec64 97}\r
98\r
31b89da4 99void\r
ef44ec64 100AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
101{\r
31b89da4 102 // See header file for documentation\r
103 fRecPointDataPtr = recPointDataPtr;\r
ef44ec64 104}\r
105\r
31b89da4 106Int_t\r
7c80a370 107AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
ef44ec64 108{\r
31b89da4 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
f92543f5 124 HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);\r
125 \r
126 if (fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold && fSortedByEnergy)\r
31b89da4 127 {\r
f92543f5 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
31b89da4 137\r
f92543f5 138 CheckArray();\r
31b89da4 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
f92543f5 155 \r
156 \r
157 //fDigitsPointerArray[i]->fEnergy = 0;\r
158 fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;\r
159 \r
160 \r
161 fDigitsInCluster++;\r
31b89da4 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
ef44ec64 179 }//end of clusterization\r
70c86cde 180\r
181 return nRecPoints;\r
ef44ec64 182}\r
183\r
9bd6f29b 184Int_t\r
ef44ec64 185AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
186{\r
31b89da4 187 //see header file for documentation\r
ef44ec64 188\r
31b89da4 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
ef44ec64 198 {\r
f92543f5 199 if (fDigitsPointerArray[j]->fAssociatedCluster == -1 && fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
31b89da4 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
f92543f5 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
31b89da4 227 fDigitsInCluster++;\r
228\r
229 // Scan for neighbours of this digit\r
230 ScanForNeighbourDigits(j, recPoint);\r
231 }\r
232 }\r
233 }\r
ef44ec64 234 }\r
31b89da4 235 return 0;\r
ef44ec64 236}\r
237\r
31b89da4 238Int_t\r
239AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,\r
240 AliHLTCaloDigitDataStruct* digit2)\r
ef44ec64 241{\r
31b89da4 242 //see header file for documentation\r
36b00077 243 if ( (digit1->fModule == digit2->fModule) || AreEdgeCells(digit1, digit2))\r
31b89da4 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
36b00077 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
31b89da4 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
ef44ec64 259 }\r
31b89da4 260 return 0;\r
ef44ec64 261}\r
7c80a370 262\r
263\r
264\r
265Int_t AliHLTCaloClusterizer::CheckArray()\r
266{\r
31b89da4 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
7c80a370 277}\r
278\r
279Int_t AliHLTCaloClusterizer::CheckBuffer()\r
280{\r
31b89da4 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
dd71d443 287 \r
288 if (tmp == NULL)\r
289 {\r
290 HLTError("Pointer error");\r
291 return(-1);\r
292 }\r
293 \r
31b89da4 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
dd71d443 300 delete [] fBuffer; //FR\r
301 fBuffer = tmp; //FR\r
302 \r
31b89da4 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
7c80a370 309}\r
310\r
31b89da4 311void AliHLTCaloClusterizer::SetSortDigitsByPosition()\r
7c80a370 312{\r
31b89da4 313 // Sort the digit pointers by position\r
314 fCompareFunction = &CompareDigitsByPosition;\r
315 fSortDigits = true;\r
316 fSortedByPosition = true;\r
7c80a370 317}\r
318\r
31b89da4 319void 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
327void 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
333Int_t\r
334AliHLTCaloClusterizer::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
340Int_t\r
341AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)\r
7c80a370 342{\r
31b89da4 343 // See header file for documentation\r
02f2218d 344 if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;\r
345 return 1;\r
2a24cbbe 346}\r
36b00077 347\r
348void 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
360Bool_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