2cb92d1e457c51f35e9b3a748ef965557777fce6
[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   fRecPointArray(0),\r
44   fRecPointDataPtr(0),\r
45   fFirstRecPointPtr(0),\r
46   fArraySize(0),\r
47   fAvailableSize(0),\r
48   fUsedSize(0),\r
49   fNRecPoints(0),\r
50   fDigitIndexPtr(0),\r
51   fEmcClusteringThreshold(0),\r
52   fEmcMinEnergyThreshold(0),\r
53   fEmcTimeGate(0),\r
54   fDigitsInCluster(0),\r
55   fDigitsPointerArray(0),\r
56   fDigitContainerPtr(0),\r
57   fMaxDigitIndexDiff(0),\r
58   fNDigits(0)\r
59 {\r
60   //See header file for documentation\r
61   fEmcClusteringThreshold = 0.2;\r
62   fEmcMinEnergyThreshold = 0.03;\r
63   fEmcTimeGate = 1.e-6 ;\r
64   \r
65   fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
66   \r
67   \r
68   fArraySize = 10;\r
69   fRecPointArray = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
70   \r
71   fAvailableSize = sizeof(AliHLTCaloRecPointDataStruct) * 20;\r
72   fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(new UChar_t[fAvailableSize]);\r
73   fFirstRecPointPtr = fRecPointDataPtr;  \r
74   printf("Start of rec point data: %x, end of rec point data: %x\n", fRecPointDataPtr, reinterpret_cast<UChar_t*>(fRecPointDataPtr) + fAvailableSize);\r
75 \r
76 }//end\r
77 \r
78 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()  \r
79 {\r
80   //See header file for documentation\r
81 }\r
82 \r
83 void \r
84 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
85 {\r
86   // See header file for documentation\r
87   fRecPointDataPtr = recPointDataPtr;\r
88 }\r
89 \r
90 Int_t \r
91 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits)\r
92 {\r
93   //see header file for documentation\r
94   Int_t nRecPoints = 0;\r
95   fNRecPoints = 0;\r
96   fUsedSize = 0;\r
97   fNDigits = nDigits;\r
98   fRecPointDataPtr = fFirstRecPointPtr;\r
99   //Clusterization starts\r
100   for(Int_t i = 0; i < nDigits; i++)\r
101     { \r
102       fDigitsInCluster = 0;\r
103       //      printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
104       if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
105         {\r
106           continue;\r
107         }\r
108       CheckArray();\r
109       CheckBuffer();\r
110       //            printf("` candidate!\n");\r
111       // First digit is placed at the fDigits member variable in the recpoint\r
112       fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
113       fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
114       \r
115       fRecPointDataPtr->fAmp = 0;\r
116       fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
117 \r
118       // Assigning the digit to this rec point\r
119       fRecPointDataPtr->fDigits = i;\r
120       printf("Clusterizier: adding digit:  index pointer: %x, index: %d\n", fDigitIndexPtr, *fDigitIndexPtr);\r
121       fUsedSize += sizeof(AliHLTCaloRecPointDataStruct);\r
122       \r
123       // Incrementing the pointer to be ready for new entry\r
124       fDigitIndexPtr++;\r
125 \r
126       fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
127       fDigitsPointerArray[i]->fEnergy = 0;\r
128       fDigitsInCluster++;\r
129       nRecPoints++;\r
130 \r
131       // Scanning for the neighbours\r
132       if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
133       {\r
134          return -1;\r
135       }\r
136 \r
137       //fUsedSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);   \r
138       \r
139       fRecPointDataPtr->fMultiplicity = fDigitsInCluster;     \r
140       printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
141       printf("Multiplicity: %d\n", fDigitsInCluster);\r
142       fRecPointArray[fNRecPoints] = fRecPointDataPtr; \r
143       \r
144       fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
145 \r
146       fNRecPoints++;\r
147       \r
148     }//end of clusterization\r
149 \r
150    return nRecPoints;\r
151 }\r
152 \r
153 Int_t\r
154 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
155 {\r
156   //see header file for documentation\r
157   Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
158   Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
159 \r
160   max = fNDigits;\r
161   min = 0;\r
162   for(Int_t j = min; j < max; j++)\r
163     {\r
164       if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
165         {\r
166           if(j != index)\r
167             {\r
168               if(AreNeighbours(fDigitsPointerArray[index],\r
169                                fDigitsPointerArray[j]))\r
170                 {\r
171 //                if((fAvailableSize - fUsedSize) < sizeof(Int_t))\r
172 //                   {\r
173 //                      UChar_t *tmp = new UChar_t[fAvailableSize*2];\r
174 //                      memcpy(tmp, fRecPointDataPtr, fAvailableSize);\r
175 //                      for(Int_t n = 0; n < fNRecPoints; n++)\r
176 //                      {\r
177 //                         fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
178 //                      }\r
179 //                      fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
180 //                      fFirstRecPointPtr = fRecPointDataPtr;\r
181 //                      fUsedSize = 0;\r
182 //                   }  \r
183                   CheckBuffer();\r
184                   // Assigning index to digit\r
185                   printf("Digit index pointer: %x\n", fDigitIndexPtr);\r
186                   *fDigitIndexPtr = j;\r
187                   fUsedSize += sizeof(Int_t);\r
188                   \r
189                   printf("Clusterizier: adding digit:  index pointer: %x, index: %d\n", fDigitIndexPtr, *fDigitIndexPtr); \r
190                   // Incrementing digit pointer to be ready for new entry\r
191                   fDigitIndexPtr++;\r
192 \r
193                   recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
194                   fDigitsPointerArray[j]->fEnergy = 0;        \r
195                   fDigitsInCluster++;\r
196                   ScanForNeighbourDigits(j, recPoint);\r
197                 }\r
198             }\r
199         }\r
200     }\r
201   return 0;\r
202 }\r
203 \r
204 Int_t \r
205 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1, \r
206                                             AliHLTCaloDigitDataStruct* digit2)\r
207 {\r
208   //see header file for documentation\r
209   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
210     { \r
211       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
212       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
213       if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
214         {\r
215 //        cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
216 //          " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
217 \r
218           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
219             {\r
220               return 1; \r
221             }\r
222         }\r
223 \r
224    /*   Float_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
225       Float_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
226       if (( coldiff <= 2.4   &&  rowdiff < 0.4 ) || ( coldiff < 0.4 &&  rowdiff <= 2.4 ))\r
227         {\r
228           //      cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
229           //        " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
230 \r
231           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
232             {\r
233               return 1; \r
234             }\r
235         }\r
236    */\r
237       else\r
238         {\r
239 //          cout << "Not neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
240 //          " is not neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
241         }\r
242     }\r
243   return 0;\r
244 }\r
245 \r
246 \r
247 \r
248 Int_t AliHLTCaloClusterizer::CheckArray()\r
249 {\r
250       printf("CheckArray: fArraySize: %d, fNRecPoints: %d\n", fArraySize, fNRecPoints);\r
251       if(fArraySize == fNRecPoints)\r
252         {\r
253            printf("Expanding array...");\r
254            fArraySize *= 2;\r
255            AliHLTCaloRecPointDataStruct **tmp = new AliHLTCaloRecPointDataStruct*[fArraySize];\r
256            memcpy(tmp, fRecPointArray, fArraySize/2 * sizeof(AliHLTCaloRecPointDataStruct*));\r
257            delete fRecPointArray;\r
258            fRecPointArray = tmp;\r
259         }\r
260    return 0;\r
261 }\r
262 \r
263 Int_t AliHLTCaloClusterizer::CheckBuffer()\r
264 {\r
265    // See header file for class documentation \r
266          printf("CheckBuffer: Used size %d, fAvailableSize: %d\n", fUsedSize, fAvailableSize);\r
267         if((fAvailableSize - fUsedSize) < sizeof(AliHLTCaloRecPointDataStruct) )\r
268         {\r
269            printf("Expanding buffer...\n");\r
270             Int_t recPointOffset = reinterpret_cast<UChar_t*>(fRecPointDataPtr) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr);\r
271             fAvailableSize *= 2;\r
272             UChar_t *tmp = new UChar_t[fAvailableSize];\r
273             memcpy(tmp, fRecPointDataPtr, fAvailableSize/2);\r
274             for(Int_t n = 0; n < fNRecPoints; n++)\r
275             {\r
276                fRecPointArray[n] = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<UChar_t*>(fRecPointArray[n]) - reinterpret_cast<UChar_t*>(fFirstRecPointPtr) + reinterpret_cast<UChar_t*>(tmp));\r
277             }\r
278             delete fRecPointDataPtr;\r
279             fFirstRecPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp);\r
280             fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(tmp + recPointOffset);\r
281             fUsedSize = 0;\r
282         }\r
283    return 0;\r
284 }\r
285 \r
286 Int_t AliHLTCaloClusterizer::CheckDigits(AliHLTCaloRecPointDataStruct** recArray, AliHLTCaloDigitDataStruct** digitArray, Int_t nRP)\r
287 {\r
288   AliHLTCaloRecPointDataStruct **recpoints = recArray;\r
289   AliHLTCaloDigitDataStruct **digits = digitArray;\r
290   Int_t nRecPoints = nRP;\r
291   \r
292   if(recArray == 0)\r
293   {\r
294      recpoints = fRecPointArray;\r
295   }\r
296   if(digitArray == 0)\r
297   {\r
298      digits = fDigitsPointerArray;\r
299   }\r
300   if(nRP == 0)\r
301   {\r
302      nRecPoints = fNRecPoints;\r
303   }\r
304   printf("CL: CheckDigits: Number of rec points: %d\n", nRecPoints);\r
305   for(Int_t i = 0; i < nRecPoints; i++)\r
306   {\r
307           \r
308      AliHLTCaloRecPointDataStruct *recPoint = recpoints[i];\r
309 \r
310      //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[0];\r
311      Int_t multiplicity = recPoint->fMultiplicity;\r
312      Int_t *digitIndexPtr = &(recPoint->fDigits);\r
313      printf("CL: Rec point with energy: %f, multiplicity: %d\n", recPoint->fAmp, recPoint->fMultiplicity);\r
314      for(Int_t j = 0; j < multiplicity; j++)\r
315      {\r
316         //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[j];\r
317         AliHLTCaloDigitDataStruct *digit = digits[*digitIndexPtr];\r
318         printf("CL: Digit ID: %d, energy: %f, index: %d, indexpointer: %x\n", digit->fID, digit->fEnergy, *digitIndexPtr, digitIndexPtr);\r
319         digitIndexPtr++;\r
320         //recPoint = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitIndexPtr);\r
321      }\r
322   }\r
323      \r
324      \r
325      \r
326      \r
327      \r
328      \r
329      \r
330 }\r
331 \r
332 Int_t AliHLTCaloClusterizer::CheckDigits(AliHLTCaloRecPointDataStruct** recArray, AliHLTCaloDigitDataStruct* digitArray, Int_t nRP)\r
333 {\r
334   AliHLTCaloRecPointDataStruct **recpoints = recArray;\r
335   AliHLTCaloDigitDataStruct *digits = digitArray;\r
336   Int_t nRecPoints = nRP;\r
337   \r
338   if(recArray == 0)\r
339   {\r
340      recpoints = fRecPointArray;\r
341   }\r
342     if(nRP == 0)\r
343   {\r
344      nRecPoints = fNRecPoints;\r
345   }\r
346   printf("CL: CheckDigits: Number of rec points: %d\n", nRecPoints);\r
347   for(Int_t i = 0; i < nRecPoints; i++)\r
348   {\r
349           \r
350      AliHLTCaloRecPointDataStruct *recPoint = recpoints[i];\r
351 \r
352      //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[0];\r
353      Int_t multiplicity = recPoint->fMultiplicity;\r
354      Int_t *digitIndexPtr = &(recPoint->fDigits);\r
355      printf("CL: Rec point with energy: %f, multiplicity: %d\n", recPoint->fAmp, recPoint->fMultiplicity);\r
356      for(Int_t j = 0; j < multiplicity; j++)\r
357      {\r
358         //AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[j];\r
359         AliHLTCaloDigitDataStruct digit = digits[*digitIndexPtr];\r
360         printf("CL: digits: %x, recpoints: %x, digitIndexPtr: %x\n", digits, recpoints, digitIndexPtr);\r
361         printf("CL: Digit ID: %d, energy: %f, index: %d, indexpointer: %x\n", digit.fID, digit.fEnergy, *digitIndexPtr, digitIndexPtr);\r
362         digitIndexPtr++;\r
363         //recPoint = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitIndexPtr);\r
364      }\r
365   }\r
366      \r
367      \r
368      \r
369      \r
370      \r
371      \r
372      \r
373 }