1 // $Id: AliHLTCaloClusterAnalyser.cxx 35107 2009-09-30 01:45:06Z phille $
3 /**************************************************************************
4 * This file is property of and copyright by the ALICE HLT Project *
5 * All rights reserved. *
7 * Primary Authors: Oystein Djuvsland *
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 **************************************************************************/
19 * @file AliHLTCaloClusterAnalyser.cxx
20 * @author Oystein Djuvsland
22 * @brief Cluster analyser for Calo HLT
25 // see header file for class documentation
27 // refer to README to build package
29 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
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"
42 #include "AliHLTCaloRecoParamHandler.h"
44 ClassImp(AliHLTCaloClusterAnalyser);
46 AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() :
52 fCaloClusterDataPtr(0),
53 fCaloClusterHeaderPtr(0),
58 fHaveDistanceToBadChannel(false),
60 #ifndef HAVE_NOT_ALIVCLUSTER // backward compatibility for r42844
61 fClusterType(AliVCluster::kPHOSNeutral),
63 fClusterType(AliESDCaloCluster::kPHOSCluster),
66 fCutOnSingleCellClusters(false),
67 fSingleCellEnergyCut(0.5)
69 //See header file for documentation
72 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser()
74 // See header file for class documentation
78 AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
80 //see header file for documentation
81 fCaloClusterDataPtr = caloClusterDataPtr;
85 AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
87 fRecPointArray = recPointDataPtr;
88 fNRecPoints = nRecPoints;
92 AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct **digits)
94 // AliHLTCaloClusterizer cl("PHOS");
95 // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
96 fDigitDataArray = digits;
97 //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
101 AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
103 //see header file for documentation
110 AliHLTCaloDigitDataStruct *digit = 0;
114 for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++)
116 AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
117 // digit = &(recPoint->fDigits);
119 Int_t *digitIndexPtr = &(recPoint->fDigits);
124 /* Float_t maxAmp = 0;
127 if (fDigitDataArray[*digitIndexPtr])
129 for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
132 digit = fDigitDataArray[*digitIndexPtr];
137 //xi = digit->fX+0.5;
138 //zi = digit->fZ+0.5;
140 if (recPoint->fAmp > 0 && digit->fEnergy > 0)
142 Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
146 /* if(digit->fEnergy > maxAmp)
148 maxAmp = digit->fEnergy;
149 maxX = digit->fX;// + 0.5;
150 maxZ = digit->fZ;// + 0.5;
158 recPoint->fX = x/wtot ;
159 recPoint->fZ = z/wtot ;
163 recPoint->fX = -9999;
165 // no good crashes depth with FP exception
166 //recPoint->fAmp = 0;
168 // 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);
175 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
177 //See header file for documentation
182 AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
184 //See header file for documentation
190 AliHLTCaloClusterAnalyser::DeconvoluteClusters()
192 //See header file for documentation
197 AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
199 //See header file for documentation
203 fLogWeight = fRecoParamsPtr->GetLogWeight();
206 fNRecPoints = nRecPoints;
210 HLTError("No geometry object is initialised, creation of clusters stopped");
214 CalculateCenterOfGravity();
216 // AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);
217 AliHLTCaloDigitDataStruct* digitPtr = 0;
219 AliHLTCaloClusterDataStruct* caloClusterPtr = 0;
224 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
226 if((availableSize - totSize) < sizeof(AliHLTCaloClusterDataStruct))
228 HLTError("Out of buffer: available size is: %d, total size used: %d", availableSize, totSize);
232 AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
234 if(fCutOnSingleCellClusters && recPointPtr->fAmp > fSingleCellEnergyCut && recPointPtr->fMultiplicity == 1) continue;
236 totSize += sizeof(AliHLTCaloClusterDataStruct);
238 caloClusterPtr = fCaloClusterDataPtr;
239 caloClusterPtr->fChi2 = 0;
240 caloClusterPtr->fClusterType = kUndef;
241 caloClusterPtr->fDispersion = 0;
242 caloClusterPtr->fDistanceToBadChannel = 0;
243 caloClusterPtr->fDistToBadChannel = 0;
244 caloClusterPtr->fEmcCpvDistance = 0;
245 caloClusterPtr->fEnergy = 0;
246 caloClusterPtr->fFitQuality = 0;
247 caloClusterPtr->fID = 0;
248 caloClusterPtr->fM02 = 0;
249 caloClusterPtr->fM20 = 0;
250 caloClusterPtr->fNCells = 0;
251 caloClusterPtr->fNExMax = 0;
252 caloClusterPtr->fTOF = 0;
253 caloClusterPtr->fTrackDx = -999;
254 caloClusterPtr->fTrackDz = -999;
256 AliHLTCaloGlobalCoordinate globalCoord;
259 fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord, 0);
261 caloClusterPtr->fModule = recPointPtr->fModule;
262 caloClusterPtr->fGlobalPos[0] = globalCoord.fX;
263 caloClusterPtr->fGlobalPos[1] = globalCoord.fY;
264 caloClusterPtr->fGlobalPos[2] = globalCoord.fZ;
266 HLTDebug("Cluster local position: x = %f, z = %f, module = %d", recPointPtr->fX, recPointPtr->fZ, recPointPtr->fModule);
267 HLTDebug("Cluster global position: x = %f, y = %f, z = %f", globalCoord.fX, globalCoord.fY, globalCoord.fZ);
269 //caloClusterPtr->fNCells = 0;//recPointPtr->fMultiplicity;
270 caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
272 caloClusterPtr->fClusterType = fClusterType;
273 // Int_t tmpSize = 0;//totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
275 //TODO remove hardcoded 10;
276 memset(caloClusterPtr->fTracksMatched, 0xff, sizeof(Int_t)*10);
278 //Int_t tmpSize = totSize + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
279 UInt_t tmpSize = (caloClusterPtr->fNCells-1)*sizeof(AliHLTCaloCellDataStruct);
281 if((availableSize - totSize) < tmpSize)
283 HLTError("Out of buffer, available size is: %d, total size used: %d, extra size needed: %d", availableSize, totSize, tmpSize);
287 Int_t *digitIndexPtr = &(recPointPtr->fDigits);
290 AliHLTCaloCellDataStruct *cellPtr = &(caloClusterPtr->fCaloCells);
291 Float_t maxTime = 0; //time of maximum amplitude cell is assigned to cluster
292 for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
294 digitPtr = fDigitDataArray[*digitIndexPtr];
295 fGeometry->GetCellAbsId(recPointPtr->fModule, digitPtr->fX, digitPtr->fZ, id);
297 cellPtr->fCellsAbsId= id;
298 cellPtr->fCellsAmpFraction = digitPtr->fEnergy/recPointPtr->fAmp;
299 if(digitPtr->fTime > maxTime) maxTime = digitPtr->fTime;
300 //printf("Cell ID pointer: %x\n", cellIDPtr);
301 //printf("Cell Amp Pointer: %x\n", cellAmpFracPtr);
302 //printf("Cell pos: x = %d, z = %d\n", digitPtr->fX, digitPtr->fZ);
303 // printf("Cell ID: %d\n", cellPtr->fCellsAbsId);
304 //printf("Cell Amp: %f, pointer: %x\n", *cellAmpFracPtr, cellAmpFracPtr);
314 caloClusterPtr->fEnergy = fRecoParamsPtr->GetCorrectedEnergy(recPointPtr->fAmp);
318 caloClusterPtr->fEnergy = recPointPtr->fAmp;
321 // Set the time of the maximum digit as cluster time
322 //caloClusterPtr->fTOF = recPointPtr->fTime;
323 caloClusterPtr->fTOF = maxTime;
325 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);
329 FitCluster(recPointPtr);
333 caloClusterPtr->fDispersion = 0;
334 caloClusterPtr->fFitQuality = 0;
335 caloClusterPtr->fM20 = 0;
336 caloClusterPtr->fM02 = 0;
341 caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
345 caloClusterPtr->fEmcCpvDistance = -1;
349 DoParticleIdentification(caloClusterPtr);
353 for(Int_t k = 0; k < AliPID::kSPECIESCN; k++)
355 caloClusterPtr->fPID[k] = 0;
358 if(fHaveDistanceToBadChannel)
360 caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
364 caloClusterPtr->fDistanceToBadChannel = -1;
367 caloClusterPtr->fClusterType = fClusterType;
368 //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);
369 //totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));
370 //printf("CaloClusterPtr: %x, energy: %f\n", caloClusterPtr, caloClusterPtr->fEnergy);
372 // caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
373 //caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
374 fCaloClusterDataPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellPtr);
375 recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
376 //digitPtr = &(recPointPtr->fDigits);