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