17c2b3b99b47edc0e0a3b5cd0870b2e24139a59c
[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   fRecPointDataPtr(0),\r
44   fDigitIndexPtr(0),\r
45   fEmcClusteringThreshold(0),\r
46   fEmcMinEnergyThreshold(0),\r
47   fEmcTimeGate(0),\r
48   fDigitsInCluster(0),\r
49   fDigitsPointerArray(0),\r
50   fDigitContainerPtr(0),\r
51   fMaxDigitIndexDiff(0),\r
52   fNDigits(0)\r
53 {\r
54   //See header file for documentation\r
55   fEmcClusteringThreshold = 0.2;\r
56   fEmcMinEnergyThreshold = 0.03;\r
57   fEmcTimeGate = 1.e-6 ;\r
58   \r
59   fMaxDigitIndexDiff = 2*fCaloConstants->GetNZROWSMOD();\r
60 \r
61 }//end\r
62 \r
63 \r
64 //BALLE how do you set the right detector?\r
65 // AliHLTCaloClusterizer::AliHLTCaloClusterizer(const AliHLTCaloClusterizer &) :\r
66 //   AliHLTCaloConstantsHandler("BALLE"),\r
67 //   fRecPointDataPtr(0),\r
68 //   fDigitDataPtr(0),\r
69 //   fEmcClusteringThreshold(0),\r
70 //   fEmcMinEnergyThreshold(0),\r
71 //   fEmcTimeGate(0),\r
72 //   fDigitsInCluster(0),\r
73 //   fDigitContainerPtr(0),\r
74 //   fMaxDigitIndexDiff(0)\r
75 // {\r
76 //   // dummy copy constructor\r
77 // }//end\r
78 \r
79 \r
80 AliHLTCaloClusterizer::~AliHLTCaloClusterizer()  \r
81 {\r
82   //See header file for documentation\r
83 }\r
84 \r
85 void \r
86 AliHLTCaloClusterizer::SetRecPointDataPtr(AliHLTCaloRecPointDataStruct* recPointDataPtr)\r
87 {\r
88   // See header file for documentation\r
89   fRecPointDataPtr = recPointDataPtr;\r
90 }\r
91 \r
92 Int_t \r
93 AliHLTCaloClusterizer::ClusterizeEvent(Int_t nDigits, UInt_t availableSize, UInt_t& totSize)\r
94 {\r
95   //see header file for documentation\r
96   Int_t nRecPoints = 0;\r
97   \r
98   fAvailableSize = availableSize;\r
99   \r
100   fNDigits = nDigits;\r
101 \r
102   UInt_t maxRecPointSize = sizeof(AliHLTCaloRecPointDataStruct) + (sizeof(AliHLTCaloDigitDataStruct) << 7); //Reasonable estimate... \r
103 \r
104   //Clusterization starts\r
105   for(Int_t i = 0; i < nDigits; i++)\r
106     { \r
107       fDigitsInCluster = 0;\r
108       //      printf("ENERGY: %f\n", fDigitsPointerArray[i]->fEnergy);\r
109       if(fDigitsPointerArray[i]->fEnergy < fEmcClusteringThreshold)\r
110         {\r
111           continue;\r
112         }\r
113         \r
114         if(fAvailableSize < (sizeof(AliHLTCaloRecPointDataStruct)))\r
115         {\r
116           HLTError("Out of buffer, stopping clusterisation");\r
117           return -1; \r
118         }\r
119         \r
120       //            printf("cluster candidate!\n");\r
121       // First digit is placed at the fDigits member variable in the recpoint\r
122       fDigitIndexPtr = &(fRecPointDataPtr->fDigits);\r
123 \r
124       fRecPointDataPtr->fAmp = 0;\r
125       fRecPointDataPtr->fModule = fDigitsPointerArray[i]->fModule;\r
126 \r
127       // Assigning the digit to this rec point\r
128       fRecPointDataPtr->fDigits = i;\r
129 \r
130       // Incrementing the pointer to be ready for new entry\r
131       fDigitIndexPtr++;\r
132 \r
133       fRecPointDataPtr->fAmp += fDigitsPointerArray[i]->fEnergy;\r
134       fDigitsPointerArray[i]->fEnergy = 0;\r
135       fDigitsInCluster++;\r
136       nRecPoints++;\r
137 \r
138       // Scanning for the neighbours\r
139       if(ScanForNeighbourDigits(i, fRecPointDataPtr) != 0)\r
140       {\r
141          return -1;\r
142       }\r
143 \r
144       totSize += sizeof(AliHLTCaloRecPointDataStruct) + (fDigitsInCluster-1)*sizeof(AliHLTCaloDigitDataStruct);   \r
145       fRecPointDataPtr->fMultiplicity = fDigitsInCluster;     \r
146       //      printf("Rec point energy: %f\n", fRecPointDataPtr->fAmp);\r
147       fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(fDigitIndexPtr);\r
148 \r
149     }//end of clusterization\r
150 \r
151    return nRecPoints;\r
152 }\r
153 \r
154 Int_t\r
155 AliHLTCaloClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTCaloRecPointDataStruct* recPoint)\r
156 {\r
157   //see header file for documentation\r
158   Int_t max = TMath::Min(fNDigits, (Int_t)fMaxDigitIndexDiff+index);\r
159   Int_t min = TMath::Max(0, (Int_t)(index - (Int_t)fMaxDigitIndexDiff));\r
160 \r
161   max = fNDigits;\r
162   min = 0;\r
163   for(Int_t j = min; j < max; j++)\r
164     {\r
165       if(fDigitsPointerArray[j]->fEnergy > fEmcMinEnergyThreshold)\r
166         {\r
167           if(j != index)\r
168             {\r
169               if(AreNeighbours(fDigitsPointerArray[index],\r
170                                fDigitsPointerArray[j]))\r
171                 {\r
172                   if(fAvailableSize < (sizeof(Int_t)))\r
173                      {\r
174                         HLTError("Out of buffer, stopping clusterisation");\r
175                         return -1; \r
176                      }  \r
177                   // Assigning index to digit\r
178                   *fDigitIndexPtr = j;\r
179                   // Incrementing digit pointer to be ready for new entry\r
180                   fDigitIndexPtr++;\r
181 \r
182                   recPoint->fAmp += fDigitsPointerArray[j]->fEnergy;\r
183                   fDigitsPointerArray[j]->fEnergy = 0;        \r
184                   fDigitsInCluster++;\r
185                   ScanForNeighbourDigits(j, recPoint);\r
186                 }\r
187             }\r
188         }\r
189     }\r
190   return 0;\r
191 }\r
192 \r
193 Int_t \r
194 AliHLTCaloClusterizer::AreNeighbours(AliHLTCaloDigitDataStruct* digit1, \r
195                                             AliHLTCaloDigitDataStruct* digit2)\r
196 {\r
197   //see header file for documentation\r
198   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module\r
199     { \r
200 //       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
201 //       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
202 //       if (( coldiff <= 1   &&  rowdiff == 0 ) || ( coldiff == 0 &&  rowdiff <= 1 ))\r
203 //      {\r
204 //        cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
205 //          " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
206 \r
207 //        if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
208 //          {\r
209 //            return 1; \r
210 //          }\r
211 //      }\r
212 \r
213       Float_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  \r
214       Float_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); \r
215       if (( coldiff <= 2.4   &&  rowdiff < 0.4 ) || ( coldiff < 0.4 &&  rowdiff <= 2.4 ))\r
216         {\r
217           //      cout << "Are neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
218           //        " is neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
219 \r
220           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)\r
221             {\r
222               return 1; \r
223             }\r
224         }\r
225       else\r
226         {\r
227           //  cout << "Not neighbours: digit (E = "  << digit1->fEnergy << ") with x = " << digit1->fX << " and z = " << digit1->fZ << \r
228           //  " is not neighbour with digit (E = " << digit2->fEnergy << ") with x = " << digit2->fX << " and z = " << digit2->fZ << endl;\r
229         }\r
230     }\r
231   return 0;\r
232 }\r