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(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   SetSignal(sig);
 118   fNmods=1;
 119   fRowdisp=0.;
 120   fColdisp=0.;
 121   m.SetClusteredSignal(0.); // mark module as used in cluster
 122   fNvetos=0;
 123   fVetos=0;
 124  }
 125  else
 126  {
 127   fRow=0;
 128   fCol=0;
 129   SetPosition(r);
 130   fNmods=0;
 131   fRowdisp=0.;
 132   fColdisp=0.;
 133   fNvetos=0;
 134   fVetos=0;
 135  }
 136  SetName("AliCalcluster [sig, sig11, sig33, sig55,...]");
 137 }
 138 ///////////////////////////////////////////////////////////////////////////
 139 Int_t AliCalcluster::GetRow()
 140 {
 141 // Provide the row number of the cluster center
 142  return fRow;
 143 }
 144 ///////////////////////////////////////////////////////////////////////////
 145 Int_t AliCalcluster::GetColumn()
 146 {
 147 // Provide the column number of the cluster center
 148  return fCol;
 149 }
 150 ///////////////////////////////////////////////////////////////////////////
 151 Int_t AliCalcluster::GetNmodules()
 152 {
 153 // Provide the number of modules in the cluster
 154  return fNmods;
 155 }
 156 ///////////////////////////////////////////////////////////////////////////
 157 Float_t AliCalcluster::GetRowDispersion()
 158 {
 159 // Provide the normalised row dispersion of the cluster.
 160  Float_t sig=GetSignal();
 161  if (sig > 0.)
 162  {
 163   return fRowdisp/sig;
 164  }
 165  else
 166  {
 167   return 0.;
 168  }
 169 }
 170 ///////////////////////////////////////////////////////////////////////////
 171 Float_t AliCalcluster::GetColumnDispersion()
 172 {
 173 // Provide the normalised column dispersion of the cluster
 174  Float_t sig=GetSignal();
 175  if (sig > 0.)
 176  {
 177   return fColdisp/sig;
 178  }
 179  else
 180  {
 181   return 0.;
 182  }
 183 }
 184 ///////////////////////////////////////////////////////////////////////////
 185 void AliCalcluster::Start(AliCalmodule& m)
 186 {
 187 // Reset cluster data and start with module m.
 188 // A module can only start a cluster when it contains a signal,
 189 // has not been used in a cluster yet, and is not declared dead.
 190 //
 191 // Note :
 192 // It is advised NOT to start a cluster with modules situated at a detector edge.
 193 // This feature is automatically checked when using the built-in clustering
 194 // of AliCalorimeter.
 195
 196  AliSignal::Reset();
 197
 198  Ali3Vector r;
 199
 200  if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0)
 201  {
 202   fRow=m.GetRow();
 203   fCol=m.GetColumn();
 204   r=m.GetPosition();
 205   SetPosition(r);
 206   SetSignal(m.GetSignal());
 207   fNmods=1;
 208   fRowdisp=0.;
 209   fColdisp=0.;
 210   m.SetClusteredSignal(0.); // mark module as used in cluster
 211  }
 212  else
 213  {
 214   fRow=0;
 215   fCol=0;
 216   SetPosition(r);
 217   fNmods=0;
 218   fRowdisp=0.;
 219   fColdisp=0.;
 220  }
 221 }
 222 ///////////////////////////////////////////////////////////////////////////
 223 void AliCalcluster::Add(AliCalmodule& m)
 224 {
 225 // Add module data to the cluster.
 226 // Dead modules are NOT added to the cluster.
 227 // According to the distance of the module w.r.t. the cluster center
 228 // the various signal values are updated.
 229
 230  if (fNmods)
 231  {
 232   Float_t sigm=m.GetClusteredSignal();
 233
 234   if (sigm>0. && m.GetDeadValue()==0) // only add unused modules
 235   {
 236    Int_t drow=int(fabs(double(GetRow()-m.GetRow())));       // row distance to center
 237    Int_t dcol=int(fabs(double(GetColumn()-m.GetColumn()))); // column distance to center
 238
 239    // Determine the ring index for this module around the cluster center
 240    Int_t jring=drow;
 241    if (dcol>drow) jring=dcol;
 242
 243    Int_t nvalues=GetNvalues();
 244
 245    if ((jring+2)<=nvalues) // Module within existing ring(s) ==> Add module signal to the enclosing ring(s)
 246    {
 247     for (Int_t i=(jring+2); i<=nvalues; i++)
 248     {
 249      AddSignal(sigm,i);
 250     }
 251    }
 252    else  // Module outside all existing rings ==> Init. new ring signals with existing enclosed signal(s)
 253    {
 254     for (Int_t j=(nvalues+1); j<=(jring+2); j++)
 255     {
 256      SetSignal(GetSignal(j-1),j);
 257     }
 258     // Add current module signal to the signal value for the corresponding ring
 259     AddSignal(sigm,(jring+2));
 260    }
 261
 262    // Update total cluster signal
 263    AddSignal(sigm);
 264
 265    fNmods+=1;
 266    fRowdisp+=sigm*float(drow*drow);
 267    fColdisp+=sigm*float(dcol*dcol);
 268    m.SetClusteredSignal(0.); // mark module as used in cluster
 269   }
 270  }
 271  else
 272  {
 273   cout << " *AliCalcluster::Add* No action. Cluster should be started first."
 274        << endl;
 275  }
 276 }
 277 ///////////////////////////////////////////////////////////////////////////
 278 void AliCalcluster::AddVetoSignal(AliSignal& s,Int_t extr)
 279 {
 280 // Associate an (extrapolated) AliSignal as veto to the cluster.
 281 // By default a straight line extrapolation is performed which extrapolates
 282 // the signal position until the length of its position vector matches that
 283 // of the position vector of the cluster.
 284 // In this extrapolation procedure the error propagation is performed
 285 // automatically.
 286 // Based on the cluster and extrapolated veto signal (x,y) positions and
 287 // position errors the confidence level of association is calculated
 288 // and stored as an additional signal value.
 289 // By means of the GetVetoSignal memberfunction the confidence level of
 290 // association can always be updated by the user.
 291 // In case the user wants to invoke a more detailed extrapolation procedure,
 292 // the automatic extrapolation can be suppressed by setting the argument
 293 // extr=0. In this case it is assumed that the AliSignal as entered via
 294 // the argument contains already the extrapolated position vector and
 295 // corresponding errors.
 296 // Note : Three additional values are added to the original AliSignal
 297 //        to hold the chi2, ndf and confidence level values of the association.
 298  if (!fVetos)
 299  {
 300   fNvetos=0;
 301   fVetos=new TObjArray();
 302   fVetos->SetOwner();
 303  }
 304
 305  Int_t nvalues=s.GetNvalues();
 306  AliSignal* sx=new AliSignal(s); // Additional values will be added
 307  TString name=s.GetName();
 308  name.Append(" + additional chi2, ndf and CL values");
 309  sx->SetName(name);
 310
 311  Double_t vecc[3],vecv[3];
 312  if (extr)
 313  {
 314   // Extrapolate the veto hit position
 315   Double_t scale=1;
 316   GetPosition(vecc,"sph");
 317   s.GetPosition(vecv,"sph");
 318   if (vecv[0]) scale=vecc[0]/vecv[0];
 319   Ali3Vector r=s*scale;
 320   sx->SetPosition(r);
 321  }
 322
 323  // Calculate the confidence level of association
 324  GetPosition(vecc,"car");
 325  sx->GetPosition(vecv,"car");
 326  Double_t dx=vecc[0]-vecv[0];
 327  Double_t dy=vecc[1]-vecv[1];
 328  GetPositionErrors(vecc,"car");
 329  sx->GetPositionErrors(vecv,"car");
 330  Double_t sxc2=vecc[0]*vecc[0];
 331  Double_t syc2=vecc[1]*vecc[1];
 332  Double_t sxv2=vecv[0]*vecv[0];
 333  Double_t syv2=vecv[1]*vecv[1];
 334  Double_t sumx2=sxc2+sxv2;
 335  Double_t sumy2=syc2+syv2;
 336  Double_t chi2=0;
 337  if (sumx2>0 && sumy2>0) chi2=(dx*dx/sumx2)+(dy*dy/sumy2);
 338  Int_t ndf=2;
 339  AliMath m;
 340  Double_t prob=m.Prob(chi2,ndf);
 341  if (chi2>0) sx->SetSignal(chi2,nvalues+1);
 342  if (ndf>0) sx->SetSignal(ndf,nvalues+2);
 343  if (prob>0) sx->SetSignal(prob,nvalues+3);
 344
 345  fVetos->Add(sx);
 346  fNvetos++;
 347 }
 348 ///////////////////////////////////////////////////////////////////////////
 349 Int_t AliCalcluster::GetNvetos()
 350 {
 351 // Provide the number of veto signals associated to the cluster
 352  return fNvetos;
 353 }
 354 ///////////////////////////////////////////////////////////////////////////
 355 AliSignal* AliCalcluster::GetVetoSignal(Int_t i)
 356 {
 357 // Provide access to the i-th veto signal of this cluster.
 358 // Note : The first hit corresponds to i=1.
 359  if (!fVetos)
 360  {
 361   cout << " *AliCalcluster::GetVetoSignal* No veto signals present." << endl;
 362   return 0;
 363  }
 364  else
 365  {
 366   if (i>0 && i<=fNvetos)
 367   {
 368    return (AliSignal*)fVetos->At(i-1);
 369   }
 370   else
 371   {
 372    cout << " *AliCalcluster::GetVetoSignal* Signal number " << i << " out of range."
 373         << " Nvetos = " << fNvetos << endl;
 374    return 0;
 375   }
 376  }
 377 }
 378 ///////////////////////////////////////////////////////////////////////////
 379 Float_t AliCalcluster::GetVetoLevel()
 380 {
 381 // Provide the confidence level of best associated veto signal.
 382  Float_t cl=0;
 383  Float_t clmax=0;
 384  AliSignal* s=0;
 385  Int_t nvalues=0;
 386  if (fVetos)
 387  {
 388   for (Int_t i=0; i<fNvetos; i++)
 389   {
 390    s=((AliSignal*)fVetos->At(i));
 391    if (s)
 392    {
 393     nvalues=s->GetNvalues();
 394     cl=s->GetSignal(nvalues);
 395     if (cl>clmax) clmax=cl;
 396    }
 397   }
 398  }
 399  return clmax;
 400 }
 401 ///////////////////////////////////////////////////////////////////////////
 402 Int_t AliCalcluster::HasVetoHit(Double_t cl)
 403 {
 404 // Investigate if cluster has an associated veto hit with conf. level > cl.
 405 // Returns 1 if there is such an associated veto hit, otherwise returns 0.
 406 // Note : This function is faster than GetVetoLevel().
 407  AliSignal* s=0;
 408  Int_t nvalues=0;
 409  if (fVetos)
 410  {
 411   for (Int_t i=0; i<fNvetos; i++)
 412   {
 413    s=((AliSignal*)fVetos->At(i));
 414    if (s)
 415    {
 416     nvalues=s->GetNvalues();
 417     if (s->GetSignal(nvalues) > cl) return 1;
 418    }
 419   }
 420  }
 421  return 0;
 422 }
 423 ///////////////////////////////////////////////////////////////////////////