VtxCluster.cxx

Go to the documentation of this file.

//                     VtxCluster.h  -  description
//                          -------------------
//    begin                : Thu June 3 2004
//    copyright            : (C) 2004 by Joshua Boehm
//    email                : boehm@physics.harvard.edu
//    Last Modified        : August 10, 2004
//    VtxCluster is a reconstruction class used by VertexClusterList
//    to generate a 1-d longitudal clustering of the energy of an event

#include "VtxCluster.h"

CVSID("$Id: VtxCluster.cxx,v 1.3 2006/12/01 20:30:53 rhatcher Exp $");

// Inline Functions Defined in VtxCluster.h
// Int_t GetEdgeTrigger() {return fEdgeTrigger;}; 
// Int_t GetNCTrigger() {return fNCTrigger;};

VtxCluster::VtxCluster()
{
  fEdgeTrigger = 0;
  fNCTrigger = 0;
  fplanes.clear();
  fenergy.clear();
}

VtxCluster::VtxCluster(Int_t plane, Float_t energy_in_pe)
{
  fEdgeTrigger = 0;
  fNCTrigger = 0;
  fplanes.clear();
  fenergy.clear();
  AddPlane(plane, energy_in_pe);
}

VtxCluster::VtxCluster(const VtxCluster& input)
{
  fEdgeTrigger = input.fEdgeTrigger;
  fNCTrigger = input.fNCTrigger;
  fplanes.clear();
  fenergy.clear();

  Int_t tempPlane;
  Float_t tempE;  

  for(int i = 0; i < input.Size(); i++)
  { 
     tempPlane = input.GetPlane(i);
     tempE = input.GetEnergy(i);

     AddPlane(tempPlane,tempE);
  }

}


//Destructor
VtxCluster::~VtxCluster()
{
  //Clear the member lists
  fplanes.clear();
  fenergy.clear();
}

Float_t VtxCluster::GetEnergy(Int_t i) const
{
  if(i < 0 || i >= Size()){ cerr<<"Array out of bounds"<<endl; return 1;}
  return fenergy[i];
}

Int_t VtxCluster::GetPlane(Int_t i) const
{
  if(i < 0 || i >= Size()){ cerr<<"Array out of bounds"<<endl; return 1;}
  return fplanes[i];
}

Float_t VtxCluster::GetTotalEnergy() const
{
  Float_t total = 0;
  for(Int_t i = 0; i < Size(); i++) total += fenergy[i];
 
  return total;
}

Int_t VtxCluster::GetFirstPlane() const
{
  if(Size() == 0) return -1;
  return fplanes[0];
}

Int_t VtxCluster::GetLastPlane() const
{
  if(Size() == 0) return -1;
  return fplanes[Size()-1];
}

void VtxCluster::AddPlane(Int_t plane, Float_t pes)
{
  fplanes.push_back(plane);
  fenergy.push_back(pes);
}

Int_t VtxCluster::Size() const
{
   if(fplanes.size() != fenergy.size()) 
        cerr<<"Size: list length mismatch "<<fplanes.size()<<"  "<<fenergy.size()<<endl;
   return fplanes.size();
}

bool VtxCluster::HasRisingEdge()
{
//Look for two cases -> starting rising edge 
// some initial activity at < 5 pe and then a rising edge
  
  if(Size() < 2) return false;

  bool done = false;
  bool rising = false;
  Int_t candidate = -1;
  Float_t subTotal = 0;
  Float_t runningTotal = 0;
  bool rescind = false;

  Float_t energyCutOne = 4.0;



  for(Int_t i = 0; i < Size()-1 && !done; i++)
  {
    rising = false;
    rescind = false;
    if(fenergy[i] > energyCutOne)
    {
      runningTotal += fenergy[i];

      if(runningTotal > GetTotalEnergy()/3){
         MSG("VtxCluster",Msg::kDebug)<<"(HRE): More than 1/3 of energy "
              <<"passed, stopping"<<endl;
         i = Size();
         continue;
      }

      if(fenergy[i+1] > TMath::Max(fenergy[i]*1.05, fenergy[i] + 2.0)
              && fplanes[i+1] == fplanes[i]+1){
         //Then we have two rising planes
         rising = true;

         if(i < Size()-3)
         if(fenergy[i]+fenergy[i+1]+fenergy[i+2] < 0.5*runningTotal)
         {
            rescind = true;
         }

        if(Size() >= 12 && i > static_cast<Float_t>(Size())/2.2 
                        && runningTotal > GetTotalEnergy()/3){
           done = false;
           i = Size();
           continue;
        }

       if(!rescind){
        if(candidate != -1){ //what is the proper rising edge?
           Float_t oldTotal = 0;
           for(Int_t j = candidate; j < TMath::Min(candidate+4, Size()-1); j++)
              oldTotal += fenergy[j];
           Float_t newTotal = 0;
           for(Int_t j = i; j < TMath::Min(i+4, Size()-1); j++)
              newTotal += fenergy[j];
           if(newTotal > oldTotal) candidate = i;
         }else
            candidate = i;
         
         if(Size()-i > 4){
           subTotal = 0;
           for(Int_t j = i; j < i+4; j++) subTotal += fenergy[j];
           if(subTotal/4 > 1.5*GetTotalEnergy()/Size()){
           //Then I am satisfied that this is the beginning of
           //    something at least locally
              done = true;
              fEdgeTrigger = candidate = i;
           }
         }
        }
      }
      Int_t temp = i;
      if((rising))
      while(i < Size()-1 && fplanes[i+1] == fplanes[i]+1 && fenergy[i] > 10)
         { i++; runningTotal += fenergy[i]; }
      i = TMath::Max(temp, i--);
      if( i != temp) runningTotal -= fenergy[TMath::Min(i++, Size()-1)];
    } 
  }
 
  if(candidate != -1 && done == false) {
          fEdgeTrigger = candidate; 
          done = true;
  }

  fNCTrigger = DetermineNCTrigger();

  return done;
}

Int_t VtxCluster::DetermineNCTrigger()
{
  Int_t i = 0;
  while(i < Size() && fenergy[i] < 5) i++;
  fNCTrigger = TMath::Min(i, Size()-1);

for(Int_t j = TMath::Max(i,1) ; j < Size() -1; j++){
    if(fplanes[j]-fplanes[j-1] != 1 && fenergy[j] > 4*fenergy[fNCTrigger])
     fNCTrigger = j;
  }

  return fNCTrigger;
}

---