04-apr-2004 NvE SetMass() invoked from AliTrack::Set3Momentum to get also the energy...
[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 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
42 // to enable EM or hadronic cluster identification.
43 //
44 //--- Author: Nick van Eijndhoven 13-jun-1997 UU-SAP Utrecht
45 //- Modified: NvE $Date$ UU-SAP Utrecht
46 ///////////////////////////////////////////////////////////////////////////
47
48 #include "AliCalcluster.h"
49 #include "Riostream.h"
50  
51 ClassImp(AliCalcluster) // Class implementation to enable ROOT I/O
52  
53 AliCalcluster::AliCalcluster() : AliSignal()
54 {
55 // Default constructor, all data is set to 0
56  fRow=0;
57  fCol=0;
58  fNmods=0;
59  fRowdisp=0.;
60  fColdisp=0.;
61  fNvetos=0;
62  fVetos=0;
63  SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
64 }
65 ///////////////////////////////////////////////////////////////////////////
66 AliCalcluster::~AliCalcluster()
67 {
68 // Destructor to delete dynamically allocated memory
69  if (fVetos)
70  {
71   delete fVetos;
72   fVetos=0;
73  }
74 }
75 ///////////////////////////////////////////////////////////////////////////
76 AliCalcluster::AliCalcluster(const AliCalcluster& c) : AliSignal(c)
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 ///////////////////////////////////////////////////////////////////////////
96 AliCalcluster::AliCalcluster(AliCalmodule& m) : AliSignal()
97 {
98 // Cluster constructor with module m as center.
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.  
106
107  Ali3Vector r;
108
109  Float_t sig=m.GetClusteredSignal();
110
111  if (sig>0. && m.GetDeadValue()==0)
112  {
113   fRow=m.GetRow();
114   fCol=m.GetColumn();
115   r=m.GetPosition();
116   SetPosition(r);
117   sig=m.GetSignal(1,1); // Use the gain etc... corrected module signal
118   SetSignal(sig);
119   fNmods=1;
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  {
128   fRow=0;
129   fCol=0;
130   SetPosition(r);
131   fNmods=0;
132   fRowdisp=0.;
133   fColdisp=0.;
134   fNvetos=0;
135   fVetos=0;
136  }
137  SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
138 }
139 ///////////////////////////////////////////////////////////////////////////
140 Int_t AliCalcluster::GetRow() const
141 {
142 // Provide the row number of the cluster center
143  return fRow;
144 }
145 ///////////////////////////////////////////////////////////////////////////
146 Int_t AliCalcluster::GetColumn() const
147 {
148 // Provide the column number of the cluster center
149  return fCol;
150 }
151 ///////////////////////////////////////////////////////////////////////////
152 Int_t AliCalcluster::GetNmodules() const
153 {
154 // Provide the number of modules in the cluster
155  return fNmods;
156 }
157 ///////////////////////////////////////////////////////////////////////////
158 Float_t AliCalcluster::GetRowDispersion() const
159 {
160 // Provide the normalised row dispersion of the cluster.
161  Float_t sig=GetSignal();
162  if (sig > 0.)
163  {
164   return fRowdisp/sig;
165  }
166  else
167  {
168   return 0.;
169  }
170 }
171 ///////////////////////////////////////////////////////////////////////////
172 Float_t AliCalcluster::GetColumnDispersion() const
173 {
174 // Provide the normalised column dispersion of the cluster
175  Float_t sig=GetSignal();
176  if (sig > 0.)
177  {
178   return fColdisp/sig;
179  }
180  else
181  {
182   return 0.;
183  }
184 }
185 ///////////////////////////////////////////////////////////////////////////
186 void AliCalcluster::Start(AliCalmodule& m)
187 {
188 // Reset cluster data and start with module m.
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.  
196
197  AliSignal::Reset();
198
199  Ali3Vector r;
200
201  if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0)
202  {
203   fRow=m.GetRow();
204   fCol=m.GetColumn();
205   r=m.GetPosition();
206   SetPosition(r);
207   SetSignal(m.GetSignal(1,1)); // Use the gain etc... corrected module signal
208   fNmods=1;
209   fRowdisp=0.;
210   fColdisp=0.;
211   m.SetClusteredSignal(0.); // mark module as used in cluster
212  }
213  else
214  {
215   fRow=0;
216   fCol=0;
217   SetPosition(r);
218   fNmods=0;
219   fRowdisp=0.;
220   fColdisp=0.;
221  }
222 }
223 ///////////////////////////////////////////////////////////////////////////
224 void AliCalcluster::Add(AliCalmodule& m)
225 {
226 // Add module data to the cluster.
227 // Dead modules are NOT added to the cluster.
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  {
233   if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0) // only add unused modules
234   {
235    Float_t sigm=m.GetSignal(1,1); // Use the gain etc... corrected module signal
236
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    }
262
263    // Update total cluster signal
264    AddSignal(sigm);
265  
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
273  {
274   cout << " *AliCalcluster::Add* No action. Cluster should be started first."
275        << endl;
276  }
277 }
278 ///////////////////////////////////////////////////////////////////////////
279 void AliCalcluster::AddVetoSignal(AliSignal& s,Int_t extr)
280 {
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. 
299  if (!fVetos)
300  {
301   fNvetos=0;
302   fVetos=new TObjArray();
303   fVetos->SetOwner();
304  } 
305
306  Int_t nvalues=s.GetNvalues();
307  AliSignal* sx=new AliSignal(s); // Additional values will be added
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];
313  if (extr)
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  }
323
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++;
348 }
349 ///////////////////////////////////////////////////////////////////////////
350 Int_t AliCalcluster::GetNvetos() const
351 {
352 // Provide the number of veto signals associated to the cluster
353  return fNvetos;
354 }
355 ///////////////////////////////////////////////////////////////////////////
356 AliSignal* AliCalcluster::GetVetoSignal(Int_t i) const
357 {
358 // Provide access to the i-th veto signal of this cluster.
359 // Note : The first hit corresponds to i=1.
360  if (!fVetos)
361  {
362   cout << " *AliCalcluster::GetVetoSignal* No veto signals present." << endl;
363   return 0;
364  }
365  else
366  {
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   }
377  }
378 }
379 ///////////////////////////////////////////////////////////////////////////
380 Float_t AliCalcluster::GetVetoLevel() const
381 {
382 // Provide the confidence level of best associated veto signal.
383  Float_t cl=0;
384  Float_t clmax=0;
385  AliSignal* s=0;
386  Int_t nvalues=0;
387  if (fVetos)
388  {
389   for (Int_t i=0; i<fNvetos; i++)
390   {
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    }
398   }
399  }
400  return clmax;
401 }
402 ///////////////////////////////////////////////////////////////////////////
403 Int_t AliCalcluster::HasVetoHit(Double_t cl) const
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().
408  AliSignal* s=0;
409  Int_t nvalues=0;
410  if (fVetos)
411  {
412   for (Int_t i=0; i<fNvetos; i++)
413   {
414    s=((AliSignal*)fVetos->At(i));
415    if (s)
416    {
417     nvalues=s->GetNvalues();
418     if (s->GetSignal(nvalues) > cl) return 1;
419    }
420   }
421  }
422  return 0;
423 }
424 ///////////////////////////////////////////////////////////////////////////