Fast online clusterizer now working (Oystein)
[u/mrichter/AliRoot.git] / HLT / PHOS / AliHLTPHOSClusterizer.cxx
CommitLineData
9cc0deb1 1//Insert copyright
aac22523 2
3
aac22523 4#include "AliHLTPHOSClusterizer.h"
9cc0deb1 5#include "AliHLTPHOSBase.h"
aac22523 6#include "TMath.h"
9cc0deb1 7#include "AliHLTPHOSRecPointContainerStruct.h"
91b95d47 8#include "AliHLTPHOSRecPointDataStruct.h"
9cc0deb1 9#include "AliHLTPHOSDigitDataStruct.h"
10#include "AliHLTPHOSDigitContainerDataStruct.h"
11#include "TClonesArray.h"
12#include "AliPHOSGeometry.h"
13#include "AliPHOSDigit.h"
14#include "AliPHOSRecoParamEmc.h"
15
16ClassImp(AliHLTPHOSClusterizer);
17
18AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():
19 AliHLTPHOSBase(),
20 fEmcClusteringThreshold(0),
21 fEmcMinEnergyThreshold(0),
22 fEmcTimeGate(0),
23 fLogWeight(0),
24 fDigitsInCluster(0),
25 fOnlineMode(true),
26 fDigitArrayPtr(0),
27 fEmcRecPointsPtr(0),
28 fDigitPtr(0),
29 fDigitContainerPtr(0),
30 fRecPointContainerPtr(0),
31 fPHOSGeometry(0),
32 fGetterPtr(0)
aac22523 33{
6e709a0d 34 //See header file for documentation
9cc0deb1 35 fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
36 fEmcClusteringThreshold = 0.2;
37 fEmcMinEnergyThreshold = 0.03;
38 fEmcTimeGate = 1.e-8 ;
39 fLogWeight = 4.5;
aac22523 40
9cc0deb1 41
aac22523 42}//end
43
9c9d15d6 44
9cc0deb1 45AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &)
46{
91b95d47 47 //Copy constructor, not implemented
9cc0deb1 48}//end
aac22523 49
9cc0deb1 50AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()
aac22523 51{
91b95d47 52 //Destructor
aac22523 53}
54
9cc0deb1 55void
56AliHLTPHOSClusterizer::SetRecoParameters(AliPHOSRecoParamEmc* params)
57{
58 fEmcClusteringThreshold = params->GetClusteringThreshold();
59 fEmcMinEnergyThreshold = params->GetMinE();
60 fLogWeight = params->GetLogWeight();
61}
aac22523 62
9cc0deb1 63void
64AliHLTPHOSClusterizer::SetOfflineMode(AliPHOSGetter* getter)
aac22523 65{
9cc0deb1 66 fRecPointContainerPtr = new AliHLTPHOSRecPointContainerStruct();
67 fDigitContainerPtr = new AliHLTPHOSDigitContainerDataStruct();
68 fGetterPtr = getter;
69 fDigitArrayPtr = fGetterPtr->Digits();
70 fEmcRecPointsPtr = fGetterPtr->EmcRecPoints();
71 fOnlineMode = false;
72}
aac22523 73
9cc0deb1 74Int_t
75AliHLTPHOSClusterizer::GetEvent(Int_t i)
76{
77 Int_t coord[4];
aac22523 78
9cc0deb1 79 fGetterPtr->Event(i, "D");
80 for(Int_t j = 0; j < fDigitArrayPtr->GetEntries(); j++)
aac22523 81 {
9cc0deb1 82 fDigitPtr = (AliPHOSDigit*)fDigitArrayPtr->At(j);
83 fPHOSGeometry->AbsToRelNumbering(fDigitPtr->GetId(), coord);
84 fDigitContainerPtr->fDigitDataStruct[j].fX = coord[3];
85 fDigitContainerPtr->fDigitDataStruct[j].fZ = coord[2];
86 fDigitContainerPtr->fDigitDataStruct[j].fModule = coord[0];
87 fDigitContainerPtr->fDigitDataStruct[j].fEnergy = fDigitPtr->GetEnergy();
88 fDigitContainerPtr->fDigitDataStruct[j].fTime = fDigitPtr->GetTime();
89 }
90 fDigitContainerPtr->fNDigits = fDigitArrayPtr->GetEntriesFast();
91 return 0;
92}
aac22523 93
9cc0deb1 94Int_t
95AliHLTPHOSClusterizer::GetNEvents()
96{
97 if(fOnlineMode)
98 {
99 printf("Number of events not available in online mode!\n");
100 return -1;
101 }
102 return fGetterPtr->MaxEvent();
103}
aac22523 104
aac22523 105
9cc0deb1 106Int_t
107AliHLTPHOSClusterizer::ClusterizeEvent()
aac22523 108{
9cc0deb1 109 // Steering method to construct the clusters stored in a list of Reconstructed Points
110 // A cluster is defined as a list of neighbour digits
aac22523 111 Int_t nRecPoints = 0;
9cc0deb1 112 UInt_t i = 0;
aac22523 113
6e709a0d 114
9cc0deb1 115 AliHLTPHOSRecPointDataStruct *recPoint = 0;
116
117 //Clusterization starts
118 for(i = 0; i < fDigitContainerPtr->fNDigits; i++)
119 {
120
121 fDigitsInCluster = 0;
aac22523 122
9cc0deb1 123 if(fDigitContainerPtr->fDigitDataStruct[i].fEnergy < fEmcClusteringThreshold)
aac22523 124 {
9cc0deb1 125 continue;
aac22523 126 }
aac22523 127
9cc0deb1 128 recPoint = &(fRecPointContainerPtr->fRecPointArray[nRecPoints]);
129 recPoint->fAmp = 0;
130 //recPoint->
131 recPoint->fDigitsList[fDigitsInCluster] = fDigitContainerPtr->fDigitDataStruct[i];
132 recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[i].fEnergy;
133 fDigitContainerPtr->fDigitDataStruct[i].fEnergy = 0;
134 fDigitsInCluster++;
135 nRecPoints++;
136 ScanForNeighbourDigits(i, recPoint);
137
138 recPoint->fMultiplicity = fDigitsInCluster;
139 }//end of clusterization
140 fRecPointContainerPtr->fNRecPoints = nRecPoints;
aac22523 141
9cc0deb1 142 return nRecPoints;
143}
aac22523 144
9cc0deb1 145void
146AliHLTPHOSClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTPHOSRecPointDataStruct* recPoint)
aac22523 147
aac22523 148{
9cc0deb1 149
150 for(UInt_t j = 0; j < fDigitContainerPtr->fNDigits; j++)
aac22523 151 {
9cc0deb1 152 if(fDigitContainerPtr->fDigitDataStruct[j].fEnergy > fEmcMinEnergyThreshold)
153 {
154 switch(AreNeighbours(&(fDigitContainerPtr->fDigitDataStruct[index]),
155 &(fDigitContainerPtr->fDigitDataStruct[j])))
156 {
157 case 0:
158 break;
159 case 1:
160 recPoint->fDigitsList[fDigitsInCluster] = fDigitContainerPtr->fDigitDataStruct[j];
161 recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[j].fEnergy;
162 fDigitContainerPtr->fDigitDataStruct[j].fEnergy = 0;
163 fDigitsInCluster++;
164 ScanForNeighbourDigits(j, recPoint);
165 break;
166 case 2:
167 break;
168 }
aac22523 169 }
9cc0deb1 170 }
171 return;
172}
6e709a0d 173
aac22523 174
9cc0deb1 175Int_t
176AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1,
177 AliHLTPHOSDigitDataStruct* digit2)
6e709a0d 178{
6e709a0d 179
9cc0deb1 180 //Int_t coord1[4];
181 //Int_t coord2[4];
6e709a0d 182
9cc0deb1 183 // fPHOSGeometry->AbsToRelNumbering(digit1->fID, coord1);
184 //fPHOSGeometry->AbsToRelNumbering(digit2->fID, coord2);
6e709a0d 185
d2b84453 186
9cc0deb1 187 if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
188 {
d2b84453 189
9cc0deb1 190 Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );
191 Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX );
192
193 if (( coldiff <= 1 ) && ( rowdiff <= 1 ))
6e709a0d 194 {
9cc0deb1 195 if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
196 return 1;
6e709a0d 197 }
6e709a0d 198 }
199 return 0;
200}
aac22523 201
9cc0deb1 202void
203AliHLTPHOSClusterizer::CalculateCenterOfGravity()
aac22523 204{
aac22523 205
9cc0deb1 206 // Calculates the center of gravity in the local PHOS-module coordinates
207 Float_t wtot = 0.;
208
209 Int_t relid[4];
aac22523 210
9cc0deb1 211 Float_t x = 0.;
212 Float_t z = 0.;
213 Float_t xi = 0.;
214 Float_t zi = 0.;
aac22523 215
9cc0deb1 216 AliHLTPHOSRecPointDataStruct *recPoint = 0;
217 AliHLTPHOSDigitDataStruct *digit = 0;
aac22523 218
9cc0deb1 219 AliPHOSGeometry * phosgeom = AliPHOSGeometry::GetInstance() ;
aac22523 220
9cc0deb1 221 Int_t iDigit = 0;
222 Int_t iRecPoint = 0;
91b95d47 223
9cc0deb1 224 for(iRecPoint=0; iRecPoint<fRecPointContainerPtr->fNRecPoints; iRecPoint++)
aac22523 225 {
9cc0deb1 226 recPoint = &(fRecPointContainerPtr->fRecPointArray[iRecPoint]);
227 for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
aac22523 228 {
9cc0deb1 229 digit = &(recPoint->fDigitsList[iDigit]);
230
231 //fPHOSGeometry->AbsToRelNumbering(digit->fID, relid) ;
232 // fPHOSGeometry->RelPosInModule(relid, xi, zi);
233 xi = digit->fX;
234 zi = digit->fZ;
235
236 if (recPoint->fAmp > 0 && digit->fEnergy > 0)
237 {
238 Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
239 x += xi * w ;
240 z += zi * w ;
241 wtot += w ;
242 }
243 }
244
245 if (wtot>0)
246 {
247 recPoint->fX = x/wtot ;
248 recPoint->fZ = z/wtot ;
249 }
250 else
251 {
252 recPoint->fAmp = 0;
aac22523 253 }
254 }
9cc0deb1 255
256}
aac22523 257
aac22523 258