FillBasic.cxx

Go to the documentation of this file.

// $Id: FillBasic.cxx,v 1.14 2009/02/28 21:46:16 gmieg Exp $

// C++
#include <algorithm>
#include <cmath>

// ROOT
#include "TClonesArray.h"

// MINOS
#include "Conventions/PlaneView.h"
#include "CandNtupleSR/NtpSREvent.h"
#include "CandNtupleSR/NtpSRShower.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "DataUtil/infid.h"
#include "MessageService/MsgService.h"
#include "StandardNtuple/NtpStRecord.h"

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

// Local
#include "FillBasic.h"

CVSID("$Id: FillBasic.cxx,v 1.14 2009/02/28 21:46:16 gmieg Exp $");

using namespace std;

//-------------------------------------------------------------------
Anp::FillBasic::FillBasic()
   :fCheck(true),
    fTrackZOffset(0.0392)
{    
}

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

//----------------------------------------------------------------------
void Anp::FillBasic::Config(const Registry &reg)
{
   //
   // Configure self
   //
   Anp::Read(reg, "FillBasicCheck", fCheck);

   reg.Get("TrackZOffset", fTrackZOffset);

   if(reg.KeyExists("PrintConfigBasic"))
   {
      cout << "FillBasic::Config" << endl
           << "   Check = " << fCheck << endl
           << "   TrackZOffset = " << fTrackZOffset << endl;
   }
}

//----------------------------------------------------------------------
bool Anp::FillBasic::Fill(const NtpStRecord &ntprec, const NtpSREvent &record)
{
   //
   // Fill Anp::Basic and Anp::Vertex data
   //

   fBasic.Clear();
   fBegVtx.Clear();
   fEndVtx.Clear();

   const bool isfid = infid(ntprec.GetHeader().GetVldContext().GetDetector(),
                            ntprec.GetHeader().GetVldContext().GetSimFlag(),
                            record.vtx.x,
                            record.vtx.y,
                            record.vtx.z);

   fBasic  = FillBasic::Get(ntprec, record.stp, record.nstrip);
   fBegVtx = FillBasic::Get(record.vtx, isfid);
   
   if(!Fill(fBasic, record.ph))
   {      
      MSG("FillAlg", Msg::kWarning) << "Fill failed for NtpSREvent \n" << record;
      return false;
   }

   if(fCheck && !Check(fBasic, record.plane))
   {
      MSG("FillAlg", Msg::kWarning) << "NtpSREvent plane check failed \n" << record.plane;
      return false;
   }

   return true;
}

//----------------------------------------------------------------------
bool Anp::FillBasic::Fill(const NtpStRecord &ntprec, const NtpSRShower &record)
{
   //
   // Fill Anp::Basic and Anp::Vertex data
   //

   fBasic.Clear();
   fBegVtx.Clear();
   fEndVtx.Clear();

   const bool isfid = infid(ntprec.GetHeader().GetVldContext().GetDetector(),
                            ntprec.GetHeader().GetVldContext().GetSimFlag(),
                            record.vtx.x,
                            record.vtx.y,
                            record.vtx.z);
   
   fBasic  = FillBasic::Get(ntprec, record.stp, record.nstrip);
   fBegVtx = FillBasic::Get(record.vtx, isfid);

   if(!Fill(fBasic, record.ph))
   {      
      MSG("FillAlg", Msg::kWarning) << "Fill failed for NtpSRShower \n" << record;
      return false;
   }

   if(fCheck && !Check(fBasic, record.plane))
   {
      MSG("FillAlg", Msg::kWarning) << "NtpSRShower plane check failed \n" << record.plane;
      return false;
   }

   return true;
}

