new classes for cosmic calibration
[u/mrichter/AliRoot.git] / HLT / PHOS / AliHLTPHOSClusterizer.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Authors: Øystein Djuvsland <oysteind@ift.uib.no>                       *
5  *                                                                        *
6  * Permission to use, copy, modify and distribute this software and its   *
7  * documentation strictly for non-commercial purposes is hereby granted   *
8  * without fee, provided that the above copyright notice appears in all   *
9  * copies and that both the copyright notice and this permission notice   *
10  * appear in the supporting documentation. The authors make no claims     *
11  * about the suitability of this software for any purpose. It is          *
12  * provided "as is" without express or implied warranty.                  *
13  **************************************************************************/
14
15  
16 /** @file   AliHLTPHOSClusterizer.cxx
17     @author Øystein Djuvsland
18     @date   
19     @brief  A temporary clusterizer for PHOS
20 */
21
22
23
24 //#include "AliHLTPHOSDefinitions.h"
25 #include "AliHLTPHOSClusterizer.h"
26 //#include "AliHLTPHOSCommonDefs.h"
27 #include "TVector3.h"
28 #include "TMath.h"
29 #include "AliHLTPHOSRcuCellEnergyDataStruct.h"
30 #include "AliHLTPHOSRecPointListDataStruct.h"
31 #include "AliHLTPHOSValidCellDataStruct.h"
32 #include "AliHLTPHOSRecPointDataStruct.h"
33 #include "AliHLTPHOSClusterDataStruct.h"
34 //#include "AliHLTPHOSConstants.h"
35
36 //using namespace PhosHLTConst;
37
38
39
40 ClassImp(AliHLTPHOSClusterizer)
41
42 /**
43 * Main constructor
44 **/
45 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():AliHLTPHOSBase(),fPHOSModule(0), fThreshold(0), fClusterThreshold(0), 
46                                                fHighGainFactor(0.005), fLowGainFactor(0.08),
47                                                fArraySize(3), fMultiplicity(fArraySize*fArraySize)
48 //AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():AliHLTPHOSProcessor(),fPHOSModule(0), fThreshold(0), fClusterThreshold(0), 
49 //                                             fHighGainFactor(0.005), fLowGainFactor(0.08),
50 //                                             fArraySize(3), fMultiplicity(fArraySize*fArraySize)
51 {
52   //See header file for documentation
53   
54 }//end
55
56 // PTH AliHLTPHOSClusterizer::AliHLTPHOSClusterizer(const AliHLTPHOSClusterizer &):AliHLTPHOSBase(),fPHOSModule(0), fThreshold(0), fClusterThreshold(0), 
57 //                                                                          fHighGainFactor(0.005), fLowGainFactor(0.08),
58 //                                                                          fArraySize(3), fMultiplicity(fArraySize*fArraySize)
59
60 //{
61   //Copy constructor, not implemented
62 //}//end
63
64 AliHLTPHOSClusterizer:: ~AliHLTPHOSClusterizer()  
65 {
66   //Destructor
67 }
68
69 /**
70 * Building a 2D array of cell energies of the PHOS detector
71 * @param cellData object containing the cell energies from one event
72 * @param recPointList list to be filled with coordinates of local maxima in the detector
73 **/
74
75 int
76 AliHLTPHOSClusterizer::BuildCellEnergyArray(AliHLTPHOSRcuCellEnergyDataStruct* cellData, 
77                                             AliHLTPHOSRecPointListDataStruct* recPointList)
78 {
79   //BuildCellEnergyArray
80
81   Int_t x = 0;
82   Int_t z = 0;
83   Int_t gain = 0;
84   Int_t xMod = 0;
85   Int_t zMod = 0;
86   Int_t index = 0;
87   Double_t energyCount = 0;
88
89   for(Int_t cell = 0; cell < cellData->fCnt; cell++)
90     {
91      
92       z = (cellData->fValidData[cell]).fZ;
93       x = (cellData->fValidData[cell]).fX;
94       gain = (cellData->fValidData[cell]).fGain;
95
96       zMod = z + (cellData->fRcuZ)*N_ZROWS_RCU;
97       xMod = x + (cellData->fRcuX)*N_XCOLUMNS_RCU;
98       
99       energyCount = (cellData->fValidData[cell]).fEnergy;
100
101       if(gain == 0 && energyCount < 1023) 
102         {
103           fEnergyArray[xMod][zMod] = fHighGainFactor * energyCount;
104          
105           if(fEnergyArray[xMod][zMod] > fClusterThreshold)
106             {         
107               recPointList[index].fZ = zMod;
108               recPointList[index].fX = xMod;
109               
110               for(Int_t j = 0; j < index; j++)
111                 {
112                   if(recPointList[j].fZ == zMod && recPointList[j].fX == xMod)
113                     recPointList[j].fZ = -1;
114                 }
115               index++;
116             }
117           
118         }
119       else if(gain == 1 && fEnergyArray[xMod][zMod] == 0) 
120         {
121           fEnergyArray[xMod][zMod] = fLowGainFactor * energyCount;
122            if(fEnergyArray[xMod][zMod] > fClusterThreshold)
123             {         
124               recPointList[index].fZ = zMod;
125               recPointList[index].fX = xMod;
126               recPointList[index].fModule = cellData->fModuleID;
127               index++;
128             }
129         }
130
131     }  
132   
133   fPHOSModule = cellData->fModuleID;
134
135   return index;
136 }//end BuildCellEnergyArray
137
138
139
140 /**
141 * Creating an array of rec points
142 * @param recPointStructArrayPtr array to store the rec points
143 * @param list list of rec points
144 * @param nPoints number of points
145 **/
146 int 
147 AliHLTPHOSClusterizer::CreateRecPointStructArray(AliHLTPHOSRecPointDataStruct* recPointStructArrayPtr, 
148                                                  AliHLTPHOSRecPointListDataStruct* list, 
149                                                  int nPoints) 
150
151 {
152   //CreateRecPointStructArray
153
154   Int_t flag = 0;
155   Int_t edgeFlagRows = 0;
156   Int_t edgeFlagCols = 0;
157   Int_t k = 0;
158   Int_t nRecPoints = 0;
159   Int_t z = 0;
160   Int_t x = 0;
161   Int_t i = 0;
162   Int_t j = 0;
163   Int_t a = 0;
164
165   Float_t* energiesList = NULL;
166
167   for(Int_t point = 0; point < nPoints; point++)
168     {
169       z = list[point].fZ;
170       x = list[point].fX;
171       if(z == -1) continue;
172
173       if((z-fArraySize/2) < 0/*= - N_ZROWS_MOD/2*/ || (z+fArraySize/2) >= N_ZROWS_MOD/*/2*/) 
174         {
175           edgeFlagRows = 1;
176           continue;
177         }
178       if((x-fArraySize/2) < 0/*= - N_COLUMNS_MOD/2*/ || (x+fArraySize/2) >= N_XCOLUMNS_MOD) 
179         {
180           edgeFlagCols = 1;
181           continue;
182         }
183       
184
185       if(!flag) energiesList = new Float_t[fMultiplicity];
186
187       flag = 0;
188       k = 0;          
189           
190       for(i = -fArraySize/2; i <= fArraySize/2; i++)
191         {
192           if(flag) break;
193           for(j = -fArraySize/2; j <= fArraySize/2; j++)
194             {                     
195               
196               if(fEnergyArray[x+i][z+j] > fEnergyArray[x][z] && abs(i) < 2 && abs(j) < 2)
197                 {
198                   flag = 1;
199                   break;
200                 }
201              
202               energiesList[k] = fEnergyArray[x+i][z+j];
203               k++;
204             }
205         }
206           
207       
208       if(!flag && k)
209         {
210           for(a = 0; a < k; a++) 
211             {
212               (recPointStructArrayPtr[nRecPoints].fEnergiesListPtr)[a] = energiesList[a];
213             }
214           recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[0] = x;
215           recPointStructArrayPtr[nRecPoints].fCoordinatesPtr[1] = z;
216           recPointStructArrayPtr[nRecPoints].fX = x;
217           recPointStructArrayPtr[nRecPoints].fZ = z;
218           // recPointStructArrayPtr[nRecPoints].fMultiplicity = fMultiplicity;
219           recPointStructArrayPtr[nRecPoints].fPHOSModule = list[point].fModule;
220           nRecPoints++;
221         }
222     }
223   if(energiesList)
224     {
225       delete [] energiesList;
226       energiesList = NULL;
227     }
228
229   return nRecPoints;
230   
231 }//end CreateRecPointStructArray
232
233
234 /**
235 * Calculating the center of gravity of a rec point
236 * Not working well at this point!
237 * @param recPointPtr pointer to the rec point
238 **/
239 int
240 AliHLTPHOSClusterizer::CalculateCenterOfGravity(AliHLTPHOSRecPointDataStruct* recPointPtr)
241 {
242   //CalculateCenterOfGravity
243   //Copied from offline code, and modified
244
245   Float_t xt = 0;
246   Float_t zt = 0;
247   Float_t xi = 0;
248   Float_t zj = 0;
249   Float_t w = 0;
250   Float_t w0 = 4.5;
251   Float_t wtot = 0;
252
253   Int_t x = recPointPtr->fCoordinatesPtr[0];
254   Int_t z = recPointPtr->fCoordinatesPtr[1];
255   
256   for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
257     {
258       for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
259         {                         
260           xi = x + i;
261           zj = z + j;
262           w = TMath::Max( (Float_t)0., (Float_t)(w0 + log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]) ) ) ;
263           //printf("log in cog: %f\n", log( fEnergyArray[x+i][z+j] / fEnergyArray[x][z]));
264           xt += xi * w ;
265           zt += zj * w ;
266           wtot += w ;     
267         }
268     }
269   xt /= wtot ;
270   zt /= wtot ;
271   
272   //printf("wtot cog: %f\n", wtot);
273
274   recPointPtr->fX = xt;
275   recPointPtr->fZ = zt;
276   
277   return 0;
278
279 }//end CalculateCenterOfGravity
280   
281 int
282 AliHLTPHOSClusterizer::CalculateMoments(AliHLTPHOSRecPointDataStruct* recPointPtr, Bool_t axisOnly)
283 {
284   //Copied from offline code, and modified
285   
286   
287   // Calculate the shower moments in the eigen reference system
288   // M2x, M2z, M3x, M4z
289   // Calculate the angle between the shower position vector and the eigen vector
290   
291   Double_t wtot = 0. ;
292   Double_t x    = 0.;
293   Double_t z    = 0.;
294   Double_t dxx  = 0.;
295   Double_t dzz  = 0.;
296   Double_t dxz  = 0.;
297   Double_t lambda0=0, lambda1=0;
298   Float_t logWeight = 4.5;
299   
300   //Int_t iDigit;
301   
302   // 1) Find covariance matrix elements:
303   //    || dxx dxz ||
304   //    || dxz dzz ||
305   
306   Int_t xi;
307   Int_t zj;
308   Double_t w;
309   
310   Int_t xc = recPointPtr->fCoordinatesPtr[0];
311   Int_t zc = recPointPtr->fCoordinatesPtr[1];
312   
313   for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
314     {
315       xi = xc + i;
316       for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
317         {
318           zj = zc + j;
319           if (fEnergyArray[xi][zj] > 0) 
320             {
321
322               //              w     = TMath::Max(0.,logWeight + log( fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
323
324               //printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
325               x    += w * xi ;
326               z    += w * zj ; 
327               dxx  += w * xi * xi ;
328               dzz  += w * zj * zj ;
329               dxz  += w * xi * zj ; 
330               wtot += w ;
331             }
332         }
333     }
334   //printf("wtot: %f\n", wtot);
335   if (wtot>0) 
336     {
337       x   /= wtot ;
338       z   /= wtot ;
339       dxx /= wtot ;
340       dzz /= wtot ;
341       dxz /= wtot ;
342       dxx -= x * x ;
343       dzz -= z * z ;
344       dxz -= x * z ;
345       
346   // 2) Find covariance matrix eigen values lambda0 and lambda1
347
348       lambda0 =  0.5 * (dxx + dzz) + TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
349       lambda1 =  0.5 * (dxx + dzz) - TMath::Sqrt( 0.25 * (dxx - dzz) * (dxx - dzz) + dxz * dxz )  ;
350     }
351     
352   recPointPtr->fM2x   = lambda0;
353   recPointPtr->fM2z   = lambda1;
354   //  printf("lambda0 = %f -- lambda1 = %f\n", lambda0, lambda1);
355   // 3) Find covariance matrix eigen vectors e0 and e1
356   if(!axisOnly)
357     {
358       TVector2 e0,e1;
359       if (dxz != 0)
360         e0.Set(1.,(lambda0-dxx)/dxz);
361       else
362         e0.Set(0.,1.);
363       
364       e0 = e0.Unit();
365       e1.Set(-e0.Y(),e0.X());
366       
367       // 4) Rotate cluster tensor from (x,z) to (e0,e1) system
368       //    and calculate moments M3x and M4z
369       
370       Float_t cosPhi = e0.X();
371       Float_t sinPhi = e0.Y();
372       
373       Float_t xiPHOS ;
374       Float_t zjPHOS ;
375       Double_t dx3, dz3, dz4;
376       wtot = 0.;
377       x    = 0.;
378       z    = 0.;
379       dxx  = 0.;
380       dzz  = 0.;
381       dxz  = 0.;
382       dx3  = 0.;
383       dz3  = 0.;
384       dz4  = 0.;
385       for(Int_t i = -fArraySize/2; i <= fArraySize/2; i++)
386         {
387           for(Int_t j = -fArraySize/2; j <= fArraySize/2; j++)
388             {
389               xi = xc + i;
390               zj = zc + j;
391               xiPHOS = (Float_t)xi*cosPhi + (Float_t)zj*sinPhi;
392               zjPHOS = (Float_t)zj*cosPhi - (Float_t)xi*sinPhi;
393               //              xi = (Float_t)xi*cosPhi + (Float_t)zj*sinPhi;
394               //              zj = (Float_t)zj*cosPhi - (Float_t)xi*sinPhi;
395               if (fEnergyArray[xi][zj] > 0) 
396                 {
397                   w     = TMath::Max(0.,logWeight+TMath::Log(fEnergyArray[xi][zj]/fEnergyArray[xc][zc] ) ) ;
398                   //  printf("log in mom: %f\n", TMath::Log( fEnergyArray[xi][zj] / fEnergyArray[xc][zc]));
399                           
400                   x    += w * xi ;
401                   z    += w * zj ; 
402                   dxx  += w * xi * xi ;
403                   dzz  += w * zj * zj ;
404                   dxz  += w * xi * zj ; 
405                   dx3  += w * xi * xi * xi;
406                   dz3  += w * zj * zj * zj ;
407                   dz4  += w * zj * zj * zj * zj ;
408                   wtot += w ;
409                 }
410             }
411         }
412       if (wtot>0) 
413         {
414           x   /= wtot ;
415           z   /= wtot ;
416           dxx /= wtot ;
417           dzz /= wtot ;
418           dxz /= wtot ;
419           dx3 /= wtot ;
420           dz3 /= wtot ;
421           dz4 /= wtot ;
422           dx3 += -3*dxx*x + 2*x*x*x;
423           dz4 += -4*dz3*z + 6*dzz*z*z -3*z*z*z*z;
424           dxx -= x * x ;
425           dzz -= z * z ;
426           dxz -= x * z ;
427         }
428   
429       // 5) Find an angle between cluster center vector and eigen vector e0
430       
431       Float_t phi = 0;
432       if ( (x*x+z*z) > 0 ) 
433         phi = TMath::ACos ((x*e0.X() + z*e0.Y()) / sqrt(x*x + z*z));
434
435       recPointPtr->fM3x   = dx3;
436       recPointPtr->fM4z   = dz4;
437       recPointPtr->fPhixe = phi;
438       //      printf("%f\n",phi);
439     }
440   return 0;
441 }
442
443
444 /**
445 * Create a simpler data struct for a rec point, containing
446 * only coordinates, energy and timing
447 * @param recPointPtr pointer to the rec point
448 * @param clusterStructPtr pointer to the cluster struct
449 **/
450 int
451 AliHLTPHOSClusterizer::ClusterizeStruct(AliHLTPHOSRecPointDataStruct* recPointPtr, AliHLTPHOSClusterDataStruct* clusterStructPtr)
452 {
453   //ClusterizeStruct
454
455   Float_t clusterEnergy = 0;
456   Float_t* energiesListPtr = recPointPtr->fEnergiesListPtr;
457   
458   for(Int_t i = 0; i < recPointPtr->fMultiplicity; i++)
459     {
460       clusterEnergy += energiesListPtr[i];
461     }
462
463   clusterStructPtr->fLocalPositionPtr[0] = recPointPtr->fX;
464   clusterStructPtr->fLocalPositionPtr[1] = recPointPtr->fZ;
465   clusterStructPtr->fClusterEnergy = clusterEnergy;
466   clusterStructPtr->fPHOSModule = recPointPtr->fPHOSModule;
467
468   return 0;
469
470 }//end ClusterizeStruct
471
472
473
474
475 /**
476 * Resets the cell energy array
477 **/ 
478 int 
479 AliHLTPHOSClusterizer::ResetCellEnergyArray()
480
481 {
482   //ResetCellEnergyArray
483
484   for(Int_t x = 0; x < N_ZROWS_MOD; x++)
485     {
486       for(Int_t z = 0; z < N_XCOLUMNS_MOD; z++)
487         {
488           fEnergyArray[x][z] = 0;
489         }
490     }
491
492   return 0;
493
494 }//end ResetCellEnergyArray
495