Deleted erroneously added simlinks
[u/mrichter/AliRoot.git] / RALICE / AliCalcluster.cxx
index afe3770..393df68 100644 (file)
  * provided "as is" without express or implied warranty.                  *
  **************************************************************************/
 
-/*
-$Log$
-Revision 1.2  1999/09/29 09:24:28  fca
-Introduction of the Copyright and cvs Log
-
-*/
+// $Id$
 
 ///////////////////////////////////////////////////////////////////////////
 // Class AliCalcluster
 // Description of a cluster of calorimeter modules.
-// A matrix geometry is assumed in which a cluster center
-// is identified by (row,col) and contains sig as signal
-// being the signal of the complete cluster.
-// Some info about cluster topology is provided in order
+// A 2D (matrix) geometry is assumed in which a cluster center is identified
+// by two integer indices (i,j), e.g. row and column indicators.
+//
+// The 1st signal value is the signal of the complete cluster.
+// This is the signal which is provided as default by invoking GetSignal().
+//
+// In case clustering/grouping of module signals was performed over several
+// rings around the center (see e.g. AliCalorimeter::Group), the following
+// additional information is provided by the various signal values :
+//
+// The 2nd signal value is the original signal of the central module.
+// The 3rd signal value is the total signal within the 1st (i.e. 3x3) ring of
+// modules around the cluster center.
+// The 4th signal value is the total signal within the 2nd (i.e. 5x5) ring of
+// modules around the cluster center.
+// Etc....
+//
+// Note : In case the cluster consists of only 1 module, then only the
+//        1st signal value will be present (for obvious reasons).
+//
+// Some dispersion info about cluster topology is provided in order
 // to enable EM or hadronic cluster identification.
 //
 //--- Author: Nick van Eijndhoven 13-jun-1997 UU-SAP Utrecht
-//- Modified: NvE 31-oct-1999 UU-SAP Utrecht 
+//- Modified: NvE $Date$ UU-SAP Utrecht
 ///////////////////////////////////////////////////////////////////////////
 
 #include "AliCalcluster.h"
+#include "Riostream.h"
  
 ClassImp(AliCalcluster) // Class implementation to enable ROOT I/O
  
-AliCalcluster::AliCalcluster()
+AliCalcluster::AliCalcluster() : AliSignal()
 {
-// Default constructer, all data is set to 0
- fCenter=0;
- fSig=0.;
+// Default constructor, all data is set to 0
+ fRow=0;
+ fCol=0;
  fNmods=0;
- fSig11=0.;
- fSig33=0.;
- fSig55=0.;
  fRowdisp=0.;
  fColdisp=0.;
  fNvetos=0;
  fVetos=0;
+ SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
 }
 ///////////////////////////////////////////////////////////////////////////
 AliCalcluster::~AliCalcluster()
