- moving the cluster parameter evalution from own component to clusterizer
[u/mrichter/AliRoot.git] / HLT / CALO / AliHLTCaloClusterAnalyser.cxx
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 "AliHLTCaloClusterizer.h"
41
42 ClassImp(AliHLTCaloClusterAnalyser);
43
44 AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() : 
45   //  AliHLTCaloBase(),
46   fLogWeight(4.5),
47   fRecPointArray(0),
48   fDigitDataArray(0),
49   fNRecPoints(0),
50   fCaloClusterDataPtr(0),
51   fCaloClusterHeaderPtr(0),
52   //fAnalyzerPtr(0),
53   fDoClusterFit(false),
54   fHaveCPVInfo(false),
55   fDoPID(false),
56   fHaveDistanceToBadChannel(false),
57   fGeometry(0)
58 {
59   //See header file for documentation
60 }
61
62 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser() 
63 {
64 }
65
66 void 
67 AliHLTCaloClusterAnalyser::SetCaloClusterData(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
68
69   //see header file for documentation
70   fCaloClusterDataPtr = caloClusterDataPtr; 
71 }
72
73 void
74 AliHLTCaloClusterAnalyser::SetRecPointArray(AliHLTCaloRecPointDataStruct **recPointDataPtr, Int_t nRecPoints)
75
76   fRecPointArray = recPointDataPtr; 
77   fNRecPoints = nRecPoints;
78 }
79
80 void 
81 AliHLTCaloClusterAnalyser::SetDigitDataArray(AliHLTCaloDigitDataStruct *digits) 
82
83 //   AliHLTCaloClusterizer cl("PHOS");
84   // cl.CheckDigits(fRecPointArray, digits, fNRecPoints);
85    fDigitDataArray = digits; 
86    //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
87 }
88
89 Int_t
90 AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
91 {
92   //see header file for documentation
93   Float_t wtot = 0.;
94   Float_t x = 0.;
95   Float_t z = 0.;
96   Float_t xi = 0.;
97   Float_t zi = 0.;
98
99   AliHLTCaloDigitDataStruct *digit = 0;
100   //AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
101    
102   
103   //AliHLTCaloClusterizer       cl("PHOS");
104   //cl.CheckDigits(fRecPointArray, fDigitDataArray, fNRecPoints);
105
106   UInt_t iDigit = 0;
107
108   for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++) 
109     {
110       AliHLTCaloRecPointDataStruct *recPoint = fRecPointArray[iRecPoint];
111       //      digit = &(recPoint->fDigits);
112       cout << "COG: Rec point multiplicity: " << recPoint->fMultiplicity << ", rec point energy: " << recPoint->fAmp << endl;
113       Int_t *digitIndexPtr = &(recPoint->fDigits);
114
115       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
116         {
117           //cout << "1. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;
118
119           digit = &(fDigitDataArray[*digitIndexPtr]);
120           cout << "COG: 2. Digit array: " << fDigitDataArray << ", Digit: " << digit << ", Index: " << *digitIndexPtr << ", index pointer: " << digitIndexPtr<<  endl;
121           cout << ", dig Energy: " << digit->fEnergy << ", Index: " << *digitIndexPtr << ", ID: " << digit->fID << endl;        
122           xi = digit->fX;
123           zi = digit->fZ;
124           //  cout << "COG digits (x:z:E:time): " << xi << " : " << zi << " : " << digit->fEnergy << " : " << digit->fTime << endl;
125           if (recPoint->fAmp > 0 && digit->fEnergy > 0) 
126             {
127               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
128               x    += xi * w ;
129               z    += zi * w ;
130               wtot += w ;
131             }
132           digitIndexPtr++;
133         }
134       //cout << endl;
135       if (wtot>0) 
136         {
137           recPoint->fX = x/wtot ;
138           recPoint->fZ = z/wtot ;
139         }
140       else
141         {
142           recPoint->fAmp = 0;
143         }
144     }
145   return 0;
146 }
147
148
149 Int_t 
150 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
151 {
152   //See header file for documentation
153   return 0;
154 }
155
156 Int_t 
157 AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
158 {
159   //See header file for documentation
160   return 0;
161 }
162
163
164 Int_t 
165 AliHLTCaloClusterAnalyser::DeconvoluteClusters()
166 {
167   //See header file for documentation
168   return 0;
169 }
170
171 Int_t 
172 AliHLTCaloClusterAnalyser::CreateClusters(Int_t nRecPoints, UInt_t availableSize, UInt_t& totSize)
173 {
174   //See header file for documentation
175
176   fNRecPoints = nRecPoints;
177
178   if(fGeometry == 0)
179   {
180      HLTError("No geometry object is initialised, creation of clusters stopped");
181   }
182
183   CalculateCenterOfGravity();
184
185   UInt_t maxClusterSize = sizeof(AliHLTCaloClusterDataStruct) + (6 << 7); //Reasonable estimate... (6 = sizeof(Short_t) + sizeof(Float_t)
186
187   //  AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);  
188   AliHLTCaloDigitDataStruct* digitPtr = 0;
189
190   AliHLTCaloClusterDataStruct* caloClusterPtr = fCaloClusterDataPtr;
191   UShort_t* cellIDPtr = &(caloClusterPtr->fCellsAbsId);
192   Float_t* cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
193   
194   Int_t id = -1;
195   TVector3 globalPos;
196
197   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
198     {
199 //      if(availableSize
200       AliHLTCaloRecPointDataStruct *recPointPtr = fRecPointArray[i];
201       cout << "CA: rec point energy: " << recPointPtr->fAmp << ", rec point multiplicity: " << recPointPtr->fMultiplicity << endl;
202       if(availableSize < (totSize + maxClusterSize)) 
203         {
204           return -1; //Might get out of buffer, exiting
205         }
206       //      cout << "Local Position (x:z:module): " << recPointPtr->fX << " : "<< recPointPtr->fZ << " : " << recPointPtr->fModule << endl;
207       
208       AliHLTCaloGlobalCoordinate globalCoord;
209       fGeometry->GetGlobalCoordinates(*recPointPtr, globalCoord);
210       // cout << "Global Position (x:y:z): " << globalPos[0] << " : "<< globalPos[1] << " : " << globalPos[2] << endl << endl;
211
212       caloClusterPtr->fGlobalPos[0] = globalCoord.fX;
213       caloClusterPtr->fGlobalPos[1] = globalCoord.fY;
214       caloClusterPtr->fGlobalPos[2] = globalCoord.fZ;
215
216       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
217   
218       cellIDPtr = &(caloClusterPtr->fCellsAbsId);
219       cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
220      
221       for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
222         {
223          // fPHOSGeometry->RelPosToAbsId((Int_t)(recPointPtr->fModule), (double)(digitPtr->fX), (double)(digitPtr->fZ), id);
224 /*        *cellIDPtr = id;
225           *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
226           digitPtr++;
227           cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t)); 
228           cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Short_t)); */
229         }
230
231       caloClusterPtr->fEnergy = recPointPtr->fAmp;
232       cout << "CA: cluster energy: " << caloClusterPtr->fEnergy << endl;
233       cout << "CA: recpoint energy: " << recPointPtr->fAmp << endl;
234       if(fDoClusterFit)
235         {
236           FitCluster(recPointPtr);
237         }
238       else
239         {
240           caloClusterPtr->fDispersion = 0;
241           caloClusterPtr->fFitQuality = 0;
242           caloClusterPtr->fM20 = 0;
243           caloClusterPtr->fM02 = 0;
244
245         }
246       if(fHaveCPVInfo)
247         {
248           caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
249         }
250       else
251         {
252           caloClusterPtr->fEmcCpvDistance = -1;
253         }
254       if(fDoPID)
255         {
256           DoParticleIdentification(caloClusterPtr);
257         }
258       else
259         {
260           for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
261             {
262               caloClusterPtr->fPID[k] = 0;
263             }
264         }
265       if(fHaveDistanceToBadChannel)
266         {
267           caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
268         }
269       else
270         {
271           caloClusterPtr->fDistanceToBadChannel = -1;
272         }
273
274       caloClusterPtr->fClusterType = (AliESDCaloCluster::kPHOSCluster);
275       //      totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);   
276       totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));   
277
278      
279       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
280       caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
281
282       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
283       ///digitPtr = &(recPointPtr->fDigits);  
284     }
285   //  cout << "CA: Energy End: " << fCaloClusterDataPtr->fEnergy << endl;
286   //cout << "CA totSize: " << totSize << endl;
287   return fNRecPoints;
288
289 }
290