FillSnarl.cxx

Go to the documentation of this file.

// $Id: FillSnarl.cxx,v 1.2 2008/03/31 17:25:35 rustem Exp $

// MINOS
#include "Registry/Registry.h"

// Local
#include "PhysicsNtuple/Default.h"
#include "PhysicsNtuple/Factory.h"

// Local
#include "FillSnarl.h"

REGISTER_ANP_OBJECT(AlgSnarl,FillSnarl)

using namespace std;

//---------------------------------------------------------------------------------------------
Anp::FillSnarl::FillSnarl()
   :fData(true),
    fBase(1000),
    fJump(80),
    fStep(40),
    fStop(400),
    fOver(8)
{
}

//---------------------------------------------------------------------------------------------
Anp::FillSnarl::~FillSnarl()
{
}

//---------------------------------------------------------------------------------------------
bool Anp::FillSnarl::Run(Record &record)
{
   //
   // Count various strips in this snarl and then add counts to record
   //
     
   //
   // Count each strip frequency in this snarl
   //
   vector<CountPlaneStrip> cvec;

   for(StripIter istrip = record.StripBeg(); istrip != record.StripEnd(); ++istrip)
   {
      //
      // Make counting object
      //
      CountPlaneStrip count(istrip -> GetPlane(), istrip -> GetStrip());
      
      //
      // Count once each plane number and strip number combination
      //
      if(std::find(cvec.begin(), cvec.end(), count) == cvec.end())
      {
         cvec.push_back(std::for_each(record.StripBeg(), record.StripEnd(), count));
      }
   }

   //
   // Record frequency count
   //

   vector<pair<short,short> > ignorevec;
   map<unsigned int, unsigned int> cmap;

   for(vector<CountPlaneStrip>::const_iterator cit = cvec.begin(); cit != cvec.end(); ++cit)
   {
      //
      // fOver is maximum count for overflow
      //
      const unsigned int count = std::min<unsigned int>(cit -> GetCount(), fOver);

      map<unsigned int, unsigned int>::iterator it = 
         cmap.insert(map<unsigned int, unsigned int>::value_type(count, 0)).first;
      
      ++(it -> second);

      if(cit -> GetCount() >= fOver)
      {
         ignorevec.push_back(std::pair<short, short>(cit -> GetPlane(), cit -> GetStrip()));
      }
   }
 
   for(map<unsigned int, unsigned int>::const_iterator it = cmap.begin(); it != cmap.end(); ++it)
   {
      record.Add(fBase + it -> first, static_cast<float>(it -> second));
   }

   //
   // Add total number of strips in this snarl
   //
   record.Add(fBase, static_cast<float>(record.GetNStrips()));

   //
   // Count strips above and below series of sigcor thresholds
   //

   for(int istep = fJump; istep <= fStop; istep += fStep)
   {
      CountAboveBelow count(static_cast<float>(istep), ignorevec);
      
      count = std::for_each(record.StripBeg(), record.StripEnd(), count);

      record.Add(fBase + istep + 1, static_cast<float>(count.GetNBelow()));
      record.Add(fBase + istep + 2, static_cast<float>(count.GetNAbove()));

      record.Add(fBase + istep + 3, static_cast<float>(count.GetMeanBelow()));
      record.Add(fBase + istep + 4, static_cast<float>(count.GetMeanAbove()));
   }

   
   //
   // Fill mean pulse height of ignored strips
   //
   double asum = 0.0, acount = 0.0;
   double isum = 0.0, icount = 0.0;
   
   for(StripIter istrip = record.StripBeg(); istrip != record.StripEnd(); ++istrip)
   {
      const pair<short, short> plane_strip(istrip -> GetPlane(), istrip -> GetStrip());
      
      if(std::find(ignorevec.begin(), ignorevec.end(), plane_strip) == ignorevec.end())
      {
         asum   += istrip -> SigCor();
         acount += 1.0;
      }
      else          
      {
         isum   += istrip -> SigCor();
         icount += 1.0;
      }
   }
   
   if(icount > 0.0) record.Add(fBase + fOver + 1, static_cast<float>(isum/icount));
   if(acount > 0.0) record.Add(fBase + fOver + 2, static_cast<float>(asum/acount));

   return true;
}

