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