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 ///////////////////////////////////////////////////////////////////////////