//-----------------------------------------------------------------------------
void Anp::FillSnarl::Config(const Registry &reg)
{
   //
   // Configure self
   //

   //
   // Fill snarl strip data
   //
   Anp::Read(reg, "FillSnarlData", fData);

   //
   // Snarl strip data base key
   //
   reg.Get("FillSnarlBase", fBase);
   reg.Get("FillSnarlJump", fJump);
   reg.Get("FillSnarlStep", fStep);
   reg.Get("FillSnarlStop", fStop);

   int value_int = 0;
   if(reg.Get("FillSnarlOver", value_int) && value_int > 0)
   {
      fOver = static_cast<unsigned int>(value_int);
   }

   assert(fJump >= 0 && fStep > 0 && fStop >= 0 && "wrong input parameters");

   if(reg.KeyExists("PrintConfig"))
   {
      cout << "FillSnarl::Config" << endl
           << "   Data = " << fData << endl
           << "   Base = " << fBase << endl
           << "   Jump = " << fJump << endl
           << "   Step = " << fStep << endl
           << "   Stop = " << fStop << endl
           << "   Over = " << fOver << endl;
   }
}

//-----------------------------------------------------------------------------
Anp::CountPlaneStrip::CountPlaneStrip()
   :plane(-1),
    strip(-1),
    count(0)    
{
}

//-----------------------------------------------------------------------------
Anp::CountPlaneStrip::CountPlaneStrip(const short plane_, const short strip_)
   :plane(plane_),
    strip(strip_),
    count(0)
{
}

//-----------------------------------------------------------------------------
void Anp::CountPlaneStrip::operator()(const Strip &strip_)
{
   assert(!(plane < 0) && !(strip < 0) && "logic error");

   if(strip_.GetPlane() == plane && strip_.GetStrip() == strip)
   {
      ++count;
   }
}

//-----------------------------------------------------------------------------
unsigned int Anp::CountPlaneStrip::GetCount() const
{
   return count;
}

//-----------------------------------------------------------------------------
short Anp::CountPlaneStrip::GetPlane() const
{
   return plane;
}

//-----------------------------------------------------------------------------
short Anp::CountPlaneStrip::GetStrip() const
{
   return strip;
}

//-----------------------------------------------------------------------------
bool Anp::operator==(const CountPlaneStrip &lhs, const CountPlaneStrip &rhs)
{
   if(lhs.GetPlane() == rhs.GetPlane() && lhs.GetStrip() == rhs.GetStrip())
   {
      return true;
   }
   
   return false;
}

//-----------------------------------------------------------------------------
Anp::CountAboveBelow::CountAboveBelow()
   :ignorevec(),
    threshold(-1.0),
    nabove(0),
    nbelow(0),
    sumabove(0.0),
    sumbelow(0.0)
{
}

//-----------------------------------------------------------------------------
Anp::CountAboveBelow::CountAboveBelow(const float threshold_, const vector<pair<short,short> > &ivec)
   :ignorevec(ivec),
    threshold(threshold_),
    nabove(0),
    nbelow(0),
    sumabove(0.0),
    sumbelow(0.0)
{
}

//-----------------------------------------------------------------------------
void Anp::CountAboveBelow::operator()(const Strip &strip)
{
   assert(!(threshold < 0.0) && "logic error");
      
   const pair<short, short> plane_strip(strip.GetPlane(), strip.GetStrip());

   if(std::find(ignorevec.begin(), ignorevec.end(), plane_strip) != ignorevec.end())
   {
      return;
   }

   if(strip.SigCor() > threshold)
   {
      ++nabove;
      sumabove += strip.SigCor();
   }
   else
   {
      ++nbelow;
      sumbelow += strip.SigCor();
   }
}

//-----------------------------------------------------------------------------
unsigned int Anp::CountAboveBelow::GetNAbove() const
{
   return nabove;
}

//-----------------------------------------------------------------------------
unsigned int Anp::CountAboveBelow::GetNBelow() const
{
   return nbelow;
}
      

//-----------------------------------------------------------------------------
double Anp::CountAboveBelow::GetMeanAbove() const
{
   if(nabove > 0)
   {
      return sumabove/static_cast<double>(nabove);
   }

   return 0.0;
}

//-----------------------------------------------------------------------------
double Anp::CountAboveBelow::GetMeanBelow() const
{
   if(nbelow > 0)
   {
      return sumbelow/static_cast<double>(nbelow);
   }

   return 0.0;
}
      

---