added switch for doc to configure.ac; added monilithic doc build; configure.ac and...
[u/mrichter/AliRoot.git] / HLT / PHOS / AliHLTPHOSClusterizer.cxx
1 /**************************************************************************
2  * This file is property of and copyright by the ALICE HLT Project        * 
3  * All rights reserved.                                                   *
4  *                                                                        *
5  * Primary Authors: Oystein Djuvsland                                     *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          * 
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /** 
17  * @file   AliHLTPHOSClusterizer.cxx
18  * @author Oystein Djuvsland
19  * @date 
20  * @brief  Clusterizer for PHOS HLT 
21  */
22
23 // see header file for class documentation
24 // or
25 // refer to README to build package
26 // or
27 // visit http://web.ift.uib.no/~kjeks/doc/alice-hlt
28
29 #include "AliHLTPHOSClusterizer.h"
30 #include "AliHLTPHOSBase.h"
31 #include "TMath.h"
32 #include "AliHLTPHOSRecPointContainerStruct.h"
33 #include "AliHLTPHOSRecPointDataStruct.h"
34 #include "AliHLTPHOSDigitDataStruct.h"
35 #include "AliHLTPHOSDigitContainerDataStruct.h"
36 #include "TClonesArray.h"
37 #include "AliPHOSGeometry.h"
38 #include "AliPHOSDigit.h"
39 #ifndef HAVE_NOT_PHOSRECOPARAMEMC // set from configure if EMC functionality not available in AliPHOSRecoParam
40 #include "AliPHOSRecoParam.h"
41 #else
42 #include "AliPHOSRecoParamEmc.h"
43 #endif
44
45 ClassImp(AliHLTPHOSClusterizer);
46
47 AliHLTPHOSClusterizer::AliHLTPHOSClusterizer():
48   AliHLTPHOSBase(),
49   fEmcClusteringThreshold(0),
50   fEmcMinEnergyThreshold(0),
51   fEmcTimeGate(0),
52   fLogWeight(0),
53   fDigitsInCluster(0),
54   fOnlineMode(true),
55   fDigitArrayPtr(0),
56   fEmcRecPointsPtr(0),
57   fDigitPtr(0),
58   fDigitContainerPtr(0),
59   fRecPointContainerPtr(0),
60   fPHOSGeometry(0),
61   fGetterPtr(0)
62 {
63   //See header file for documentation
64   fPHOSGeometry = AliPHOSGeometry::GetInstance("noCPV");
65   fEmcClusteringThreshold = 0.2;
66   fEmcMinEnergyThreshold = 0.03;
67   fEmcTimeGate = 1.e-6 ;
68   fLogWeight = 4.5;
69   
70
71 }//end
72
73
74 AliHLTPHOSClusterizer::~AliHLTPHOSClusterizer()  
75 {
76   //See header file for documentation
77 }
78
79 void 
80 AliHLTPHOSClusterizer::SetRecPointContainer(AliHLTPHOSRecPointContainerStruct* recPointContainerPtr)
81   { 
82     fRecPointContainerPtr = recPointContainerPtr; 
83     fRecPointContainerPtr->fNRecPoints = 0;
84   }
85
86 void
87 AliHLTPHOSClusterizer::SetRecoParameters(AliPHOSRecoParam* params)
88 {
89   //see header file for documentation
90 #ifndef HAVE_NOT_PHOSRECOPARAMEMC // set from configure if EMC functionality not available in AliPHOSRecoParam
91   // the new AliPHOSRecoParam functions, available from revision
92   fEmcClusteringThreshold = params->GetEMCClusteringThreshold();
93   fEmcMinEnergyThreshold = params->GetEMCMinE();
94   fLogWeight = params->GetEMCLogWeight();
95 #else
96   fEmcClusteringThreshold = params->GetClusteringThreshold();
97   fEmcMinEnergyThreshold = params->GetMinE();
98   fLogWeight = params->GetLogWeight();
99 #endif
100 }
101
102 void
103 AliHLTPHOSClusterizer::SetOfflineMode(AliPHOSLoader* getter)
104 {
105   //see header file for documentation
106   fRecPointContainerPtr = new AliHLTPHOSRecPointContainerStruct();
107   fDigitContainerPtr = new AliHLTPHOSDigitContainerDataStruct();
108   fGetterPtr = getter;
109   fDigitArrayPtr = fGetterPtr->Digits();
110   fEmcRecPointsPtr = fGetterPtr->EmcRecPoints();
111   fOnlineMode = false;
112 }
113
114 Int_t
115 AliHLTPHOSClusterizer::GetEvent(Int_t i)
116 {
117   //see header file for documentation
118   Int_t coord[4];
119
120   //  fGetterPtr->Event(i, "D");
121   fGetterPtr->GetEvent();
122   for(Int_t j = 0; j < fDigitArrayPtr->GetEntries(); j++)
123     {
124       fDigitPtr = (AliPHOSDigit*)fDigitArrayPtr->At(j);
125       fPHOSGeometry->AbsToRelNumbering(fDigitPtr->GetId(), coord);
126       fDigitContainerPtr->fDigitDataStruct[j].fX = coord[3];
127       fDigitContainerPtr->fDigitDataStruct[j].fZ = coord[2];
128       fDigitContainerPtr->fDigitDataStruct[j].fModule = coord[0];
129       fDigitContainerPtr->fDigitDataStruct[j].fEnergy = fDigitPtr->GetEnergy();
130       fDigitContainerPtr->fDigitDataStruct[j].fTime = fDigitPtr->GetTime();
131     }
132   fDigitContainerPtr->fNDigits = fDigitArrayPtr->GetEntriesFast();
133   return 0;
134 }
135
136 Int_t 
137 AliHLTPHOSClusterizer::GetNEvents()
138 {
139   //see header file for documentation
140   if(fOnlineMode)
141     {
142       //     Logging(kHLTLogWarning, __FILE__ , "information not available" , "GetNEvents()  Number of events not available in online mod");
143       return -1;
144     }
145   return fGetterPtr->MaxEvent();
146 }
147   
148
149 Int_t 
150 AliHLTPHOSClusterizer::ClusterizeEvent()
151 {
152   //see header file for documentation
153   Int_t nRecPoints = 0;
154   UInt_t i = 0;
155
156   AliHLTPHOSRecPointDataStruct *recPoint = 0;
157
158   //Clusterization starts
159   for(i = 0; i < fDigitContainerPtr->fNDigits; i++)
160     { 
161       
162       fDigitsInCluster = 0;
163       if(fDigitContainerPtr->fDigitDataStruct[i].fEnergy < fEmcClusteringThreshold)
164         {
165           continue;
166         }
167       recPoint = &(fRecPointContainerPtr->fRecPointArray[nRecPoints]);
168       recPoint->fAmp = 0;
169       //TODO!!!!!!!
170       recPoint->fModule = fDigitContainerPtr->fDigitDataStruct[i].fModule;
171       //TODO!!!!!!!
172       //recPoint->
173       recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[i];
174       recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[i].fEnergy;
175       fDigitContainerPtr->fDigitDataStruct[i].fEnergy = 0;
176       fDigitsInCluster++;
177       nRecPoints++;
178       ScanForNeighbourDigits(i, recPoint);
179
180       recPoint->fMultiplicity = fDigitsInCluster;
181       
182     }//end of clusterization
183   fRecPointContainerPtr->fNRecPoints = nRecPoints;
184   
185   return nRecPoints;
186 }
187
188 void
189 AliHLTPHOSClusterizer::ScanForNeighbourDigits(Int_t index, AliHLTPHOSRecPointDataStruct* recPoint)
190
191 {
192   //see header file for documentation
193
194   for(UInt_t j = 0; j < fDigitContainerPtr->fNDigits; j++)
195     {
196       if(fDigitContainerPtr->fDigitDataStruct[j].fEnergy > fEmcMinEnergyThreshold)
197         {
198           switch(AreNeighbours(&(fDigitContainerPtr->fDigitDataStruct[index]),
199                                &(fDigitContainerPtr->fDigitDataStruct[j])))
200             {
201             case 0:
202               break; 
203             case 1:      
204               recPoint->fDigitsList[fDigitsInCluster] =  fDigitContainerPtr->fDigitDataStruct[j];
205               recPoint->fAmp += fDigitContainerPtr->fDigitDataStruct[j].fEnergy;
206               fDigitContainerPtr->fDigitDataStruct[j].fEnergy = 0;
207               fDigitsInCluster++;
208               ScanForNeighbourDigits(j, recPoint);
209               break;
210             case 2:
211               break;
212             }
213         }
214     } 
215   return;
216 }
217
218
219 Int_t 
220 AliHLTPHOSClusterizer::AreNeighbours(AliHLTPHOSDigitDataStruct* digit1, 
221                                             AliHLTPHOSDigitDataStruct* digit2)
222 {
223   //see header file for documentation
224
225   if ( (digit1->fModule == digit2->fModule) /*&& (coord1[1]==coord2[1])*/ ) // inside the same PHOS module
226     { 
227       Int_t rowdiff = TMath::Abs( digit1->fZ - digit2->fZ );  
228       Int_t coldiff = TMath::Abs( digit1->fX - digit2->fX ); 
229           
230       if (( coldiff <= 1 )  && ( rowdiff <= 1 ))
231         {
232           if(TMath::Abs(digit1->fTime - digit2->fTime ) < fEmcTimeGate)
233             {
234               return 1; 
235             }
236         }
237     }
238   return 0;
239 }
240
241 void
242 AliHLTPHOSClusterizer::CalculateCenterOfGravity()
243 {
244   //see header file for documentation
245   Float_t wtot = 0.;
246   Float_t x = 0.;
247   Float_t z = 0.;
248   Float_t xi = 0.;
249   Float_t zi = 0.;
250
251   AliHLTPHOSRecPointDataStruct *recPoint = 0;
252   AliHLTPHOSDigitDataStruct *digit = 0;
253   UInt_t iDigit = 0;
254   UInt_t iRecPoint = 0;
255
256   for(iRecPoint=0; iRecPoint<fRecPointContainerPtr->fNRecPoints; iRecPoint++) 
257     {
258       recPoint = &(fRecPointContainerPtr->fRecPointArray[iRecPoint]);
259       for(iDigit = 0; iDigit < recPoint->fMultiplicity; iDigit++)
260         {
261           digit = &(recPoint->fDigitsList[iDigit]);
262           
263           //fPHOSGeometry->AbsToRelNumbering(digit->fID, relid) ;
264           //  fPHOSGeometry->RelPosInModule(relid, xi, zi);
265           xi = digit->fX;
266           zi = digit->fZ;
267           
268           if (recPoint->fAmp > 0 && digit->fEnergy > 0) 
269             {
270               Float_t w = TMath::Max( 0., fLogWeight + TMath::Log( digit->fEnergy / recPoint->fAmp ) ) ;
271               x    += xi * w ;
272               z    += zi * w ;
273               wtot += w ;
274             }
275         }
276       
277       if (wtot>0) 
278         {
279           recPoint->fX = x/wtot ;
280           recPoint->fZ = z/wtot ;
281         }
282       else
283         {
284           recPoint->fAmp = 0;
285         }
286     }
287       
288 }
289
290