Fast online clusterizer now working (Oystein)
[u/mrichter/AliRoot.git] / HLT / PHOS / AliHLTPHOSClusterizer.cxx
1 //Insert copyright
2
3
4 #include "AliHLTPHOSClusterizer.h"
5 #include "AliHLTPHOSBase.h"
6 #include "TMath.h"
7 #include "AliHLTPHOSRecPointContainerStruct.h"
8 #include "AliHLTPHOSRecPointDataStruct.h"
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
16 ClassImp(AliHLTPHOSClusterizer);
17
18 AliHLTPHOSClusterizer::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)
33 {
34   //See header file for documentation
35   fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
36   fEmcClusteringThreshold = 0.2;
37   fEmcMinEnergyThreshold = 0.03;
38   fEmcTimeGate = 1.e-8 ;
39   fLogWeight = 4.5;
40   
41
42 }//end
43
44
45 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &)
46 {
47   //Copy constructor, not implemented
48 }//end
49
50 AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()  
51 {
52   //Destructor
53 }
54
55 void
56 AliHLTPHOSClusterizer::SetRecoParameters(AliPHOSRecoParamEmc* params)
57 {
58   fEmcClusteringThreshold = params->GetClusteringThreshold();
59   fEmcMinEnergyThreshold = params->GetMinE();
60   fLogWeight = params->GetLogWeight();
61 }
62
63 void
64 AliHLTPHOSClusterizer::SetOfflineMode(AliPHOSGetter* getter)
65 {
66   fRecPointContainerPtr = new AliHLTPHOSRecPointContainerStruct();
67   fDigitContainerPtr = new AliHLTPHOSDigitContainerDataStruct();
68   fGetterPtr = getter;
69   fDigitArrayPtr = fGetterPtr->Digits();
70   fEmcRecPointsPtr = fGetterPtr->EmcRecPoints();
71   fOnlineMode = false;
72 }
73
74 Int_t
75 AliHLTPHOSClusterizer::GetEvent(Int_t i)
76 {
77   Int_t coord[4];
78
79   fGetterPtr->Event(i, "D");
80   for(Int_t j = 0; j < fDigitArrayPtr->GetEntries(); j++)
81     {
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 }
93
94 Int_t 
95 AliHLTPHOSClusterizer::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 }
104   
105
106 Int_t 
107 AliHLTPHOSClusterizer::ClusterizeEvent()
108 {
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
111   Int_t nRecPoints = 0;
112   UInt_t i = 0;
113
114
115   AliHLTPHOSRecPointDataStruct *recPoint = 0;
116   
117   //Clusterization starts
118   for(i = 0; i < fDigitContainerPtr->fNDigits; i++)
119     { 
120    
121       fDigitsInCluster = 0;
122       
123       if(fDigitContainerPtr->fDigitDataStruct[i].fEnergy < fEmcClusteringThreshold)
124         {
125           continue;
126         }
127
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;
141   
142   return nRecPoints;
143 }
144
145 void
146 AliHLTPHOSClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTPHOSRecPointDataStruct* recPoint)
147
148 {
149
150   for(UInt_t j = 0; j < fDigitContainerPtr->fNDigits; j++)
151     {
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             }
169         }
170     } 
171   return;
172 }
173
174
175 Int_t 
176 AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1, 
177                                             AliHLTPHOSDigitDataStruct* digit2)
178 {
179   
180   //Int_t coord1[4];
181   //Int_t coord2[4];
182   
183   // fPHOSGeometry->AbsToRelNumbering(digit1->fID, coord1);
184   //fPHOSGeometry->AbsToRelNumbering(digit2->fID, coord2);
185   
186
187   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
188     { 
189
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 ))
194         {
195           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
196             return 1; 
197         }
198     }
199   return 0;
200 }
201
202 void
203 AliHLTPHOSClusterizer::CalculateCenterOfGravity()
204 {
205
206   // Calculates the center of gravity in the local PHOS-module coordinates 
207   Float_t wtot = 0.;
208  
209   Int_t relid[4];
210
211   Float_t x = 0.;
212   Float_t z = 0.;
213   Float_t xi = 0.;
214   Float_t zi = 0.;
215
216   AliHLTPHOSRecPointDataStruct *recPoint = 0;
217   AliHLTPHOSDigitDataStruct *digit = 0;
218
219   AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
220
221   Int_t iDigit = 0;
222   Int_t iRecPoint = 0;
223
224   for(iRecPoint=0; iRecPoint<fRecPointContainerPtr->fNRecPoints; iRecPoint++) 
225     {
226       recPoint = &(fRecPointContainerPtr->fRecPointArray[iRecPoint]);
227       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
228         {
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;
253         }
254     }
255       
256 }
257
258