defects fixed
[u/mrichter/AliRoot.git] / RALICE / AliCalcluster.cxx
CommitLineData
4c039060 1/**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3 * *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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
f531a546 16// $Id$
4c039060 17
959fbac5 18///////////////////////////////////////////////////////////////////////////
19// Class AliCalcluster
20// Description of a cluster of calorimeter modules.
dafe31a2 21// A 2D (matrix) geometry is assumed in which a cluster center is identified
22// by two integer indices (i,j), e.g. row and column indicators.
23//
24// The 1st signal value is the signal of the complete cluster.
25// This is the signal which is provided as default by invoking GetSignal().
26//
27// In case clustering/grouping of module signals was performed over several
28// rings around the center (see e.g. AliCalorimeter::Group), the following
29// additional information is provided by the various signal values :
30//
31// The 2nd signal value is the original signal of the central module.
32// The 3rd signal value is the total signal within the 1st (i.e. 3x3) ring of
33// modules around the cluster center.
34// The 4th signal value is the total signal within the 2nd (i.e. 5x5) ring of
35// modules around the cluster center.
36// Etc....
37//
38// Note : In case the cluster consists of only 1 module, then only the
39// 1st signal value will be present (for obvious reasons).
40//
41// Some dispersion info about cluster topology is provided in order
959fbac5 42// to enable EM or hadronic cluster identification.
43//
44//--- Author: Nick van Eijndhoven 13-jun-1997 UU-SAP Utrecht
f531a546 45//- Modified: NvE $Date$ UU-SAP Utrecht
959fbac5 46///////////////////////////////////////////////////////////////////////////
47
d88f97cc 48#include "AliCalcluster.h"
c72198f1 49#include "Riostream.h"
d88f97cc 50
51ClassImp(AliCalcluster) // Class implementation to enable ROOT I/O
52
c72198f1 53AliCalcluster::AliCalcluster() : AliSignal()
d88f97cc 54{
c72198f1 55// Default constructor, all data is set to 0
56 fRow=0;
57 fCol=0;
d88f97cc 58 fNmods=0;
d88f97cc 59 fRowdisp=0.;
60 fColdisp=0.;
61 fNvetos=0;
62 fVetos=0;
c72198f1 63 SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
d88f97cc 64}
65///////////////////////////////////////////////////////////////////////////
66AliCalcluster::~AliCalcluster()
67{
68// Destructor to delete dynamically allocated memory
69 if (fVetos)
70 {
d88f97cc 71 delete fVetos;
72 fVetos=0;
73 }
74}
75///////////////////////////////////////////////////////////////////////////
261c0caf 76AliCalcluster::AliCalcluster(const AliCalcluster& c) : AliSignal(c)
c72198f1 77{
78// Copy constructor
79 fRow=c.fRow;
80 fCol=c.fCol;
81 fNmods=c.fNmods;
82 fRowdisp=c.fRowdisp;
83 fColdisp=c.fColdisp;
84 fNvetos=c.fNvetos;
85
86 fVetos=new TObjArray();
87 fVetos->SetOwner();
88
89 for (Int_t i=1; i<=fNvetos; i++)
90 {
91 AliSignal* sx=c.GetVetoSignal(i);
92 fVetos->Add(new AliSignal(*sx));
93 }
94}
95///////////////////////////////////////////////////////////////////////////
96AliCalcluster::AliCalcluster(AliCalmodule& m) : AliSignal()
d88f97cc 97{
98// Cluster constructor with module m as center.
f40f8fbd 99// Module data is only entered for a module which contains a signal,
100// has not been used in a cluster yet, and is not declared dead.
101//
102// Note :
103// It is advised NOT to start a cluster with modules situated at a detector edge.
104// This feature is automatically checked when using the built-in clustering
105// of AliCalorimeter.
d88f97cc 106
f40f8fbd 107 Ali3Vector r;
d88f97cc 108
dafe31a2 109 Float_t sig=m.GetClusteredSignal();
110
111 if (sig>0. && m.GetDeadValue()==0)
d88f97cc 112 {
c72198f1 113 fRow=m.GetRow();
114 fCol=m.GetColumn();
f40f8fbd 115 r=m.GetPosition();
116 SetPosition(r);
1fbffa23 117 sig=m.GetSignal(1,1); // Use the gain etc... corrected module signal
dafe31a2 118 SetSignal(sig);
d88f97cc 119 fNmods=1;
d88f97cc 120 fRowdisp=0.;
121 fColdisp=0.;
122 m.SetClusteredSignal(0.); // mark module as used in cluster
123 fNvetos=0;
124 fVetos=0;
125 }
126 else
127 {
c72198f1 128 fRow=0;
129 fCol=0;
f40f8fbd 130 SetPosition(r);
d88f97cc 131 fNmods=0;
d88f97cc 132 fRowdisp=0.;
133 fColdisp=0.;
134 fNvetos=0;
135 fVetos=0;
8e8e6c7f 136 }
c72198f1 137 SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
d88f97cc 138}
139///////////////////////////////////////////////////////////////////////////
261c0caf 140Int_t AliCalcluster::GetRow() const
d88f97cc 141{
142// Provide the row number of the cluster center
c72198f1 143 return fRow;
d88f97cc 144}
145///////////////////////////////////////////////////////////////////////////
261c0caf 146Int_t AliCalcluster::GetColumn() const
d88f97cc 147{
148// Provide the column number of the cluster center
c72198f1 149 return fCol;
d88f97cc 150}
151///////////////////////////////////////////////////////////////////////////
261c0caf 152Int_t AliCalcluster::GetNmodules() const
d88f97cc 153{
154// Provide the number of modules in the cluster
155 return fNmods;
156}
157///////////////////////////////////////////////////////////////////////////
261c0caf 158Float_t AliCalcluster::GetRowDispersion() const
d88f97cc 159{
dafe31a2 160// Provide the normalised row dispersion of the cluster.
161 Float_t sig=GetSignal();
162 if (sig > 0.)
d88f97cc 163 {
dafe31a2 164 return fRowdisp/sig;
d88f97cc 165 }
166 else
167 {
168 return 0.;
169 }
170}
171///////////////////////////////////////////////////////////////////////////
261c0caf 172Float_t AliCalcluster::GetColumnDispersion() const
d88f97cc 173{
174// Provide the normalised column dispersion of the cluster
dafe31a2 175 Float_t sig=GetSignal();
176 if (sig > 0.)
d88f97cc 177 {
dafe31a2 178 return fColdisp/sig;
d88f97cc 179 }
180 else
181 {
182 return 0.;
183 }
184}
185///////////////////////////////////////////////////////////////////////////
186void AliCalcluster::Start(AliCalmodule& m)
187{
959fbac5 188// Reset cluster data and start with module m.
f40f8fbd 189// A module can only start a cluster when it contains a signal,
190// has not been used in a cluster yet, and is not declared dead.
191//
192// Note :
193// It is advised NOT to start a cluster with modules situated at a detector edge.
194// This feature is automatically checked when using the built-in clustering
195// of AliCalorimeter.
d88f97cc 196
dafe31a2 197 AliSignal::Reset();
198
f40f8fbd 199 Ali3Vector r;
d88f97cc 200
f40f8fbd 201 if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0)
d88f97cc 202 {
c72198f1 203 fRow=m.GetRow();
204 fCol=m.GetColumn();
f40f8fbd 205 r=m.GetPosition();
206 SetPosition(r);
1fbffa23 207 SetSignal(m.GetSignal(1,1)); // Use the gain etc... corrected module signal
d88f97cc 208 fNmods=1;
d88f97cc 209 fRowdisp=0.;
210 fColdisp=0.;
211 m.SetClusteredSignal(0.); // mark module as used in cluster
212 }
213 else
214 {
c72198f1 215 fRow=0;
216 fCol=0;
f40f8fbd 217 SetPosition(r);
d88f97cc 218 fNmods=0;
d88f97cc 219 fRowdisp=0.;
220 fColdisp=0.;
221 }
222}
223///////////////////////////////////////////////////////////////////////////
224void AliCalcluster::Add(AliCalmodule& m)
225{
f40f8fbd 226// Add module data to the cluster.
227// Dead modules are NOT added to the cluster.
dafe31a2 228// According to the distance of the module w.r.t. the cluster center
229// the various signal values are updated.
230
231 if (fNmods)
232 {
1fbffa23 233 if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0) // only add unused modules
dafe31a2 234 {
1fbffa23 235 Float_t sigm=m.GetSignal(1,1); // Use the gain etc... corrected module signal
236
dafe31a2 237 Int_t drow=int(fabs(double(GetRow()-m.GetRow()))); // row distance to center
238 Int_t dcol=int(fabs(double(GetColumn()-m.GetColumn()))); // column distance to center
239
240 // Determine the ring index for this module around the cluster center
241 Int_t jring=drow;
242 if (dcol>drow) jring=dcol;
243
244 Int_t nvalues=GetNvalues();
245
246 if ((jring+2)<=nvalues) // Module within existing ring(s) ==> Add module signal to the enclosing ring(s)
247 {
248 for (Int_t i=(jring+2); i<=nvalues; i++)
249 {
250 AddSignal(sigm,i);
251 }
252 }
253 else // Module outside all existing rings ==> Init. new ring signals with existing enclosed signal(s)
254 {
255 for (Int_t j=(nvalues+1); j<=(jring+2); j++)
256 {
257 SetSignal(GetSignal(j-1),j);
258 }
259 // Add current module signal to the signal value for the corresponding ring
260 AddSignal(sigm,(jring+2));
261 }
d88f97cc 262
dafe31a2 263 // Update total cluster signal
264 AddSignal(sigm);
d88f97cc 265
dafe31a2 266 fNmods+=1;
267 fRowdisp+=sigm*float(drow*drow);
268 fColdisp+=sigm*float(dcol*dcol);
269 m.SetClusteredSignal(0.); // mark module as used in cluster
270 }
271 }
272 else
d88f97cc 273 {
dafe31a2 274 cout << " *AliCalcluster::Add* No action. Cluster should be started first."
275 << endl;
d88f97cc 276 }
277}
278///////////////////////////////////////////////////////////////////////////
8e8e6c7f 279void AliCalcluster::AddVetoSignal(AliSignal& s,Int_t extr)
d88f97cc 280{
8e8e6c7f 281// Associate an (extrapolated) AliSignal as veto to the cluster.
282// By default a straight line extrapolation is performed which extrapolates
283// the signal position until the length of its position vector matches that
284// of the position vector of the cluster.
285// In this extrapolation procedure the error propagation is performed
286// automatically.
287// Based on the cluster and extrapolated veto signal (x,y) positions and
288// position errors the confidence level of association is calculated
289// and stored as an additional signal value.
290// By means of the GetVetoSignal memberfunction the confidence level of
291// association can always be updated by the user.
292// In case the user wants to invoke a more detailed extrapolation procedure,
293// the automatic extrapolation can be suppressed by setting the argument
294// extr=0. In this case it is assumed that the AliSignal as entered via
295// the argument contains already the extrapolated position vector and
296// corresponding errors.
297// Note : Three additional values are added to the original AliSignal
298// to hold the chi2, ndf and confidence level values of the association.
d88f97cc 299 if (!fVetos)
300 {
301 fNvetos=0;
302 fVetos=new TObjArray();
6516b62d 303 fVetos->SetOwner();
d88f97cc 304 }
305
8e8e6c7f 306 Int_t nvalues=s.GetNvalues();
c72198f1 307 AliSignal* sx=new AliSignal(s); // Additional values will be added
8e8e6c7f 308 TString name=s.GetName();
309 name.Append(" + additional chi2, ndf and CL values");
310 sx->SetName(name);
311
312 Double_t vecc[3],vecv[3];
c72198f1 313 if (extr)
8e8e6c7f 314 {
315 // Extrapolate the veto hit position
316 Double_t scale=1;
317 GetPosition(vecc,"sph");
318 s.GetPosition(vecv,"sph");
319 if (vecv[0]) scale=vecc[0]/vecv[0];
320 Ali3Vector r=s*scale;
321 sx->SetPosition(r);
322 }
d88f97cc 323
8e8e6c7f 324 // Calculate the confidence level of association
325 GetPosition(vecc,"car");
326 sx->GetPosition(vecv,"car");
327 Double_t dx=vecc[0]-vecv[0];
328 Double_t dy=vecc[1]-vecv[1];
329 GetPositionErrors(vecc,"car");
330 sx->GetPositionErrors(vecv,"car");
331 Double_t sxc2=vecc[0]*vecc[0];
332 Double_t syc2=vecc[1]*vecc[1];
333 Double_t sxv2=vecv[0]*vecv[0];
334 Double_t syv2=vecv[1]*vecv[1];
335 Double_t sumx2=sxc2+sxv2;
336 Double_t sumy2=syc2+syv2;
337 Double_t chi2=0;
338 if (sumx2>0 && sumy2>0) chi2=(dx*dx/sumx2)+(dy*dy/sumy2);
339 Int_t ndf=2;
340 AliMath m;
341 Double_t prob=m.Prob(chi2,ndf);
342 if (chi2>0) sx->SetSignal(chi2,nvalues+1);
343 if (ndf>0) sx->SetSignal(ndf,nvalues+2);
344 if (prob>0) sx->SetSignal(prob,nvalues+3);
345
346 fVetos->Add(sx);
347 fNvetos++;
d88f97cc 348}
349///////////////////////////////////////////////////////////////////////////
261c0caf 350Int_t AliCalcluster::GetNvetos() const
d88f97cc 351{
352// Provide the number of veto signals associated to the cluster
353 return fNvetos;
354}
355///////////////////////////////////////////////////////////////////////////
261c0caf 356AliSignal* AliCalcluster::GetVetoSignal(Int_t i) const
d88f97cc 357{
959fbac5 358// Provide access to the i-th veto signal of this cluster.
359// Note : The first hit corresponds to i=1.
6516b62d 360 if (!fVetos)
d88f97cc 361 {
6516b62d 362 cout << " *AliCalcluster::GetVetoSignal* No veto signals present." << endl;
363 return 0;
d88f97cc 364 }
365 else
366 {
6516b62d 367 if (i>0 && i<=fNvetos)
368 {
369 return (AliSignal*)fVetos->At(i-1);
370 }
371 else
372 {
373 cout << " *AliCalcluster::GetVetoSignal* Signal number " << i << " out of range."
374 << " Nvetos = " << fNvetos << endl;
375 return 0;
376 }
d88f97cc 377 }
378}
379///////////////////////////////////////////////////////////////////////////
261c0caf 380Float_t AliCalcluster::GetVetoLevel() const
959fbac5 381{
382// Provide the confidence level of best associated veto signal.
383 Float_t cl=0;
384 Float_t clmax=0;
8e8e6c7f 385 AliSignal* s=0;
386 Int_t nvalues=0;
959fbac5 387 if (fVetos)
388 {
389 for (Int_t i=0; i<fNvetos; i++)
390 {
8e8e6c7f 391 s=((AliSignal*)fVetos->At(i));
392 if (s)
393 {
394 nvalues=s->GetNvalues();
395 cl=s->GetSignal(nvalues);
396 if (cl>clmax) clmax=cl;
397 }
959fbac5 398 }
399 }
400 return clmax;
401}
402///////////////////////////////////////////////////////////////////////////
261c0caf 403Int_t AliCalcluster::HasVetoHit(Double_t cl) const
959fbac5 404{
405// Investigate if cluster has an associated veto hit with conf. level > cl.
406// Returns 1 if there is such an associated veto hit, otherwise returns 0.
407// Note : This function is faster than GetVetoLevel().
8e8e6c7f 408 AliSignal* s=0;
409 Int_t nvalues=0;
959fbac5 410 if (fVetos)
411 {
412 for (Int_t i=0; i<fNvetos; i++)
413 {
8e8e6c7f 414 s=((AliSignal*)fVetos->At(i));
415 if (s)
416 {
417 nvalues=s->GetNvalues();
418 if (s->GetSignal(nvalues) > cl) return 1;
419 }
959fbac5 420 }
421 }
422 return 0;
423}
424///////////////////////////////////////////////////////////////////////////