]> git.uio.no Git - u/mrichter/AliRoot.git/blame - HLT/CALO/AliHLTCaloClusterAnalyser.cxx
Bug fix: AliHLTTPCRawSpacePointContainer
[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"
c1e4a18c 42#include "AliHLTCaloRecoParamHandler.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),
c8fe2783 60 fClusterType(AliVCluster::kPHOSNeutral),
38a3b2ad 61 fRecoParamsPtr(0),
62 fCutOnSingleCellClusters(false),
63 fSingleCellEnergyCut(0.5)
ea367f6a 64{
65 //See header file for documentation
ea367f6a 66}
67
68AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser()
69{
98baf84d 70 // See header file for class documentation
ea367f6a 71}
72
73void
7c80a370 74AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
ea367f6a 75{
76 //see header file for documentation
77 fCaloClusterDataPtr = caloClusterDataPtr;
78}
7c80a370 79
ea367f6a 80void
7c80a370 81AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
ea367f6a 82{
7c80a370 83 fRecPointArray = recPointDataPtr;
84 fNRecPoints = nRecPoints;
85}
86
87void
f92543f5 88AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct **digits)
7c80a370 89{
90// AliHLTCaloClusterizer cl("PHOS");
91 // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
92 fDigitDataArray = digits;
93 //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
ea367f6a 94}
95
96Int_t
97AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
98{
99 //see header file for documentation
100 Float_t wtot = 0.;
101 Float_t x = 0.;
102 Float_t z = 0.;
103 Float_t xi = 0.;
104 Float_t zi = 0.;
105
106 AliHLTCaloDigitDataStruct *digit = 0;
ea367f6a 107
108 UInt_t iDigit = 0;
109
110 for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++)
111 {
7c80a370 112 AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
ad44d760 113 // digit = &(recPoint->fDigits);
98baf84d 114
7c80a370 115 Int_t *digitIndexPtr = &(recPoint->fDigits);
256f1bf0 116 wtot = 0;
117 x = 0;
118 z = 0;
119
c1e4a18c 120/* Float_t maxAmp = 0;
256f1bf0 121 Int_t maxX = 0;
c1e4a18c 122 Int_t maxZ = 0;*/
7a12c14d 123 if (fDigitDataArray[*digitIndexPtr])
124
125 for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
126 {
127
128 digit = fDigitDataArray[*digitIndexPtr];
129
130 xi = digit->fX;
131 zi = digit->fZ;
132
133 //xi = digit->fX+0.5;
134 //zi = digit->fZ+0.5;
135
136 if (recPoint->fAmp > 0 && digit->fEnergy > 0)
137 {
138 Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
139 x += xi * w ;
140 z += zi * w ;
141 wtot += w ;
142 /* if(digit->fEnergy > maxAmp)
143 {
144 maxAmp = digit->fEnergy;
145 maxX = digit->fX;// + 0.5;
146 maxZ = digit->fZ;// + 0.5;
147 }*/
148 }
149 digitIndexPtr++;
150 }
151
152 if (wtot>0)
153 {
154 recPoint->fX = x/wtot ;
155 recPoint->fZ = z/wtot ;
156 }
157 else
158 {
159 recPoint->fX = -9999;
160 recPoint->fZ =-9999;
161 // no good crashes depth with FP exception
162 //recPoint->fAmp = 0;
163 }
256f1bf0 164// printf("Max digit: E = %f, x = %d, z= %d, cluster: E = %f, x = %f, z = %f\n" , maxAmp, maxX, maxZ, recPoint->fAmp, recPoint->fX, recPoint->fZ);
ea367f6a 165 }
166 return 0;
167}
168
169
170Int_t
171AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
cd9e2797 172{
ea367f6a 173 //See header file for documentation
174 return 0;
175}
176
177Int_t
178AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
179{
180 //See header file for documentation
181 return 0;
182}
183
184
185Int_t
186AliHLTCaloClusterAnalyser::DeconvoluteClusters()
187{
188 //See header file for documentation
189 return 0;
190}
191
192Int_t
7c80a370 193AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
ea367f6a 194{
195 //See header file for documentation
196
c1e4a18c 197 if(fRecoParamsPtr)
198 {
199 fLogWeight = fRecoParamsPtr->GetLogWeight();
200 }
73d6f579 201
7c80a370 202 fNRecPoints = nRecPoints;
203
204 if(fGeometry == 0)
205 {
206 HLTError("No geometry object is initialised, creation of clusters stopped");
73cdb797 207 return -1;
7c80a370 208 }
209
210 CalculateCenterOfGravity();
211
ad44d760 212 // AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);
213 AliHLTCaloDigitDataStruct* digitPtr = 0;
214
73d6f579 215 AliHLTCaloClusterDataStruct* caloClusterPtr = 0;
7fdb2519 216
217 // Int_t id = -1;
ea367f6a 218 TVector3 globalPos;
219
220 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
221 {
73d6f579 222 if((availableSize - totSize) < sizeof(AliHLTCaloClusterDataStruct))
223 {
9a0cbaba 224 HLTError("Out of buffer: available size is: %d, total size used: %d", availableSize, totSize);
73d6f579 225 return -ENOBUFS;
226 }
38a3b2ad 227
228 AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
229
230 if(fCutOnSingleCellClusters && recPointPtr->fAmp > fSingleCellEnergyCut && recPointPtr->fMultiplicity == 1) continue;
231
7fdb2519 232 totSize += sizeof(AliHLTCaloClusterDataStruct);
73d6f579 233
234 caloClusterPtr = fCaloClusterDataPtr;
34eb3256 235 caloClusterPtr->fChi2 = 0;
236 caloClusterPtr->fClusterType = kUndef;
237 caloClusterPtr->fDispersion = 0;
238 caloClusterPtr->fDistanceToBadChannel = 0;
239 caloClusterPtr->fDistToBadChannel = 0;
240 caloClusterPtr->fEmcCpvDistance = 0;
241 caloClusterPtr->fEnergy = 0;
242 caloClusterPtr->fFitQuality = 0;
243 caloClusterPtr->fID = 0;
244 caloClusterPtr->fM02 = 0;
245 caloClusterPtr->fM20 = 0;
246 caloClusterPtr->fNCells = 0;
247 caloClusterPtr->fNExMax = 0;
248 caloClusterPtr->fTOF = 0;
fe1a0255 249 caloClusterPtr->fTrackDx = -999;
250 caloClusterPtr->fTrackDz = -999;
7c80a370 251
252 AliHLTCaloGlobalCoordinate globalCoord;
03a0ff8a 253
254 // 0 = assume photon
255 fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord, 0);
ea367f6a 256
353c33d3 257 caloClusterPtr->fModule = recPointPtr->fModule;
a6db7bcd 258 caloClusterPtr->fGlobalPos[0] = globalCoord.fX;
cd9e2797 259 caloClusterPtr->fGlobalPos[1] = globalCoord.fY;
a6db7bcd 260 caloClusterPtr->fGlobalPos[2] = globalCoord.fZ;
ea367f6a 261
5259027d 262 HLTDebug("Cluster local position: x = %f, z = %f, module = %d", recPointPtr->fX, recPointPtr->fZ, recPointPtr->fModule);
263 HLTDebug("Cluster global position: x = %f, y = %f, z = %f", globalCoord.fX, globalCoord.fY, globalCoord.fZ);
ee72b4e5 264
7fdb2519 265 //caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
266 caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
cd9e2797 267
a6db7bcd 268 caloClusterPtr->fClusterType = fClusterType;
7fdb2519 269// Int_t tmpSize = 0;//totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
a6db7bcd 270
2a24cbbe 271 //TODO remove hardcoded 10;
272 memset(caloClusterPtr->fTracksMatched, 0xff, sizeof(Int_t)*10);
273
9a0cbaba 274 //Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
cefad3b2 275 UInt_t tmpSize = (caloClusterPtr->fNCells-1)*sizeof(AliHLTCaloCellDataStruct);
cd9e2797 276
73d6f579 277 if((availableSize - totSize) < tmpSize)
278 {
9a0cbaba 279 HLTError("Out of buffer, available size is: %d, total size used: %d, extra size needed: %d", availableSize, totSize, tmpSize);
73d6f579 280 return -ENOBUFS;
281 }
c1e4a18c 282
cd9e2797 283 Int_t *digitIndexPtr = &(recPointPtr->fDigits);
284 Int_t id = 0;
285
7fdb2519 286 AliHLTCaloCellDataStruct *cellPtr = &(caloClusterPtr->fCaloCells);
0a899b5c 287 Float_t maxTime = 0; //time of maximum amplitude cell is assigned to cluster
ea367f6a 288 for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
289 {
f92543f5 290 digitPtr = fDigitDataArray[*digitIndexPtr];
936dbe15 291 fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
7fdb2519 292
293 cellPtr->fCellsAbsId= id;
294 cellPtr->fCellsAmpFraction = digitPtr->fEnergy/recPointPtr->fAmp;
0a899b5c 295 if(digitPtr->fTime > maxTime) maxTime = digitPtr->fTime;
cd9e2797 296 //printf("Cell ID pointer: %x\n", cellIDPtr);
297 //printf("Cell Amp Pointer: %x\n", cellAmpFracPtr);
9a0cbaba 298 //printf("Cell pos: x = %d, z = %d\n", digitPtr->fX, digitPtr->fZ);
34eb3256 299// printf("Cell ID: %d\n", cellPtr->fCellsAbsId);
9a0cbaba 300 //printf("Cell Amp: %f, pointer: %x\n", *cellAmpFracPtr, cellAmpFracPtr);
7fdb2519 301 cellPtr++;
cd9e2797 302 digitIndexPtr++;
303
ea367f6a 304 }
305
73d6f579 306 totSize += tmpSize;
307
c1e4a18c 308 if(fRecoParamsPtr)
309 {
310 caloClusterPtr->fEnergy = fRecoParamsPtr->GetCorrectedEnergy(recPointPtr->fAmp);
311 }
312 else
313 {
314 caloClusterPtr->fEnergy = recPointPtr->fAmp;
315 }
34eb3256 316
afcce849 317 // Set the time of the maximum digit as cluster time
0a899b5c 318 //caloClusterPtr->fTOF = recPointPtr->fTime;
319 caloClusterPtr->fTOF = maxTime;
afcce849 320
321 HLTDebug("Cluster global position: x = %f, y = %f, z = %f, energy: %f, time: %f, number of cells: %d, cluster pointer: %x", globalCoord.fX, globalCoord.fY, globalCoord.fZ, caloClusterPtr->fEnergy, caloClusterPtr->fTOF, caloClusterPtr->fNCells, caloClusterPtr);
7fdb2519 322
ea367f6a 323 if(fDoClusterFit)
324 {
325 FitCluster(recPointPtr);
326 }
327 else
328 {
329 caloClusterPtr->fDispersion = 0;
330 caloClusterPtr->fFitQuality = 0;
331 caloClusterPtr->fM20 = 0;
332 caloClusterPtr->fM02 = 0;
333
334 }
335 if(fHaveCPVInfo)
336 {
337 caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
338 }
339 else
340 {
341 caloClusterPtr->fEmcCpvDistance = -1;
342 }
343 if(fDoPID)
344 {
345 DoParticleIdentification(caloClusterPtr);
346 }
347 else
348 {
00a38d07 349 for(Int_t k = 0; k < AliPID::kSPECIESCN; k++)
ea367f6a 350 {
351 caloClusterPtr->fPID[k] = 0;
352 }
353 }
354 if(fHaveDistanceToBadChannel)
355 {
356 caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
357 }
358 else
359 {
360 caloClusterPtr->fDistanceToBadChannel = -1;
361 }
362
73d6f579 363 caloClusterPtr->fClusterType = fClusterType;
cd9e2797 364 //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);
73d6f579 365 //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
cd9e2797 366 //printf("CaloClusterPtr: %x, energy: %f\n", caloClusterPtr, caloClusterPtr->fEnergy);
367
ea367f6a 368 // caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
cd9e2797 369 //caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
7fdb2519 370 fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellPtr);
ea367f6a 371 recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
73d6f579 372 //digitPtr = &(recPointPtr->fDigits);
ea367f6a 373 }
98baf84d 374
375return fNRecPoints;
ea367f6a 376
377}
378