//----------------------------------------------------------------------
bool Anp::FillBasic::Fill(const NtpStRecord &ntprec, const NtpSRTrack &record)
{
   //
   // Fill Anp::Basic and Anp::Vertex data
   //

   fBasic.Clear();
   fBegVtx.Clear();
   fEndVtx.Clear();

   const bool isfid = infid(ntprec.GetHeader().GetVldContext().GetDetector(),
                            ntprec.GetHeader().GetVldContext().GetSimFlag(),
                            record.vtx.x,
                            record.vtx.y,
                            record.vtx.z - fTrackZOffset);
   

   fBasic  = FillBasic::Get(ntprec, record.stp, record.nstrip);
   fBegVtx = FillBasic::Get(record.vtx, isfid);
   fEndVtx = FillBasic::Get(record.end, false);

   if(!Fill(fBasic, record.ph))
   {      
      MSG("FillAlg", Msg::kWarning) << "Fill failed for NtpSREvent \n" << record;
      return false;
   }

   if(fCheck && !Check(fBasic, record.plane))
   {
      MSG("FillAlg", Msg::kWarning) << "NtpSRTrack plane check failed \n" << record.plane;
      return false;
   }

   return true;
}

//----------------------------------------------------------------------
const Anp::Vertex Anp::FillBasic::Get(const NtpSRVertex &vtx, const bool isfid_) const
{
   //
   // This is really just a Vertex constructor, put it here to avoid NtpSRVertex dependency 
   //

   return Vertex(vtx.u,
                 vtx.v,
                 vtx.z,
                 vtx.dcosu,
                 vtx.dcosv,
                 vtx.dcosz,
                 isfid_);
}

//----------------------------------------------------------------------
bool Anp::FillBasic::Fill(Basic &basic, const NtpSRStripPulseHeight &ph) const
{
   //
   // Fill Basic with pulse height information and check consistency between standard
   // pulse height class NtpSRStripPulseHeight and pulse height computed by iterating over
   // NtpSRStrip objects that belong to reconstructed object.
   //

   basic.siglin = ph.siglin;
   basic.pe     = ph.pe;
   basic.sigmap = ph.sigmap;
   basic.mip    = ph.mip;

   if(std::fabs(basic.adcu + basic.adcv - ph.raw) > 1.001)
   {
      static unsigned int count = 0;
      MSG("FillAlg", Msg::kWarning) << "Error " << ++count << " in NtpSRStripPulseHeight adc "
                                    << basic.adcu + basic.adcv << " - " << ph.raw << " = "
                                    << basic.adcu + basic.adcv - ph.raw << std::endl;
      return false;
   }
   if(std::fabs(basic.sigcoru + basic.sigcorv - ph.sigcor) > 1.001)
   {
      static unsigned int count = 0;
      MSG("FillAlg", Msg::kWarning) << "Error " << ++count << " in NtpSRStripPulseHeight sigcor " 
                                    << basic.sigcoru + basic.sigcorv << " - " << ph.sigcor << " = "
                                    << basic.sigcoru + basic.sigcorv - ph.sigcor << std::endl;
      return false;
   }

   return true;
}

//----------------------------------------------------------------------
bool Anp::FillBasic::Check(const Basic &basic, const NtpSRPlane &plane) const
{
   //
   // Check consistency of Basic object, this is useful bug check.
   //

   bool result = true;

   if(plane.nu + plane.nv != plane.n)
   {
      MSG("FillAlg", Msg::kWarning) << "Number of planes does not match" << endl;
      result = false;
   }
   if(basic.nuplane != short(plane.nu))
   {
      MSG("FillAlg", Msg::kWarning) << "Mismatched number of U planes" << endl;
      result = false;
   }
   if(basic.nvplane != short(plane.nv))
   {
      MSG("FillAlg", Msg::kWarning) << "Mismatched number of V planes" << endl;
      result = false;
   }
   
   if(plane.beg < plane.end)
   {
      if(basic.beg_uplane != short(plane.begu))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched beg U plane" << endl;
         result = false;
      }
      if(basic.beg_vplane != short(plane.begv))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched beg V plane" << endl;
         result = false;
      }
      if(basic.end_uplane != short(plane.endu))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched end U plane" << endl;
         result = false;
      }  
      if(basic.end_vplane != short(plane.endv))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched end V plane" << endl;
         result = false;
      }
   }
   else if(plane.beg > plane.end)
   {    
      if(basic.beg_uplane != short(plane.endu))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched beg U plane" << endl;
         result = false;
      }
      if(basic.beg_vplane != short(plane.endv))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched beg V plane" << endl;
         result = false;
      }
      if(basic.end_uplane != short(plane.begu))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched end U plane" << endl;
         result = false;
      }  
      if(basic.end_vplane != short(plane.begv))
      {
         MSG("FillAlg", Msg::kWarning) << "Mismatched end V plane" << endl;
         result = false;
      }
   }
   else
   {
      MSG("FillAlg", Msg::kWarning) << "beg and end plane are equal" << endl;
      result = false;
   }

   return result;
}

