d074694d9ad7667740fb0eae1100ae857c6cf727
[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 "AliPHOSGeoUtils.h"
37 #include "AliESDCaloCluster.h"
38 #include "TMath.h"
39 #include "TVector3.h"
40
41 ClassImp(AliHLTCaloClusterAnalyser);
42
43 AliHLTCaloClusterAnalyser::AliHLTCaloClusterAnalyser() :
44   //  AliHLTCaloBase(),
45   fLogWeight(0),
46   fRecPointDataPtr(0),
47   fNRecPoints(0),
48   fCaloClusterDataPtr(0),
49   fCaloClusterHeaderPtr(0),
50   fPHOSGeometry(0),
51   //fAnalyzerPtr(0),
52   fDoClusterFit(false),
53   fHaveCPVInfo(false),
54   fDoPID(false),
55   fHaveDistanceToBadChannel(false)
56 {
57   //See header file for documentation
58   fLogWeight = 4.5;
59
60   //  fAnalyzerPtr = new AliHLTCaloPhysicsAnalyzer();
61   //  fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
62   fPHOSGeometry = new AliPHOSGeoUtils("PHOS", "noCPV");
63 }
64
65 AliHLTCaloClusterAnalyser::~AliHLTCaloClusterAnalyser() 
66 {
67 }
68
69 void 
70 AliHLTCaloClusterAnalyser::SetCaloClusterDataPtr(AliHLTCaloClusterDataStruct *caloClusterDataPtr)
71
72   //see header file for documentation
73   fCaloClusterDataPtr = caloClusterDataPtr; 
74 }
75 void
76 AliHLTCaloClusterAnalyser::SetRecPointDataPtr(AliHLTCaloRecPointHeaderStruct *recPointDataPtr)
77
78   fNRecPoints = recPointDataPtr->fNRecPoints;
79   fRecPointDataPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(reinterpret_cast<Char_t*>(recPointDataPtr)+sizeof(AliHLTCaloRecPointHeaderStruct)); 
80 }
81
82 Int_t
83 AliHLTCaloClusterAnalyser::CalculateCenterOfGravity()
84 {
85   //see header file for documentation
86   Float_t wtot = 0.;
87   Float_t x = 0.;
88   Float_t z = 0.;
89   Float_t xi = 0.;
90   Float_t zi = 0.;
91
92   AliHLTCaloDigitDataStruct *digit = 0;
93   //AliPHOSGeometry * phosgeom =  AliPHOSGeometry::GetInstance() ;
94
95   AliHLTCaloRecPointDataStruct *recPoint = fRecPointDataPtr;
96
97   UInt_t iDigit = 0;
98
99   for(Int_t iRecPoint=0; iRecPoint < fNRecPoints; iRecPoint++) 
100     {
101       //      digit = &(recPoint->fDigits);
102       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
103         {
104           
105           xi = digit->fX;
106           zi = digit->fZ;
107           //  cout << "COG digits (x:z:E:time): " << xi << " : " << zi << " : " << digit->fEnergy << " : " << digit->fTime << endl;
108           if (recPoint->fAmp > 0 && digit->fEnergy > 0) 
109             {
110               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
111               x    += xi * w ;
112               z    += zi * w ;
113               wtot += w ;
114             }
115           digit++;
116         }
117       //cout << endl;
118       if (wtot>0) 
119         {
120           recPoint->fX = x/wtot ;
121           recPoint->fZ = z/wtot ;
122         }
123       else
124         {
125           recPoint->fAmp = 0;
126         }
127       recPoint = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digit);
128     }
129   return 0;
130 }
131
132
133 Int_t 
134 AliHLTCaloClusterAnalyser::CalculateRecPointMoments()
135 {
136   //See header file for documentation
137   return 0;
138 }
139
140 Int_t 
141 AliHLTCaloClusterAnalyser::CalculateClusterMoments(AliHLTCaloRecPointDataStruct */*recPointPtr*/, AliHLTCaloClusterDataStruct* /*clusterPtr*/)
142 {
143   //See header file for documentation
144   return 0;
145 }
146
147
148 Int_t 
149 AliHLTCaloClusterAnalyser::DeconvoluteClusters()
150 {
151   //See header file for documentation
152   return 0;
153 }
154
155 Int_t 
156 AliHLTCaloClusterAnalyser::CreateClusters(UInt_t availableSize, UInt_t& totSize)
157 {
158   //See header file for documentation
159
160   UInt_t maxClusterSize = sizeof(AliHLTCaloClusterDataStruct) + (6 << 7); //Reasonable estimate... (6 = sizeof(Short_t) + sizeof(Float_t)
161
162   AliHLTCaloRecPointDataStruct* recPointPtr = fRecPointDataPtr;
163   //  AliHLTCaloDigitDataStruct* digitPtr = &(recPointPtr->fDigits);  
164   AliHLTCaloDigitDataStruct* digitPtr = 0;
165
166   AliHLTCaloClusterDataStruct* caloClusterPtr = fCaloClusterDataPtr;
167   UShort_t* cellIDPtr = &(caloClusterPtr->fCellsAbsId);
168   Float_t* cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
169   
170   Int_t id = -1;
171   TVector3 globalPos;
172
173   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
174     {
175       
176       if(availableSize < (totSize + maxClusterSize)) 
177         {
178           return -1; //Might get out of buffer, exiting
179         }
180       //      cout << "Local Position (x:z:module): " << recPointPtr->fX << " : "<< recPointPtr->fZ << " : " << recPointPtr->fModule << endl;
181       fPHOSGeometry->Local2Global(recPointPtr->fModule, recPointPtr->fX, recPointPtr->fZ, globalPos);
182       // cout << "Global Position (x:y:z): " << globalPos[0] << " : "<< globalPos[1] << " : " << globalPos[2] << endl << endl;
183
184       caloClusterPtr->fGlobalPos[0] = globalPos[0];
185       caloClusterPtr->fGlobalPos[1] = globalPos[1];
186       caloClusterPtr->fGlobalPos[2] = globalPos[2];
187
188       caloClusterPtr->fNCells = recPointPtr->fMultiplicity;
189   
190       cellIDPtr = &(caloClusterPtr->fCellsAbsId);
191       cellAmpFracPtr = &(caloClusterPtr->fCellsAmpFraction);
192      
193       for(UInt_t j = 0; j < caloClusterPtr->fNCells; j++)
194         {
195           fPHOSGeometry->RelPosToAbsId((Int_t)(recPointPtr->fModule), (double)(digitPtr->fX), (double)(digitPtr->fZ), id);
196           *cellIDPtr = id;
197           *cellAmpFracPtr = digitPtr->fEnergy/recPointPtr->fAmp;
198           digitPtr++;
199           cellIDPtr = reinterpret_cast<UShort_t*>(reinterpret_cast<char*>(cellAmpFracPtr) + sizeof(Float_t)); 
200           cellAmpFracPtr = reinterpret_cast<Float_t*>(reinterpret_cast<char*>(cellIDPtr) + sizeof(Short_t)); 
201         }
202
203       caloClusterPtr->fEnergy = recPointPtr->fAmp;
204
205       if(fDoClusterFit)
206         {
207           FitCluster(recPointPtr);
208         }
209       else
210         {
211           caloClusterPtr->fDispersion = 0;
212           caloClusterPtr->fFitQuality = 0;
213           caloClusterPtr->fM20 = 0;
214           caloClusterPtr->fM02 = 0;
215
216         }
217       if(fHaveCPVInfo)
218         {
219           caloClusterPtr->fEmcCpvDistance = GetCPVDistance(recPointPtr);
220         }
221       else
222         {
223           caloClusterPtr->fEmcCpvDistance = -1;
224         }
225       if(fDoPID)
226         {
227           DoParticleIdentification(caloClusterPtr);
228         }
229       else
230         {
231           for(Int_t k = 0; k < AliPID::kSPECIESN; k++)
232             {
233               caloClusterPtr->fPID[k] = 0;
234             }
235         }
236       if(fHaveDistanceToBadChannel)
237         {
238           caloClusterPtr->fDistanceToBadChannel = GetDistanceToBadChannel(caloClusterPtr);
239         }
240       else
241         {
242           caloClusterPtr->fDistanceToBadChannel = -1;
243         }
244
245       caloClusterPtr->fClusterType = (AliESDCaloCluster::kPHOSCluster);
246       //      totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells)*(sizeof(Short_t) +sizeof(Float_t)-1);   
247       totSize += sizeof(AliHLTCaloClusterDataStruct) + (caloClusterPtr->fNCells-1)*(sizeof(Short_t) + sizeof(Float_t));   
248
249       //      caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellAmpFracPtr);
250       caloClusterPtr = reinterpret_cast<AliHLTCaloClusterDataStruct*>(cellIDPtr);
251       recPointPtr = reinterpret_cast<AliHLTCaloRecPointDataStruct*>(digitPtr);
252       ///digitPtr = &(recPointPtr->fDigits);  
253     }
254   //  cout << "CA: Energy End: " << fCaloClusterDataPtr->fEnergy << endl;
255   //cout << "CA totSize: " << totSize << endl;
256   return fNRecPoints;
257
258 }
259