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