//----------------------------------------------------------------------
const Anp::Basic Anp::FillBasic::Get(const NtpStRecord &record,
                                     const int *index_array, const int nstrip) const
{  
   //
   // Iterate over strips with indexes stored in index_array and
   // compute pulse height and time and plane spans.
   //

   Basic data;

   if(nstrip < 1 || !index_array)
   {
      return data;
   }

   //
   // Initialize variables that use += operator
   //
   data.nustrip = 0;
   data.nvstrip = 0;   
   data.adcu = 0.0;
   data.adcv = 0.0;
   data.sigcoru = 0.0;
   data.sigcorv = 0.0;
   
   const TClonesArray *strip_array = record.stp;
   if(!strip_array)
   {
      MSG("FillAlg", Msg::kWarning) << "NtpStRecord does not have valid strip array." << endl; 
      return data;
   }

   vector<short> uplane_vec, vplane_vec;

   const int entries = strip_array -> GetEntries();

   const Detector::Detector_t detector = record.GetHeader().GetVldContext().GetDetector();

   bool init_time = false;
   for(int i = 0; i < nstrip; ++i)
   {            
      const short index = index_array[i];
      if(index < 0 || index >= entries)
      {
         MSG("FillAlg", Msg::kWarning) << "NtpSRStrip index is out of range: " << index << endl;
         continue; 
      }

      NtpSRStrip *ntpstp = Anp::GetSRStrip(record, index);
      if(!ntpstp)
      {
         continue;
      }
    
      if(index != ntpstp -> index)
      {  
         MSG("FillAlg", Msg::kWarning) << "Mismatched strip and TClonesArray index" << endl;
         continue;
      }

      const short plane = ntpstp -> plane;      

      if(ntpstp -> planeview == PlaneView::kU)
      {
         if(std::find(uplane_vec.begin(), uplane_vec.end(), plane) == uplane_vec.end())
         {
            uplane_vec.push_back(plane);
         }
      }
      else if(ntpstp -> planeview == PlaneView::kV)
      {
         if(std::find(vplane_vec.begin(), vplane_vec.end(), plane) == vplane_vec.end())
         {
            vplane_vec.push_back(plane);
         }
      }
      else
      {
         MSG("FillAlg", Msg::kWarning) << "Unknown planeview " << int(ntpstp -> planeview) << endl;
         continue;
      }

      if(detector == Detector::kNear)
      {
         if(ntpstp -> pmtindex0 > 0)
         {
            MSG("FillAlg", Msg::kWarning) << "In ND East PMT index is positive" << std::endl;
            continue;   
         }

         if(ntpstp -> planeview == PlaneView::kU)
         {
            ++data.nustrip;
            data.adcu += ntpstp -> ph1.raw; 
            data.sigcoru += ntpstp -> ph1.sigcor; 
         }
         else if(ntpstp -> planeview == PlaneView::kV)
         {
            ++data.nvstrip;         
            data.adcv += ntpstp -> ph1.raw;
            data.sigcorv += ntpstp -> ph1.sigcor; 
         }

         if(!init_time)
         {
            data.max_time = ntpstp -> time1;
            data.min_time = ntpstp -> time1;
            init_time = true;
         }
         else
         {
            data.max_time = std::max(ntpstp -> time1, data.max_time);
            data.min_time = std::min(ntpstp -> time1, data.min_time);
         }
      }
      else if(detector == Detector::kFar)      
      {
         if(ntpstp -> planeview == PlaneView::kU)
         {
            ++data.nustrip;

            if(ntpstp -> pmtindex0 > 0 ||
               (ntpstp -> pmtindex0 == 0 && ntpstp -> ph0.raw > 0.0))
            {
               data.adcu += ntpstp -> ph0.raw;
               data.sigcoru += ntpstp -> ph0.sigcor; 
            }
            
            if(ntpstp -> pmtindex1 > 0 ||
               (ntpstp -> pmtindex1 == 0 && ntpstp -> ph1.raw > 0.0))
            {
               data.adcu += ntpstp -> ph1.raw;
               data.sigcoru += ntpstp -> ph1.sigcor; 
            }
         }
         else if(ntpstp -> planeview == PlaneView::kV)
         {
            ++data.nvstrip;

            if(ntpstp -> pmtindex0 > 0 ||
               (ntpstp -> pmtindex0 == 0 && ntpstp -> ph0.raw > 0.0))
            {
               data.adcv += ntpstp -> ph0.raw;
               data.sigcorv += ntpstp -> ph0.sigcor; 
            }

            if(ntpstp -> pmtindex1 > 0 ||
               (ntpstp -> pmtindex1 == 0 && ntpstp -> ph1.raw > 0.0))
            {
               data.adcv += ntpstp -> ph1.raw;
               data.sigcorv += ntpstp -> ph1.sigcor; 
            }
         }


         if(ntpstp -> pmtindex0 > 0 ||
            (ntpstp -> pmtindex0 == 0 && ntpstp -> ph0.raw > 0.0))
         {
            if(!init_time)
            {
               data.max_time = ntpstp -> time0;
               data.min_time = ntpstp -> time0;
               init_time = true;
            }
            else
            {
               data.max_time = std::max(data.max_time, ntpstp -> time0);
               data.min_time = std::min(data.min_time, ntpstp -> time0);
            }
         }
         
         if(ntpstp -> pmtindex1 > 0 ||
            (ntpstp -> pmtindex1 == 0 && ntpstp -> ph1.raw > 0.0))
         {
            if(!init_time)
            {
               data.max_time = ntpstp -> time1;
               data.min_time = ntpstp -> time1;
               init_time = true;
            }
            else
            {
               data.max_time = std::max(data.max_time, ntpstp -> time1);
               data.min_time = std::min(data.min_time, ntpstp -> time1);
            }
         }

         if(!init_time)
         {
            MSG("FillAlg", Msg::kWarning) << "Failed to find initial time" << endl;
            continue;
         }
      }
      else
      {
         MSG("FillAlg", Msg::kWarning) << "Unknown detector type" << std::endl;
      }
   } 

   data.nuplane = uplane_vec.size();
   data.nvplane = vplane_vec.size();

   if(!uplane_vec.empty())
   {
      data.beg_uplane = *std::min_element(uplane_vec.begin(), uplane_vec.end());
      data.end_uplane = *std::max_element(uplane_vec.begin(), uplane_vec.end());
   }

   if(!vplane_vec.empty())
   {
      data.beg_vplane = *std::min_element(vplane_vec.begin(), vplane_vec.end());
      data.end_vplane = *std::max_element(vplane_vec.begin(), vplane_vec.end());
   }

   if(detector == Detector::kNear)
   {
      short min_plane = -1, max_plane = -1, count = 0;

      for(vector<short>::const_iterator it = uplane_vec.begin(); it != uplane_vec.end(); ++it)
      {
         const short plane = *it;
         if(plane > 120)
         {
            break;
         }
         
         ++count;

         if(min_plane < 0 || min_plane > plane)
         {
            min_plane = plane;
         }
         if(max_plane < 0 || max_plane < plane)
         {
            max_plane = plane;
         }
      }
      
      for(vector<short>::const_iterator it = vplane_vec.begin(); it != vplane_vec.end(); ++it)
      {
         const short plane = *it;
         if(plane > 120)
         {
            break;
         }
         
         ++count;        

         if(min_plane < 0 || min_plane > plane)
         {
            min_plane = plane;
         }
         if(max_plane < 0 || max_plane < plane)
         {
            max_plane = plane;
         }
      }
      
      if(max_plane - min_plane > 0)
      {
         data.active_frac = double(count)/double(1 + max_plane - min_plane);
      }
   }
   else
   {
      short min_plane = data.beg_uplane;
      short max_plane = data.end_uplane;

      if(min_plane > data.beg_vplane)
      {
         min_plane = data.beg_vplane;
      }
      if(max_plane < data.end_vplane)
      {
         max_plane = data.end_vplane;
      }
      
      const double nactive = uplane_vec.size() + vplane_vec.size();
      if(max_plane - min_plane > 0)
      {
         data.active_frac = nactive/double(max_plane - min_plane);
      }
   }

   return data;
}

