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