@@ -57,31 +68,55 @@ AliCalcluster::~AliCalcluster()
 // Destructor to delete dynamically allocated memory
  if (fVetos)
  {
-  fVetos->Delete();
   delete fVetos;
   fVetos=0;
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-AliCalcluster::AliCalcluster(AliCalmodule& m)
+AliCalcluster::AliCalcluster(const AliCalcluster& c) : AliSignal(c)
+{
+// Copy constructor
+ fRow=c.fRow;
+ fCol=c.fCol;
+ fNmods=c.fNmods;
+ fRowdisp=c.fRowdisp;
+ fColdisp=c.fColdisp;
+ fNvetos=c.fNvetos;
+
+ fVetos=new TObjArray();
+ fVetos->SetOwner();
+
+ for (Int_t i=1; i<=fNvetos; i++)
+ {
+  AliSignal* sx=c.GetVetoSignal(i);
+  fVetos->Add(new AliSignal(*sx));
+ }
+}
+///////////////////////////////////////////////////////////////////////////
+AliCalcluster::AliCalcluster(AliCalmodule& m) : AliSignal()
 {
 // Cluster constructor with module m as center.
-// Module data is only entered for a module which contains
-// a signal and has not been used in a cluster yet.
-// Modules situated at a detector edge are not allowed to start a cluster.
+// Module data is only entered for a module which contains a signal,
+// has not been used in a cluster yet, and is not declared dead.
+//
+// Note :
+// It is advised NOT to start a cluster with modules situated at a detector edge.
+// This feature is automatically checked when using the built-in clustering
+// of AliCalorimeter.  
+
+ Ali3Vector r;
 
- Float_t pos[3]={0,0,0};
+ Float_t sig=m.GetClusteredSignal();
 
- if ((m.GetClusteredSignal() > 0.) && (m.GetEdgeValue() == 0))
+ if (sig>0. && m.GetDeadValue()==0)
  {
-  fCenter=&m;
-  m.GetPosition(pos,"sph");
-  SetPosition(pos,"sph");
-  fSig=m.GetClusteredSignal();
+  fRow=m.GetRow();
+  fCol=m.GetColumn();
+  r=m.GetPosition();
+  SetPosition(r);
+  sig=m.GetSignal(1,1); // Use the gain etc... corrected module signal
+  SetSignal(sig);
   fNmods=1;
-  fSig11=m.GetClusteredSignal();
-  fSig33=m.GetClusteredSignal();
-  fSig55=m.GetClusteredSignal();
   fRowdisp=0.;
   fColdisp=0.;
   m.SetClusteredSignal(0.); // mark module as used in cluster
@@ -90,84 +125,43 @@ AliCalcluster::AliCalcluster(AliCalmodule& m)
  }
  else
  {
-  fCenter=0;
-  SetPosition(pos,"sph");
-  fSig=0.;
+  fRow=0;
+  fCol=0;
+  SetPosition(r);
   fNmods=0;
-  fSig11=0.;
-  fSig33=0.;
-  fSig55=0.;
   fRowdisp=0.;
   fColdisp=0.;
   fNvetos=0;
   fVetos=0;
-}
+ }
+ SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCalcluster::GetRow()
+Int_t AliCalcluster::GetRow() const
 {
 // Provide the row number of the cluster center
- if (fCenter)
- {
-  return fCenter->GetRow();
- }
- else
- {
-  return 0;
- }
+ return fRow;
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCalcluster::GetColumn()
+Int_t AliCalcluster::GetColumn() const
 {
 // Provide the column number of the cluster center
- if (fCenter)
- {
-  return fCenter->GetColumn();
- }
- else
- {
-  return 0;
- }
+ return fCol;
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliCalcluster::GetSignal(Int_t n)
-{
-// Provide the total signal of a module matrix of n modules around
-// the cluster center.
-// Example : n=9 --> total signal in the 3x3 matrix
-//             1 --> signal of central module
-// Note : n=0 provides the total cluster signal (Default)
- Float_t signal=-1;
- if (n == 0)  signal=fSig;
- if (n == 1)  signal=fSig11;
- if (n == 9)  signal=fSig33;
- if (n == 25) signal=fSig55;
- if (signal > 0.)
- {
-  return signal;
- }
- else
- {
-  cout << " *AliCalcluster::GetSignal* Invalid argument n = " << n << endl;
-  return 0;
- }
-}
-///////////////////////////////////////////////////////////////////////////
-Int_t AliCalcluster::GetNmodules()
+Int_t AliCalcluster::GetNmodules() const
 {
 // Provide the number of modules in the cluster
  return fNmods;
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliCalcluster::GetRowDispersion()
+Float_t AliCalcluster::GetRowDispersion() const
 {
-// Provide the normalised row dispersion of the cluster
- if (fSig > 0.)
+// Provide the normalised row dispersion of the cluster.
+ Float_t sig=GetSignal();
+ if (sig > 0.)
  {
-  return fRowdisp/fSig;
+  return fRowdisp/sig;
  }
  else
  {
@@ -175,12 +169,13 @@ Float_t AliCalcluster::GetRowDispersion()
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliCalcluster::GetColumnDispersion()
+Float_t AliCalcluster::GetColumnDispersion() const
 {
 // Provide the normalised column dispersion of the cluster
- if (fSig > 0.)
+ Float_t sig=GetSignal();
+ if (sig > 0.)
  {
-  return fColdisp/fSig;
+  return fColdisp/sig;
  }
  else
  {
@@ -191,35 +186,36 @@ Float_t AliCalcluster::GetColumnDispersion()
 void AliCalcluster::Start(AliCalmodule& m)
 {
 // Reset cluster data and start with module m.
-// A module can only start a cluster when it contains
-// a signal, has not been used in a cluster yet, is not
-// situated at a detector edge and is not declared dead.
+// A module can only start a cluster when it contains a signal,
+// has not been used in a cluster yet, and is not declared dead.
+//
+// Note :
+// It is advised NOT to start a cluster with modules situated at a detector edge.
+// This feature is automatically checked when using the built-in clustering
+// of AliCalorimeter.  
+
+ AliSignal::Reset();
 
- Float_t pos[3]={0,0,0};
+ Ali3Vector r;
 
- if (m.GetClusteredSignal()>0. && m.GetEdgeValue()==0 && m.GetDeadValue()==0)
+ if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0)
  {
-  fCenter=&m;
-  m.GetPosition(pos,"sph");
-  SetPosition(pos,"sph");
-  fSig=m.GetSignal();
+  fRow=m.GetRow();
+  fCol=m.GetColumn();
+  r=m.GetPosition();
+  SetPosition(r);
+  SetSignal(m.GetSignal(1,1)); // Use the gain etc... corrected module signal
   fNmods=1;
-  fSig11=m.GetSignal();
-  fSig33=m.GetSignal();
-  fSig55=m.GetSignal();
   fRowdisp=0.;
   fColdisp=0.;
   m.SetClusteredSignal(0.); // mark module as used in cluster
  }
  else
  {
-  fCenter=0;
-  SetPosition(pos,"sph");
-  fSig=0.;
+  fRow=0;
+  fCol=0;
+  SetPosition(r);
   fNmods=0;
-  fSig11=0.;
-  fSig33=0.;
-  fSig55=0.;
   fRowdisp=0.;
   fColdisp=0.;
  }
@@ -227,93 +223,200 @@ void AliCalcluster::Start(AliCalmodule& m)
 ///////////////////////////////////////////////////////////////////////////
 void AliCalcluster::Add(AliCalmodule& m)
 {
-// Add module data to the cluster
+// Add module data to the cluster.
+// Dead modules are NOT added to the cluster.
+// According to the distance of the module w.r.t. the cluster center
+// the various signal values are updated.
 
- Float_t signal=m.GetClusteredSignal();
+ if (fNmods)
+ {
+  if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0) // only add unused modules
+  {
+   Float_t sigm=m.GetSignal(1,1); // Use the gain etc... corrected module signal
+
+   Int_t drow=int(fabs(double(GetRow()-m.GetRow())));       // row distance to center
+   Int_t dcol=int(fabs(double(GetColumn()-m.GetColumn()))); // column distance to center 
+
+   // Determine the ring index for this module around the cluster center
+   Int_t jring=drow;
+   if (dcol>drow) jring=dcol;
+
+   Int_t nvalues=GetNvalues();
+
+   if ((jring+2)<=nvalues) // Module within existing ring(s) ==> Add module signal to the enclosing ring(s)
+   {
+    for (Int_t i=(jring+2); i<=nvalues; i++)
+    {
+     AddSignal(sigm,i);
+    }
+   }
+   else  // Module outside all existing rings ==> Init. new ring signals with existing enclosed signal(s)
+   {
+    for (Int_t j=(nvalues+1); j<=(jring+2); j++)
+    {
+     SetSignal(GetSignal(j-1),j);
+    }
+    // Add current module signal to the signal value for the corresponding ring
+    AddSignal(sigm,(jring+2));
+   }
+
+   // Update total cluster signal
+   AddSignal(sigm);
  
- if (signal > 0.) // only add unused modules
+   fNmods+=1;
+   fRowdisp+=sigm*float(drow*drow);
+   fColdisp+=sigm*float(dcol*dcol);
+   m.SetClusteredSignal(0.); // mark module as used in cluster
+  }
+ }
+ else
  {
-  fSig+=signal;
-  fNmods+=1;
-  Int_t drow=int(fabs(double(GetRow()-m.GetRow())));       // row distance to center
-  Int_t dcol=int(fabs(double(GetColumn()-m.GetColumn()))); // column distance to center
-  if ((drow < 2) && (dcol < 2)) fSig33+=signal;
-  if ((drow < 3) && (dcol < 3)) fSig55+=signal;
-  fRowdisp+=signal*float(drow*drow);
-  fColdisp+=signal*float(dcol*dcol);
-  m.SetClusteredSignal(0.); // mark module as used in cluster
+  cout << " *AliCalcluster::Add* No action. Cluster should be started first."
+       << endl;
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-void AliCalcluster::AddVetoSignal(Float_t* r,TString f,Float_t s)
+void AliCalcluster::AddVetoSignal(AliSignal& s,Int_t extr)
 {
-// Associate an (extrapolated) AliSignal at location r as veto to the cluster.
-// The signal value s indicates the confidence level of hit association
-// and has to be provided by the user.
-// Note : The default signal value (s) is 0
+// Associate an (extrapolated) AliSignal as veto to the cluster.
+// By default a straight line extrapolation is performed which extrapolates
+// the signal position until the length of its position vector matches that
+// of the position vector of the cluster.
+// In this extrapolation procedure the error propagation is performed 
+// automatically.  
+// Based on the cluster and extrapolated veto signal (x,y) positions and
+// position errors the confidence level of association is calculated
+// and stored as an additional signal value.
+// By means of the GetVetoSignal memberfunction the confidence level of
+// association can always be updated by the user.
+// In case the user wants to invoke a more detailed extrapolation procedure,
+// the automatic extrapolation can be suppressed by setting the argument
+// extr=0. In this case it is assumed that the AliSignal as entered via
+// the argument contains already the extrapolated position vector and
+// corresponding errors. 
+// Note : Three additional values are added to the original AliSignal
+//        to hold the chi2, ndf and confidence level values of the association. 
  if (!fVetos)
  {
   fNvetos=0;
   fVetos=new TObjArray();
+  fVetos->SetOwner();
  } 
 
- fVetos->Add(new AliSignal(1));
- fNvetos++;
+ Int_t nvalues=s.GetNvalues();
+ AliSignal* sx=new AliSignal(s); // Additional values will be added
+ TString name=s.GetName();
+ name.Append(" + additional chi2, ndf and CL values");
+ sx->SetName(name);
+
+ Double_t vecc[3],vecv[3];
+ if (extr)
+ {
+  // Extrapolate the veto hit position
+  Double_t scale=1;
+  GetPosition(vecc,"sph");
+  s.GetPosition(vecv,"sph"); 
+  if (vecv[0]) scale=vecc[0]/vecv[0];
+  Ali3Vector r=s*scale;
+  sx->SetPosition(r);
+ }
+
+ // Calculate the confidence level of association
+ GetPosition(vecc,"car");
+ sx->GetPosition(vecv,"car"); 
+ Double_t dx=vecc[0]-vecv[0];
+ Double_t dy=vecc[1]-vecv[1];
+ GetPositionErrors(vecc,"car");
+ sx->GetPositionErrors(vecv,"car"); 
+ Double_t sxc2=vecc[0]*vecc[0];
+ Double_t syc2=vecc[1]*vecc[1];
+ Double_t sxv2=vecv[0]*vecv[0];
+ Double_t syv2=vecv[1]*vecv[1];
+ Double_t sumx2=sxc2+sxv2;
+ Double_t sumy2=syc2+syv2;
+ Double_t chi2=0;
+ if (sumx2>0 && sumy2>0) chi2=(dx*dx/sumx2)+(dy*dy/sumy2);
+ Int_t ndf=2;
+ AliMath m;
+ Double_t prob=m.Prob(chi2,ndf);
+ if (chi2>0) sx->SetSignal(chi2,nvalues+1);
+ if (ndf>0) sx->SetSignal(ndf,nvalues+2);
+ if (prob>0) sx->SetSignal(prob,nvalues+3);
 
- ((AliSignal*)fVetos->At(fNvetos-1))->SetPosition(r,f);
- ((AliSignal*)fVetos->At(fNvetos-1))->SetSignal(s);
+ fVetos->Add(sx);
+ fNvetos++;
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCalcluster::GetNvetos()
+Int_t AliCalcluster::GetNvetos() const
 {
 // Provide the number of veto signals associated to the cluster
  return fNvetos;
 }
 ///////////////////////////////////////////////////////////////////////////
-AliSignal* AliCalcluster::GetVetoSignal(Int_t i)
+AliSignal* AliCalcluster::GetVetoSignal(Int_t i) const
 {
 // Provide access to the i-th veto signal of this cluster.
 // Note : The first hit corresponds to i=1.
-
- if (i>0 && i<=fNvetos)
+ if (!fVetos)
  {
-  return (AliSignal*)fVetos->At(i-1);
+  cout << " *AliCalcluster::GetVetoSignal* No veto signals present." << endl;
+  return 0;
  }
  else
  {
-  cout << " *AliCalcluster::GetVetoSignal* Signal number " << i
-       << " out of range." << endl;
-  cout << " --- First signal (if any) returned." << endl;
-  return (AliSignal*)fVetos->At(0);
+  if (i>0 && i<=fNvetos)
+  {
+   return (AliSignal*)fVetos->At(i-1);
+  }
+  else
+  {
+   cout << " *AliCalcluster::GetVetoSignal* Signal number " << i << " out of range."
+        << " Nvetos = " << fNvetos << endl;
+   return 0;
+  }
  }
 }
 ///////////////////////////////////////////////////////////////////////////
-Float_t AliCalcluster::GetVetoLevel()
+Float_t AliCalcluster::GetVetoLevel() const
 {
 // Provide the confidence level of best associated veto signal.
  Float_t cl=0;
  Float_t clmax=0;
+ AliSignal* s=0;
+ Int_t nvalues=0;
  if (fVetos)
  {
   for (Int_t i=0; i<fNvetos; i++)
   {
-   cl=((AliSignal*)fVetos->At(i))->GetSignal();
-   if (cl>clmax) clmax=cl;
+   s=((AliSignal*)fVetos->At(i));
+   if (s)
+   {
+    nvalues=s->GetNvalues();
+    cl=s->GetSignal(nvalues);
+    if (cl>clmax) clmax=cl;
+   }
   }
  }
  return clmax;
 }
 ///////////////////////////////////////////////////////////////////////////
-Int_t AliCalcluster::HasVetoHit(Double_t cl)
+Int_t AliCalcluster::HasVetoHit(Double_t cl) const
 {
 // Investigate if cluster has an associated veto hit with conf. level > cl.
 // Returns 1 if there is such an associated veto hit, otherwise returns 0.
 // Note : This function is faster than GetVetoLevel().
+ AliSignal* s=0;
+ Int_t nvalues=0;
  if (fVetos)
  {
   for (Int_t i=0; i<fNvetos; i++)
   {
-   if (((AliSignal*)fVetos->At(i))->GetSignal() > cl) return 1;
+   s=((AliSignal*)fVetos->At(i));
+   if (s)
+   {
+    nvalues=s->GetNvalues();
+    if (s->GetSignal(nvalues) > cl) return 1;
+   }
   }
  }
  return 0;