Line data Source code
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 :
39 6 : ClassImp(AliHLTCaloClusterizer);
40 :
41 0 : AliHLTCaloClusterizer::AliHLTCaloClusterizer(TString det):
42 0 : AliHLTCaloConstantsHandler(det),
43 0 : fCompareFunction(CompareDigitsByPosition),
44 0 : fRecPointArray(0),
45 0 : fRecPointDataPtr(0),
46 0 : fFirstRecPointPtr(0),
47 0 : fArraySize(0),
48 0 : fAvailableSize(0),
49 0 : fUsedSize(0),
50 0 : fNRecPoints(0),
51 0 : fDigitIndexPtr(0),
52 0 : fEmcClusteringThreshold(0),
53 0 : fEmcMinEnergyThreshold(0),
54 0 : fEmcTimeGate(0),
55 0 : fDigitsInCluster(0),
56 0 : fDigitsPointerArray(0),
57 0 : fDigitContainerPtr(0),
58 0 : fMaxDigitIndexDiff(0),
59 0 : fNDigits(0),
60 0 : fSortedByPosition(false),
61 0 : fSortedByEnergy(false),
62 0 : fSortDigits(false),
63 0 : fIsEMCAL(false),
64 0 : fBuffer(0)
65 :
66 0 : {
67 : //See header file for documentation
68 : //fEmcClusteringThreshold = 0.2;
69 : //fEmcMinEnergyThreshold = 0.03;
70 :
71 0 : fEmcClusteringThreshold = 0.1;
72 0 : fEmcMinEnergyThreshold = 0.01;
73 0 : fEmcTimeGate = 1.e-6 ;
74 :
75 0 : fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();
76 :
77 :
78 0 : fArraySize = 10;
79 0 : fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];
80 :
81 0 : fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;
82 0 : fBuffer = new UChar_t[fAvailableSize]; //FR
83 :
84 0 : fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fBuffer);
85 0 : fRecPointDataPtr = fFirstRecPointPtr;
86 :
87 0 : }//end
88 :
89 0 : AliHLTCaloClusterizer::~AliHLTCaloClusterizer()
90 0 : {
91 : //See header file for documentation
92 0 : delete [] fBuffer; //FR
93 0 : fBuffer = NULL; //FR
94 0 : fAvailableSize = 0; //FR
95 :
96 0 : delete [] fRecPointArray;
97 0 : }
98 :
99 : void
100 : AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)
101 : {
102 : // See header file for documentation
103 0 : fRecPointDataPtr = recPointDataPtr;
104 0 : }
105 :
106 : Int_t
107 : AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)
108 : {
109 : //see header file for documentation
110 : Int_t nRecPoints = 0;
111 0 : fNRecPoints = 0;
112 0 : fUsedSize = 0;
113 0 : fNDigits = nDigits;
114 0 : fRecPointDataPtr = fFirstRecPointPtr;
115 :
116 : // Sort our digits
117 0 : SortDigits();
118 :
119 : //Clusterization starts
120 0 : for (Int_t i = 0; i < nDigits; i++)
121 : {
122 0 : fDigitsInCluster = 0;
123 :
124 : HLTDebug("Digit with energy: %f", fDigitsPointerArray[i]->fEnergy);
125 :
126 0 : 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 0 : return fNRecPoints;
130 : }
131 :
132 0 : if(fDigitsPointerArray[i]->fAssociatedCluster != -1)
133 : {
134 : // The digit is added to a previous cluster, continue
135 : continue;
136 : }
137 :
138 0 : CheckArray();
139 0 : CheckBuffer();
140 :
141 : // First digit is placed at the fDigits member variable in the recpoint
142 0 : fDigitIndexPtr = &(fRecPointDataPtr->fDigits);
143 :
144 0 : fRecPointDataPtr->fAmp = 0;
145 0 : fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;
146 :
147 : // Assigning the digit to this rec point
148 0 : fRecPointDataPtr->fDigits = i;
149 0 : fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);
150 :
151 : // Incrementing the pointer to be ready for new entry
152 0 : fDigitIndexPtr++;
153 :
154 0 : fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;
155 :
156 :
157 : //fDigitsPointerArray[i]->fEnergy = 0;
158 0 : fDigitsPointerArray[i]->fAssociatedCluster = fNRecPoints;
159 :
160 :
161 0 : fDigitsInCluster++;
162 0 : nRecPoints++;
163 :
164 : // Scanning for the neighbours
165 0 : if (ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)
166 : {
167 0 : return -1;
168 : }
169 :
170 : //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);
171 :
172 0 : fRecPointDataPtr->fMultiplicity = fDigitsInCluster;
173 0 : fRecPointArray[fNRecPoints] = fRecPointDataPtr;
174 :
175 0 : fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);
176 :
177 0 : fNRecPoints++;
178 :
179 0 : }//end of clusterization
180 :
181 0 : return nRecPoints;
182 0 : }
183 :
184 : Int_t
185 : AliHLTCaloClusterizer::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 0 : Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);
191 0 : Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));
192 :
193 : // All digits for now
194 0 : max = fNDigits;
195 : min = 0;
196 :
197 0 : for (Int_t j = min; j < max; j++)
198 : {
199 0 : if (fDigitsPointerArray[j]->fAssociatedCluster == -1 && fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)
200 : {
201 0 : if (j != index)
202 : {
203 0 : 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 0 : CheckBuffer();
208 :
209 : // Assigning index to digit
210 0 : *fDigitIndexPtr = j;
211 0 : fUsedSize += sizeof(Int_t);
212 :
213 : // Incrementing digit pointer to be ready for new entry
214 0 : fDigitIndexPtr++;
215 :
216 : // Adding the digit energy to the rec point
217 0 : fRecPointDataPtr->fAmp += fDigitsPointerArray[j]->fEnergy;
218 :
219 : // Setting energy to 0
220 : //fDigitsPointerArray[j]->fEnergy = 0;
221 :
222 : // Setting the associated cluster
223 0 : 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 0 : fDigitsInCluster++;
228 :
229 : // Scan for neighbours of this digit
230 0 : ScanForNeighbourDigits(j, recPoint);
231 0 : }
232 : }
233 : }
234 : }
235 0 : return 0;
236 : }
237 :
238 : Int_t
239 : AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1,
240 : AliHLTCaloDigitDataStruct* digit2)
241 : {
242 : //see header file for documentation
243 0 : if ( (digit1->fModule == digit2->fModule) || AreEdgeCells(digit1, digit2))
244 : {
245 0 : Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
246 0 : 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 0 : if (( coldiff <= 1 && rowdiff <= 1 ))
252 : {
253 : // Check also for time
254 0 : if (TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
255 : {
256 0 : return 1;
257 : }
258 : }
259 0 : }
260 0 : return 0;
261 0 : }
262 :
263 :
264 :
265 : Int_t AliHLTCaloClusterizer::CheckArray()
266 : {
267 : // See header file for class documentation
268 0 : if (fArraySize == fNRecPoints)
269 : {
270 0 : fArraySize *= 2;
271 0 : AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];
272 0 : memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));
273 0 : delete [] fRecPointArray;
274 0 : fRecPointArray = tmp;
275 0 : }
276 0 : return 0;
277 : }
278 :
279 : Int_t AliHLTCaloClusterizer::CheckBuffer()
280 : {
281 : // See header file for class documentation
282 0 : if ((fAvailableSize - fUsedSize) < (Int_t)sizeof(AliHLTCaloRecPointDataStruct))
283 : {
284 0 : Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);
285 0 : Int_t digitIndexOffset = reinterpret_cast<UChar_t*>(fDigitIndexPtr) - reinterpret_cast<UChar_t*>(fRecPointDataPtr);
286 0 : UChar_t *tmp = new UChar_t[fAvailableSize*2];
287 :
288 0 : if (tmp == NULL)
289 : {
290 0 : HLTError("Pointer error");
291 0 : return(-1);
292 : }
293 :
294 0 : memcpy(tmp, fFirstRecPointPtr, fUsedSize);
295 0 : fAvailableSize *= 2;
296 0 : for (Int_t n = 0; n < fNRecPoints; n++)
297 : {
298 0 : fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));
299 : }
300 0 : delete [] fBuffer; //FR
301 0 : fBuffer = tmp; //FR
302 :
303 0 : fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);
304 0 : fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);
305 0 : fDigitIndexPtr = reinterpret_cast<Int_t*>(reinterpret_cast<UChar_t*>(fRecPointDataPtr) + digitIndexOffset);
306 : //fUsedSize = 0;
307 0 : }
308 0 : return 0;
309 0 : }
310 :
311 : void AliHLTCaloClusterizer::SetSortDigitsByPosition()
312 : {
313 : // Sort the digit pointers by position
314 0 : fCompareFunction = &CompareDigitsByPosition;
315 0 : fSortDigits = true;
316 0 : fSortedByPosition = true;
317 0 : }
318 :
319 : void AliHLTCaloClusterizer::SetSortDigitsByEnergy()
320 : {
321 : // See header file for class documentation
322 0 : fCompareFunction = &CompareDigitsByEnergy;
323 0 : fSortDigits = true;
324 0 : fSortedByEnergy = true;
325 0 : }
326 :
327 : void AliHLTCaloClusterizer::SortDigits()
328 : {
329 : // See header file for class documentation
330 0 : if (fSortDigits) qsort(fDigitsPointerArray, fNDigits, sizeof(AliHLTCaloDigitDataStruct*), fCompareFunction);
331 0 : }
332 :
333 : Int_t
334 : AliHLTCaloClusterizer::CompareDigitsByPosition(const void *dig0, const void *dig1)
335 : {
336 : // See header file for documentation
337 0 : return (*((AliHLTCaloDigitDataStruct**)(dig0)))->fID - (*((AliHLTCaloDigitDataStruct**)(dig1)))->fID;
338 : }
339 :
340 : Int_t
341 : AliHLTCaloClusterizer::CompareDigitsByEnergy(const void *dig0, const void *dig1)
342 : {
343 : // See header file for documentation
344 0 : if ( ((*((AliHLTCaloDigitDataStruct**)(dig1)))->fEnergy - (*((AliHLTCaloDigitDataStruct**)(dig0)))->fEnergy) < 0) return -1;
345 0 : return 1;
346 0 : }
347 :
348 : void AliHLTCaloClusterizer::SetDetector(TString det)
349 : {
350 0 : if(det.CompareTo("EMCAL"))
351 : {
352 0 : fIsEMCAL = true;
353 0 : }
354 : else
355 : {
356 0 : fIsEMCAL = false;
357 : }
358 0 : }
359 :
360 : Bool_t AliHLTCaloClusterizer::AreEdgeCells(AliHLTCaloDigitDataStruct *digit0, AliHLTCaloDigitDataStruct *digit1)
361 : {
362 0 : if(fIsEMCAL)
363 : {
364 0 : Int_t modDiff = digit0->fModule - digit1->fModule;
365 0 : if(TMath::Abs(modDiff) > 1) return kFALSE;
366 0 : if(digit0->fModule > digit1->fModule && digit1->fModule%2 == 0)
367 : {
368 0 : if(digit0->fZ == 0 && digit1->fZ == (fCaloConstants->GetNZROWSMOD()-1))
369 0 : return kTRUE;
370 : }
371 0 : if(digit1->fModule > digit0->fModule && digit0->fModule%2 == 0)
372 : {
373 0 : if(digit1->fZ == 0 && digit0->fZ == (fCaloConstants->GetNZROWSMOD()-1))
374 0 : return kTRUE;
375 : }
376 0 : }
377 :
378 0 : return false;
379 :
380 0 : }
|
