NcCalcluster.cxx

Go to the documentation of this file.

 
#include "NcCalcluster.h"
#include "Riostream.h"
 
ClassImp(NcCalcluster); // Class implementation to enable ROOT I/O
 
NcCalcluster::NcCalcluster() : NcSignal()
{
 
 fRow=0;
 fCol=0;
 fNmods=0;
 fRowdisp=0.;
 fColdisp=0.;
 fNvetos=0;
 fVetos=0;
 SetName("NcCalcluster [sig, sig11, sig33, sig55,...]");
}

NcCalcluster::~NcCalcluster()
{
 
 if (fVetos)
 {
  delete fVetos;
  fVetos=0;
 }
}

NcCalcluster::NcCalcluster(const NcCalcluster& c) : NcSignal(c)
{
 
 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++)
 {
  NcSignal* sx=c.GetVetoSignal(i);
  fVetos->Add(new NcSignal(*sx));
 }
}

NcCalcluster::NcCalcluster(NcCalmodule& m) : NcSignal()
{
 
 Nc3Vector r;
 
 Float_t sig=m.GetClusteredSignal();
 
 if (sig>0. && m.GetDeadValue()==0)
 {
  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;
  fRowdisp=0.;
  fColdisp=0.;
  m.SetClusteredSignal(0.); // mark module as used in cluster
  fNvetos=0;
  fVetos=0;
 }
 else
 {
  fRow=0;
  fCol=0;
  SetPosition(r);
  fNmods=0;
  fRowdisp=0.;
  fColdisp=0.;
  fNvetos=0;
  fVetos=0;
 }
 SetName("NcCalcluster [sig, sig11, sig33, sig55,...]");
}

Int_t NcCalcluster::GetRow() const
{
 
 return fRow;
}

Int_t NcCalcluster::GetColumn() const
{
 
 return fCol;
}

Int_t NcCalcluster::GetNmodules() const
{
 
 return fNmods;
}

Float_t NcCalcluster::GetRowDispersion() const
{
 
 Float_t sig=GetSignal();
 if (sig > 0.)
 {
  return fRowdisp/sig;
 }
 else
 {
  return 0.;
 }
}

Float_t NcCalcluster::GetColumnDispersion() const
{
 
 Float_t sig=GetSignal();
 if (sig > 0.)
 {
  return fColdisp/sig;
 }
 else
 {
  return 0.;
 }
}

void NcCalcluster::Start(NcCalmodule& m)
{
 
 NcSignal::Reset();
 
 Nc3Vector r;
 
 if (m.GetClusteredSignal()>0. && m.GetDeadValue()==0)
 {
  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;
  fRowdisp=0.;
  fColdisp=0.;
  m.SetClusteredSignal(0.); // mark module as used in cluster
 }
 else
 {
  fRow=0;
  fCol=0;
  SetPosition(r);
  fNmods=0;
  fRowdisp=0.;
  fColdisp=0.;
 }
}

void NcCalcluster::Add(NcCalmodule& m)
{
 
 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);
 
   fNmods+=1;
   fRowdisp+=sigm*float(drow*drow);
   fColdisp+=sigm*float(dcol*dcol);
   m.SetClusteredSignal(0.); // mark module as used in cluster
  }
 }
 else
 {
  cout << " *NcCalcluster::Add* No action. Cluster should be started first."
       << endl;
 }
}

void NcCalcluster::AddVetoSignal(NcSignal& s,Int_t extr)
{
 
 if (!fVetos)
 {
  fNvetos=0;
  fVetos=new TObjArray();
  fVetos->SetOwner();
 } 
 
 Int_t nvalues=s.GetNvalues();
 NcSignal* sx=new NcSignal(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];
  Nc3Vector 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;
 NcMath 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);
 
 fVetos->Add(sx);
 fNvetos++;
}

Int_t NcCalcluster::GetNvetos() const
{
 
 return fNvetos;
}

NcSignal* NcCalcluster::GetVetoSignal(Int_t i) const
{
 
 if (!fVetos)
 {
  cout << " *NcCalcluster::GetVetoSignal* No veto signals present." << endl;
  return 0;
 }
 else
 {
  if (i>0 && i<=fNvetos)
  {
   return (NcSignal*)fVetos->At(i-1);
  }
  else
  {
   cout << " *NcCalcluster::GetVetoSignal* Signal number " << i << " out of range."
        << " Nvetos = " << fNvetos << endl;
   return 0;
  }
 }
}

Float_t NcCalcluster::GetVetoLevel() const
{
 
 Float_t cl=0;
 Float_t clmax=0;
 NcSignal* s=0;
 Int_t nvalues=0;
 if (fVetos)
 {
  for (Int_t i=0; i<fNvetos; i++)
  {
   s=((NcSignal*)fVetos->At(i));
   if (s)
   {
    nvalues=s->GetNvalues();
    cl=s->GetSignal(nvalues);
    if (cl>clmax) clmax=cl;
   }
  }
 }
 return clmax;
}

Int_t NcCalcluster::HasVetoHit(Double_t cl) const
{
 
 NcSignal* s=0;
 Int_t nvalues=0;
 if (fVetos)
 {
  for (Int_t i=0; i<fNvetos; i++)
  {
   s=((NcSignal*)fVetos->At(i));
   if (s)
   {
    nvalues=s->GetNvalues();
    if (s->GetSignal(nvalues) > cl) return 1;
   }
  }
 }
 return 0;
}