- fixes in the handling of the AliHLTCaloClusterStructs
[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterAnalyser.cxx
CommitLineData
ea367f6a 1// $Id: AliHLTCaloClusterAnalyser.cxx 35107 2009-09-30 01:45:06Z phille $
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 AliHLTCaloClusterAnalyser.cxx
20 * @author Oystein Djuvsland
21 * @date
22 * @brief Cluster analyser for Calo 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 "AliHLTCaloClusterAnalyser.h"
32#include "AliHLTCaloRecPointHeaderStruct.h"
33#include "AliHLTCaloRecPointDataStruct.h"
34#include "AliHLTCaloClusterDataStruct.h"
35//#include "AliHLTCaloPhysicsAnalyzer.h"
7c80a370 36#include "AliHLTCaloGeometry.h"
ea367f6a 37#include "AliESDCaloCluster.h"
38#include "TMath.h"
39#include "TVector3.h"
73d6f579 40#include "TH1F.h"
41#include "TFile.h"
7c80a370 42#include "AliHLTCaloClusterizer.h"
ea367f6a 43
44ClassImp(AliHLTCaloClusterAnalyser);
45
7c80a370 46AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() :
ea367f6a 47 // AliHLTCaloBase(),
7c80a370 48 fLogWeight(4.5),
49 fRecPointArray(0),
50 fDigitDataArray(0),
ea367f6a 51 fNRecPoints(0),
52 fCaloClusterDataPtr(0),
53 fCaloClusterHeaderPtr(0),
ea367f6a 54 //fAnalyzerPtr(0),
55 fDoClusterFit(false),
56 fHaveCPVInfo(false),
57 fDoPID(false),
7c80a370 58 fHaveDistanceToBadChannel(false),
73d6f579 59 fGeometry(0),
60 fClusterType(AliESDCaloCluster::kPHOSCluster)
ea367f6a 61{
62 //See header file for documentation
73d6f579 63 fHist = new TH1F("cluster_energies", "cluster_energies", 200, 0, 20);
64
ea367f6a 65}
66
67AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser()
68{
73d6f579 69 TFile *file = TFile::Open("debug.root", "RECREATE");
70 fHist->Write();
71 file->Close();
ea367f6a 72}
73
74void
7c80a370 75AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
ea367f6a 76{
77 //see header file for documentation
78 fCaloClusterDataPtr = caloClusterDataPtr;
79}
7c80a370 80
ea367f6a 81void
7c80a370 82AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
ea367f6a 83{
7c80a370 84 fRecPointArray = recPointDataPtr;
85 fNRecPoints = nRecPoints;
86}
87
88void
89AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct *digits)
90{
91// AliHLTCaloClusterizer cl("PHOS");
92 // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
93 fDigitDataArray = digits;
94 //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
ea367f6a 95}
96
97Int_t
98AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
99{
100 //see header file for documentation
101 Float_t wtot = 0.;
102 Float_t x = 0.;
103 Float_t z = 0.;
104 Float_t xi = 0.;
105 Float_t zi = 0.;
106
107 AliHLTCaloDigitDataStruct *digit = 0;
108 //AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
7c80a370 109
110
111 //AliHLTCaloClusterizer cl("PHOS");
112 //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
ea367f6a 113
114 UInt_t iDigit = 0;
115
116 for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++)
117 {
7c80a370 118 AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
ad44d760 119 // digit = &(recPoint->fDigits);
7c80a370 120 cout << "COG: Rec point multiplicity: " << recPoint->fMultiplicity << ", rec point energy: " << recPoint->fAmp << endl;
121 Int_t *digitIndexPtr = &(recPoint->fDigits);
122
ea367f6a 123 for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
124 {
7c80a370 125 //cout << "1. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;
126
127 digit = &(fDigitDataArray[*digitIndexPtr]);
128 cout << "COG: 2. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", index pointer: " << digitIndexPtr<< endl;
129 cout << ", dig Energy: " << digit->fEnergy << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;
ea367f6a 130 xi = digit->fX;
131 zi = digit->fZ;
132 // cout << "COG digits (x:z:E:time): " << xi << " : " << zi << " : " << digit->fEnergy << " : " << digit->fTime << endl;
133 if (recPoint->fAmp > 0 && digit->fEnergy > 0)
134 {
135 Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
136 x += xi * w ;
137 z += zi * w ;
138 wtot += w ;
139 }
7c80a370 140 digitIndexPtr++;
ea367f6a 141 }
142 //cout << endl;
143 if (wtot>0)
144 {
145 recPoint->fX = x/wtot ;
146 recPoint->fZ = z/wtot ;
147 }
148 else
149 {
150 recPoint->fAmp = 0;
151 }
ea367f6a 152 }
153 return 0;
154}
155
156
157Int_t
158AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
159{
160 //See header file for documentation
161 return 0;
162}
163
164Int_t
165AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
166{
167 //See header file for documentation
168 return 0;
169}
170
171
172Int_t
173AliHLTCaloClusterAnalyser::DeconvoluteClusters()
174{
175 //See header file for documentation
176 return 0;
177}
178
179Int_t
7c80a370 180AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
ea367f6a 181{
182 //See header file for documentation
183
73d6f579 184
b4479a87 185 totSize += sizeof(AliHLTCaloClusterDataStruct);
7c80a370 186 fNRecPoints = nRecPoints;
187
188 if(fGeometry == 0)
189 {
190 HLTError("No geometry object is initialised, creation of clusters stopped");
191 }
192
193 CalculateCenterOfGravity();
194
ad44d760 195 // AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);
196 AliHLTCaloDigitDataStruct* digitPtr = 0;
197
73d6f579 198 AliHLTCaloClusterDataStruct* caloClusterPtr = 0;
199 UShort_t* cellIDPtr = 0;
200 Float_t* cellAmpFracPtr = 0;;
ea367f6a 201
202 Int_t id = -1;
203 TVector3 globalPos;
204
205 for(Int_t i = 0; i < fNRecPoints; i++) //TODO needs fix when we start unfolding (number of clusters not necessarily same as number of recpoints gotten from the clusterizer
206 {
73d6f579 207 if((availableSize - totSize) < sizeof(AliHLTCaloClusterDataStruct))
208 {
209 HLTError("Out of buffer");
210 return -ENOBUFS;
211 }
212
213 caloClusterPtr = fCaloClusterDataPtr;
214
215 cellIDPtr = &(caloClusterPtr->fCellsAbsId);
216 cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
217
7c80a370 218 AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
219 cout << "CA: rec point energy: " << recPointPtr->fAmp << ", rec point multiplicity: " << recPointPtr->fMultiplicity << endl;
ea367f6a 220 // cout << "Local Position (x:z:module): " << recPointPtr->fX << " : "<< recPointPtr->fZ << " : " << recPointPtr->fModule << endl;
7c80a370 221
222 AliHLTCaloGlobalCoordinate globalCoord;
223 fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord);
ea367f6a 224 // cout << "Global Position (x:y:z): " << globalPos[0] << " : "<< globalPos[1] << " : " << globalPos[2] << endl << endl;
225
7c80a370 226 caloClusterPtr->fGlobalPos[0] = globalCoord.fX;
227 caloClusterPtr->fGlobalPos[1] = globalCoord.fY;
228 caloClusterPtr->fGlobalPos[2] = globalCoord.fZ;
ea367f6a 229
73d6f579 230 caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
231
232 Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
233
234 if((availableSize - totSize) < tmpSize)
235 {
236 HLTError("Out of buffer");
237 return -ENOBUFS;
238 }
239
ea367f6a 240 for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
241 {
73d6f579 242/* fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
243 *cellIDPtr = id;
ea367f6a 244 *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
245 digitPtr++;
246 cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t));
7c80a370 247 cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Short_t)); */
ea367f6a 248 }
249
73d6f579 250 totSize += tmpSize;
251
ea367f6a 252 caloClusterPtr->fEnergy = recPointPtr->fAmp;
73d6f579 253
254 fHist->Fill(caloClusterPtr->fEnergy);
255
7c80a370 256 cout << "CA: cluster energy: " << caloClusterPtr->fEnergy << endl;
257 cout << "CA: recpoint energy: " << recPointPtr->fAmp << endl;
73d6f579 258
259
ea367f6a 260 if(fDoClusterFit)
261 {
262 FitCluster(recPointPtr);
263 }
264 else
265 {
266 caloClusterPtr->fDispersion = 0;
267 caloClusterPtr->fFitQuality = 0;
268 caloClusterPtr->fM20 = 0;
269 caloClusterPtr->fM02 = 0;
270
271 }
272 if(fHaveCPVInfo)
273 {
274 caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
275 }
276 else
277 {
278 caloClusterPtr->fEmcCpvDistance = -1;
279 }
280 if(fDoPID)
281 {
282 DoParticleIdentification(caloClusterPtr);
283 }
284 else
285 {
286 for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
287 {
288 caloClusterPtr->fPID[k] = 0;
289 }
290 }
291 if(fHaveDistanceToBadChannel)
292 {
293 caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
294 }
295 else
296 {
297 caloClusterPtr->fDistanceToBadChannel = -1;
298 }
299
73d6f579 300 caloClusterPtr->fClusterType = fClusterType;
ea367f6a 301 // totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);
73d6f579 302 //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
ea367f6a 303
304 // caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
305 caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
7c80a370 306
ea367f6a 307 recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
73d6f579 308 //digitPtr = &(recPointPtr->fDigits);
ea367f6a 309 }
310 // cout << "CA: Energy End: " << fCaloClusterDataPtr->fEnergy << endl;
311 //cout << "CA totSize: " << totSize << endl;
312 return fNRecPoints;
313
314}
315