//----------------------------------------------------------------------
NtpSRStrip* Anp::GetSRStrip(const NtpStRecord &record, const int index)
{
   //
   // Return pointer to NtpSRStrip with "index"
   //
   const TClonesArray *strip_array = record.stp;
   if(!strip_array)
   {
      MSG("FillAlg", Msg::kWarning) << "NtpStRecord does not have valid strip array." << endl; 
      return 0;
   }

   NtpSRStrip *ntpstp  = dynamic_cast<NtpSRStrip *>(strip_array -> At(index));
   if(!ntpstp)
   {
      MSG("FillAlg", Msg::kWarning) << "Could not find NtpSRStrip at " << index << endl;
      return 0;
   }      
   
   return ntpstp;
}

//----------------------------------------------------------------------
void Anp::PrintStrips(const NtpStRecord &record, const int *index_array, int nstrip)
{
   if(nstrip < 1 || !index_array)
   {
      return;
   }

   const TClonesArray *strip_array = record.stp;
   if(!strip_array)
   {
      cerr << "PrintStrips - NtpStRecord does not have valid strip array." << endl; 
      return;
   }

   const int entries = strip_array -> GetEntries();

   double raw_east = 0.0, raw_west = 0.0;

   cout << "Printing " << nstrip << " strips of NtpSRStrip type" << endl;

   for(int i = 0; i < nstrip; ++i)
   {            
      const short index = index_array[i];
      if(index < 0 || index >= entries)
      {
         cerr << "PrintStrips - NtpSRStrip index is out of range: " << index << endl;
         continue; 
      }

      NtpSRStrip *ntpstp = Anp::GetSRStrip(record, index);
      if(!ntpstp)
      {
         continue;
      }
    
      if(index != ntpstp -> index)
      {  
         cerr << "PrintStrips - mismatched strip and TClonesArray index" << std::endl;
         continue;
      }
      
      if(ntpstp -> pmtindex0 >= 0  && ntpstp -> ph0.raw > 0.0)
      {
         raw_east += ntpstp -> ph0.raw; 
      }
      if(ntpstp -> pmtindex1 >= 0 && ntpstp -> ph1.raw > 0.0)
      {
         raw_west += ntpstp -> ph1.raw; 
      }

      cout << "----------------------------------------------------------------------------" << endl;
      cout << "NtpSRStrip #" << i + 1 << " raw (east, west) = (" 
           <<  ntpstp -> ph0.raw << ", " << ntpstp -> ph1.raw << ")" << endl;
      ntpstp -> Print();      
   }

   cout << "total raw (east, west, total) = (" << raw_east << ", " << raw_west << ", "
        << raw_east + raw_west << ")" << endl;
}

---