FillMuonId.cxx

Go to the documentation of this file.

// $Id: FillMuonId.cxx,v 1.25 2009/07/31 16:28:11 jyuko Exp $

//C++
#include <iostream>
#include <numeric>
#include <sstream>

//MINOS
#include "Calibrator/Calibrator.h"
#include "Plex/PlexStripEndId.h"
#include "Registry/Registry.h"
#include "UgliGeometry/UgliGeomHandle.h"
#include "Validity/VldContext.h"

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

//Local
#include "FillMuonId.h"

REGISTER_ANP_OBJECT(AlgSnarl,FillMuonId)
REGISTER_ANP_OBJECT(AlgEvent,FillMuonId)

using namespace std;

//---------------------------------------------------------------------------------------
Anp::FillMuonId::FillMuonId()
   :fStripWindow(4),
    fMinViewNPlane(5),
    fKeyBase(7000),
    fKeyPass(-1),
    fTimeWindow((2.0*18.83 + 0.1)*1.0e-9),
    fMinStripAdc(0.0),
    fMinStripMip(0.0),
    fADCShiftMuon(-14.0),
    fADCShiftMuonSpect(-79.0),
    fADCShiftOther(10.0),
    fADCShiftOtherSpect(0.0),
    fIsStudy(false),
    fFillDetec(true),
    fFillEvent(true),
    fFillShift(false),
    fFillShort(false)
{
}

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

//---------------------------------------------------------------------------------------
Anp::FillMuonId::PlaneData::PlaneData()
   :min_strip(-1),
    max_strip(-1),
    track_vec(),
    other_vec(),
    tpos(0.0)
{
}

//------------------------------------------------------------------------------------------
bool Anp::FillMuonId::Run(Record &record)
{
   map<short, const DataMap> track_dmap;

   for(TrackIterator track = record.TrackBegIterator(); track != record.TrackEndIterator(); ++track)
   {
      //
      // Check if KeyPass exists
      //
      if(fKeyPass > 0 && !(track -> KeyExists(fKeyPass)))
      {
         continue;
      }

      const DataMap dmap = FillMuonId::Run(*track, record);

      for(DataMap::const_iterator it = dmap.begin(); it != dmap.end(); ++it)
      {
         track -> Add(it -> first, it -> second);
      }

      if(!track_dmap.insert(map<short, const DataMap>::value_type(track -> TrackIndex(), dmap)).second)
      {
         cerr << "FillMuonId::Run - duplicate track index " << track -> TrackIndex() << endl;
      }
   }

   if(!fFillEvent || track_dmap.empty())
   {
      return true;
   }

   for(EventIterator event = record.EventBegIterator(); event != record.EventEndIterator(); ++event)
   {
      const TrackIter track = Anp::LongestTrack(*event, record);
      if(track == record.TrackEnd())
      {
         continue;
      }

      map<short, const DataMap>::const_iterator dit = track_dmap.find(track -> TrackIndex());
      if(dit == track_dmap.end())
      {
         cerr << "FillMuonId::Run - failed to find track index " << track -> TrackIndex() << endl;
         continue;
      }
      
      const DataMap &dmap = dit -> second;

      for(DataMap::const_iterator it = dmap.begin(); it != dmap.end(); ++it)
      {
         event -> Add(it -> first, it -> second);
      }     
   } 

   return true;
}

//----------------------------------------------------------------------
bool Anp::FillMuonId::Run(Event &event, const Record& record, const bool pass)
{ 
   if(!pass)
   {
      return true;
   }
   
   const TrackIter track = Anp::LongestTrack(event, record);
   if(track == record.TrackEnd())
   {
      return true;
   }

   const DataMap dmap = Run(*track, record);
   
   for(DataMap::const_iterator it = dmap.begin(); it != dmap.end(); ++it)
   {
      event.Add(it -> first, it -> second);
   }

   return true;
}

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

   int value_int = 0;
   if(reg.Get("FillMuonIdStripWindow", value_int) && value_int >= 0)
   {
      fStripWindow = value_int;
   }

   value_int = 0;
   if(reg.Get("FillMuonIdMinViewNPlane", value_int) && value_int >= 0)
   {
      fMinViewNPlane = value_int;
   }

   value_int = 0;
   if(reg.Get("FillMuonIdKeyBase", value_int) && value_int >= 0)
   {
      fKeyBase = value_int;
   }

   double value_double = 0;
   if(reg.Get("FillMuonIdTimeWindow", value_double) && value_double > 0.0)
   {
      fTimeWindow = value_double;
   }

   //
   // Configure minimum strip ADC or MIP and ADC shift for tuning
   //
   reg.Get("FillMuonIdKeyPass",            fKeyPass);
   reg.Get("FillMuonIdMinStripAdc",        fMinStripAdc);
   reg.Get("FillMuonIdMinStripMip",        fMinStripMip);
   reg.Get("FillMuonIdADCShiftMuon",       fADCShiftMuon);
   reg.Get("FillMuonIdADCShiftMuonSpect",  fADCShiftMuonSpect);
   reg.Get("FillMuonIdADCShiftOther",      fADCShiftOther);
   reg.Get("FillMuonIdADCShiftOtherSpect", fADCShiftOtherSpect);

   Anp::Read(reg, "FillMuonIdStudy", fIsStudy);
   Anp::Read(reg, "FillMuonIdDetec", fFillDetec);
   Anp::Read(reg, "FillMuonIdEvent", fFillEvent);
   Anp::Read(reg, "FillMuonIdShift", fFillShift);
   Anp::Read(reg, "FillMuonIdShort", fFillShort);

   if(reg.KeyExists("PrintConfig"))
   {
      cout << "FillMuonId::Config" << endl
           << "   StripWindow = "        << fStripWindow << endl
           << "   MinViewNPlane = "      << fMinViewNPlane << endl
           << "   KeyBase = "            << fKeyBase << endl
           << "   KeyPass = "            << fKeyPass << endl
           << "   TimeWindow = "         << fTimeWindow << endl
           << "   MinStripAdc = "        << fMinStripAdc << endl
           << "   MinStripMip = "        << fMinStripMip << endl
           << "   ADCShiftMuon = "       << fADCShiftMuon << endl
           << "   ADCShiftMuonSpect = "  << fADCShiftMuonSpect << endl
           << "   ADCShiftOther = "      << fADCShiftOther << endl
           << "   ADCShiftOtherSpect = " << fADCShiftOtherSpect << endl
           << "   IsStudy = "            << fIsStudy << endl
           << "   FillDetec = "          << fFillDetec << endl
           << "   FillEvent = "          << fFillEvent << endl
           << "   FillShift = "          << fFillShift << endl
           << "   FillShort = "          << fFillShort << endl;
   }
}

//------------------------------------------------------------------------------------------
const Anp::FillMuonId::DataMap Anp::FillMuonId::Run(const Track &track, const Record &record) const
{
   // fill variables used for muon identification
   // there are 4 variables filled:
   // 1) remove first few track planes and compute mean strip sigcor
   //    for remaining track strips
   // 2) removed first few planes from track, then
   //    and then sort these selected strips
   //    approximately in half and take ratio of lower half over 
   //    upper half
   // 3) removed first few planes from track,then for 
   //    remaining track strips compute mean strip sigcor
   //    for whole track and for last 20% of track
   //    take ratio of 20% mean over whole mean
   // 4) remove first 30% of track planes, add strips
   //    within 4 strip window and 40 ns around track
   //    take ratio of track strips' sigcor over all plane's
   //     strips' sigcor in window
   //
   // This function has hard coded key base values
   //
   

   DataMap dmap;
   
   //
   // Fill truth variable used by training algorithm(s)
   //
   short track_is_muon = -1;
   if(!record.GetHeader().IsData() && record.StdHepBeg() != record.StdHepEnd())
   {
      StdHepIter stdhep = record.FindStdHep(track);
      if(stdhep != record.StdHepEnd())
      {
         if(stdhep -> IsMuon()) track_is_muon = 1;
         else                   track_is_muon = 0;
         
         if(!dmap.insert(DataMap::value_type(fKeyBase, static_cast<float>(track_is_muon))).second)
         {
            cerr << "FillMuonId::Run - failed to insert value at " << fKeyBase << " key" << endl;
         }
      }
      else
      {
         cerr << "FillMuonId::Run - failed to find StdHep record" << endl;
         return dmap;
      }
   }

   //
   // Collect strips that we will use to compute variables...
   //
   const PMap pall = FillMuonId::Select(track, record);
   
   //
   // Track does not pass selection requirement
   //
   if(pall.empty()) return dmap;

   //
   // Fill variables using all of the detector planes
   //
   if(fFillDetec) FillMuonId::Fill(pall, dmap, 0);

   //
   // Fill number of track planes variable (once, using all detector planes)
   //
   FillMuonId::Fill0(pall.begin(), pall.end(), dmap);
   //
   // Fill tuned (ADC added or subtracted) variables
   //
   if(fFillShift)
   {      
      PMap pnew;
      if     (track_is_muon == 1) pnew = AddCharge(pall, fADCShiftMuon,  fADCShiftMuonSpect);
      else if(track_is_muon == 0) pnew = AddCharge(pall, fADCShiftOther, fADCShiftOtherSpect);

      if(!pnew.empty()) FillMuonId::Fill(pnew, dmap, 100);
   }

   //
   // Fill variables using calorimeter only strips
   // 
   const PMap pcal = FillMuonId::Select(track, record, "calor");
   if(!pcal.empty())
   {
      FillMuonId::Fill(pcal, dmap, 200);
   }   

   //
   // Fill variables using MIP calibrated strip charge
   //
   if(fFillDetec)
   {
      const PMap pmip = FillMuonId::Select(track, record, "mip");
      if(!pmip.empty())
      {
         FillMuonId::Fill(pmip, dmap, 300);
      }
   }

   //
   // Fill variables using MIP calibrated calorimeter strip charge
   //
   const PMap pmipcal = FillMuonId::Select(track, record, "calor mip");
   if(!pmipcal.empty())
   {
      FillMuonId::Fill(pmipcal, dmap, 400);
   }

   return dmap;
}

//------------------------------------------------------------------------------------------
void Anp::FillMuonId::Fill(const PMap &pmap, DataMap &dmap, const short base) const

{
   // fill variables used for muon identification
   // there are 4 variables filled:
   // 1) remove first couple track planes and compute mean strip sigcor
   //    for remaining track strips
   // 2) removed first couple planes from track, then
   //    and then sort these selected strips
   //    approximately in half and take ratio of lower half over 
   //    upper half
   // 3) removed first couple planes from track,then for 
   //    remaining track strips compute mean strip sigcor
   //    for whole track and for last 20% of track
   //    take ratio of 20% mean over whole mean
   // 4) remove first 30% of track planes, add strips
   //    within 4 strip window and 40 ns around track
   //    take ratio of track strips' sigcor over all plane's
   //     strips' sigcor in window

   PIter end_22_it = pmap.end();
   PIter end_24_it = pmap.end();

   if(fFillShort)
   {
      end_22_it = pmap.begin();
      end_24_it = pmap.begin();
   }
   else
   {
      short i_frac = 0, pdistance = 0;
      for(PIter pit = pmap.begin(); pit != pmap.end(); ++pit)
      {
         const short frac = (10 * pdistance)/pmap.size();      
         
         ++pdistance;
      
         if(pmap.size() > 9 && frac == i_frac)
         {
            continue;
         }
         
         // !!! hardcoded number !!!
         if(i_frac == 2 && end_22_it == pmap.end())
         {
            end_22_it = pit;
      }
         
         // !!! hardcoded number !!!
         if(i_frac == 4 && end_24_it == pmap.end())
         {
            end_24_it = pit;
         }
         
         ++i_frac;
      }
   }

   FillMuonId::Fill1(end_22_it, pmap.end(), dmap, base);
   FillMuonId::Fill2(end_22_it, pmap.end(), dmap, base);
   FillMuonId::Fill3(end_22_it, pmap.end(), dmap, base);
   FillMuonId::Fill4(end_24_it, pmap.end(), dmap, base);
}

//------------------------------------------------------------------------------------------
const Anp::FillMuonId::PMap Anp::FillMuonId::Select(const Track &track,
                                                    const Record &record, 
                                                    const string &opt) const
{
   //
   // Select strips matching this track, subject to various cuts. This function 
   // returns map: key is plane number and value is PlaneData struct.
   //

   PMap pmap;

   short option_detec = 0;
   if(opt.find("calor") != string::npos)
   {
      option_detec = 1;
   }
   else if(opt.find("spect") != string::npos)
   {
      option_detec = 2;
   }

   short option_view = 0;
   if(opt.find("Uview") != string::npos)
   {
      option_view = 1;
   }
   else if(opt.find("Vview") != string::npos)
   {
      option_view = 2;
   }

   short option_strip = 0;
   if(opt.find("sigmap") != string::npos)
   {
      option_strip = 1;
   }
   else if(opt.find("mip") != string::npos)
   {
      option_strip = 2;
   }

   const VldContext vldc = GetVldc(record.GetHeader());
   if(option_strip > 0)
   { 
      Calibrator::Instance().ReInitialise(vldc);
   }

   // count number of planes in each view
   int nuplane = 0, nvplane = 0;

   //
   // First loop: collect all track strips 
   //
   for(StripIter sit = record.StripBeg(); sit != record.StripEnd(); ++sit)
   {
      if(!sit -> MatchTrk(track.TrackIndex()))
      {
         continue;
      }

      if(option_detec > 0 && record.GetHeader().IsNear())
      {
         if(option_detec == 1 && sit -> GetPlane() > 120)
         {
            continue;
         }
         if(option_detec == 2 && sit -> GetPlane() < 121)
         {
            continue;
         }
      }

      if(option_view > 0)
      {
         if(option_view == 1 && !sit -> IsUview())
         {
            continue;
         }
         if(option_view == 2 && !sit -> IsVview())
         {
            continue;
         }
      }

      double strip_charge = -1.0;

      if(option_strip == 0)
      {
         // use sigcor pulse height
         if(!(sit -> SigCor() < fMinStripAdc))
         {
            strip_charge = sit -> SigCor();
         }
      }
      else
      {
         Strip::TrackInfoIter info_it = sit -> GetInfo(track.TrackIndex());
         if(info_it != sit -> InfoEndIter())
         {
            const Strip::TrackInfo &info = info_it -> second;   

            if(option_strip == 1)
            {
               // use sigmap pulse height
               if(!(info.sigmap_east + info.sigmap_west < fMinStripAdc))
               {
                  strip_charge = info.sigmap_east + info.sigmap_west;
               }
            }
            else if(option_strip == 2)
            {
               // use mip pulse height
               if(!(info.mip_east + info.mip_west < fMinStripMip))
               {
                  strip_charge = info.mip_east + info.mip_west;
               }
            }
            else
            {
               cerr << "FillMuonId::Select - unknown strip option" << endl;
            }
         }
         else
         {
            cerr << "FillMuonId::Select - failed to find Strip::TrackInfo" << endl;
         }
      }

      // strip pulse height did not pass selection cut
      if(strip_charge < 0.0)
      {
         continue;
      }

      // find plane data for this plane
      // for first call insert new PlaneData object
      PMap::iterator plane_it = pmap.find(sit -> GetPlane());
      if(plane_it == pmap.end())
      { 
         // count number of included active planes
         if(sit -> IsUview())
         {
            ++nuplane;
         }
         else
         {
            ++nvplane;
         }

         plane_it = pmap.insert(PMap::value_type(sit -> GetPlane(), PlaneData())).first;
      }

      PlaneData &data = plane_it -> second;
      
      data.tpos += sit -> TPos();
      data.track_vec.push_back(strip_charge);

      if(data.min_strip < 0 || data.min_strip > sit -> GetStrip())
      {
         data.min_strip = sit -> GetStrip();
      }
      if(data.max_strip < 0 || data.max_strip < sit -> GetStrip())
      {
         data.max_strip = sit -> GetStrip();
      }
   }

   if(option_view == 0)
   {      
      if(fFillShort)
      {  
         // check that track has minimum required number of hits in each view
         if(nuplane < fMinViewNPlane || nvplane < fMinViewNPlane)
         {
            pmap.clear();
         }

         // require at least 2 hits in each view
         if(nuplane < 2 || nvplane < 2)
         {
            pmap.clear();
         }
      }
      else
      {
         // check that track has minimum required number of hits in each view
         if(nuplane < fMinViewNPlane || nvplane < fMinViewNPlane)
         {
            pmap.clear();
         }
      }
   }
   else
   {
      if(nuplane < 2*fMinViewNPlane && nvplane < 2*fMinViewNPlane)
      {
         pmap.clear();
      }
   }

   if(pmap.empty())
   {
      return pmap;
   }

   if(option_strip > 0)
   { 
      for(PMap::iterator pit = pmap.begin(); pit != pmap.end(); ++pit)
      {
         PlaneData &data = pit -> second;
         if(data.track_vec.empty())
         {
            cerr << "FillMuonId::Select - no track strips in plane " << pit -> first << endl;
            continue;
         }

         data.tpos /= float(data.track_vec.size());
      }
   }

   const double min_time = track.GetBasic().MinTime() - fTimeWindow;
   const double max_time = track.GetBasic().MaxTime() + fTimeWindow;

   //
   // Second loop: collect all "nearby" non-track strips 
   //
   for(StripIter sit = record.StripBeg(); sit != record.StripEnd(); ++sit)
   {
      PMap::iterator pit = pmap.find(sit -> GetPlane());
      if(pit == pmap.end())
      {
         continue;
      }

      if(sit -> MatchTrk(track.TrackIndex()))
      {
         continue;
      }

      const double time = sit -> Time();
      if(time < 0.0)
      {
         cerr << "FillMuonId::Select - negative strip time: " << time << endl;
         continue;
      }

      if(time < min_time || time > max_time)
      {
         continue;
      }

      PlaneData &data = pit -> second;

      if(sit -> GetStrip() < data.min_strip - fStripWindow ||
         sit -> GetStrip() > data.max_strip + fStripWindow)
      {
         continue;
      }

      double strip_charge = -1.0;

      // use sigcor pulse height
      if(option_strip == 0)
      {
         if(!(sit -> SigCor() < fMinStripAdc))
         {
            strip_charge = sit -> SigCor();
         }
      }
      else
      {
         const float charge = Get(*sit, vldc, GetLPos(sit -> GetPlane(), pmap), option_strip);
         
         if(option_strip == 1)
         {
            // use sigmap pulse height
            if(!(charge < fMinStripAdc))
            {
               strip_charge = charge;
            }
         }
         else if(option_strip == 2)
         {
            // use mip pulse height
            if(!(charge < fMinStripMip))
            {
               strip_charge = charge;
            }
         }
      }

      // add charge for strips with pulse height above threshold
      if(!(strip_charge < 0.0))
      {
         data.other_vec.push_back(strip_charge);
      }
   }

   return pmap;
}

//------------------------------------------------------------------------------------------
const Anp::FillMuonId::PMap Anp::FillMuonId::AddCharge(const PMap &pmap,
                                                       const double calor,
                                                       const double spect) const
{
   //
   // Subtract adc values for calorimeter and spectrometer strips
   //

   PMap newmap = pmap;
   
   for(PMap::iterator pit = newmap.begin(); pit != newmap.end(); ++pit)
   {
      PlaneData &data = pit -> second;
      
      if(pit -> first > 120)
      {
         std::for_each(data.track_vec.begin(), data.track_vec.end(), Anp::Fill::AddCharge(spect));
         std::for_each(data.other_vec.begin(), data.other_vec.end(), Anp::Fill::AddCharge(spect));
      }
      else
      {
         std::for_each(data.track_vec.begin(), data.track_vec.end(), Anp::Fill::AddCharge(calor));
         std::for_each(data.other_vec.begin(), data.other_vec.end(), Anp::Fill::AddCharge(calor)); 
      }
   }

   return newmap;
}

//------------------------------------------------------------------------------------------
void Anp::FillMuonId::Fill0(PIter beg, PIter end, DataMap &dmap) const
{
  // this function computes track length in scintillator planes

   const float nplane_scint = std::distance(beg, end);

   if(nplane_scint < 1.0)
   {
      cerr << "FillMuonId::Fill0() - zero number of scintillator planes" << endl;
   }

   --end;

   const float nplane_steel = end -> first - beg -> first + 1;

   // this function computes track length in scintillator and steel planes
   if(!dmap.insert(DataMap::value_type(fKeyBase + 1, nplane_scint)).second)
   {
      cerr << "FillMuonId::Select() - failed to add data" << endl;
   }
   if(!dmap.insert(DataMap::value_type(fKeyBase + 2, nplane_steel)).second)
   {
      cerr << "FillMuonId::Select() - failed to add data" << endl;
   }  
}

//------------------------------------------------------------------------------------------
void Anp::FillMuonId::Fill1(PIter beg, PIter end, DataMap &dmap, const short base) const
{
   // this function computes mean strip sigcor per strip
   // perhaps should use per plane?

   double sum = 0.0;
   unsigned short count = 0;

   for(PIter pit = beg; pit != end; ++pit)
   {
      const PlaneData &data = pit -> second;

      sum += std::accumulate(data.track_vec.begin(), data.track_vec.end(), 0.0);

      count += data.track_vec.size();
   }

   if(count < 1)
   {
      cerr << "FillMuonId::Fill1() - track has no strips" << endl;
      return;
   }
   
   const float par = sum/double(count);

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 10, par)).second)
   {
      cerr << "FillMuonId::Fill1() - failed to add data" << endl;
   }
}

//------------------------------------------------------------------------------------------
void Anp::FillMuonId::Fill2(PIter beg, PIter end, DataMap &dmap, const short base) const
{
   // this function finds all strips within 4 strip window of track
   // these strips are sorted in increasing sigcor values
   // sigcor sum of lower ?% of strips is divided by sigcor sum 
   // of remaining strips
   
   ChargeVec svec;
   float plane_nstrip = 0.0, other_nstrip = 0.0;
   
   for(PIter pit = beg; pit != end; ++pit)
   {
      const PlaneData &data = pit -> second;

      plane_nstrip += data.track_vec.size() + data.track_vec.size();
      other_nstrip += data.other_vec.size();

      svec.insert(svec.end(), data.track_vec.begin(), data.track_vec.end());
      svec.insert(svec.end(), data.other_vec.begin(), data.other_vec.end());
   }
   
   //
   // Sort strip by pulse height
   //
   std::sort(svec.begin(), svec.end());
   
   const short size = svec.size();   
   if(size < 1)
   {
      cerr << "FillMuonId::Fill2() - strip vector is too short" << endl;
      return;
   }
   
   double lower_charge = 0.0, upper_charge = 0.0;
   double lower_nstrip = 0.0, upper_nstrip = 0.0;
   
   unsigned short count = 0;
   for(ChargeVec::const_iterator sit = svec.begin(); sit != svec.end(); ++sit)
   {
      // !!! hardcoded number !!!
      if((10 * count)/size < 4)
      {
         lower_charge += *sit;
         lower_nstrip += 1.0;
      }
      else
      {
         upper_charge += *sit;
         upper_nstrip += 1.0;
      }
      
      ++count;
   }
   
   if(upper_nstrip < 1.0 || lower_nstrip < 1.0 || short(count) != size)
   {
      cerr << "FillMuonId::Fill2() - failed to split track" << endl;
      return;
   }
   
   const double lower_mean = lower_charge/lower_nstrip;
   const double upper_mean = upper_charge/upper_nstrip;
   
   if(!(upper_mean > 0.0))
   {
      cerr << "FillMuonId::Fill2() -  upper sum is not positive" << endl;
      return;
   }
  
   const float par = lower_mean/upper_mean;

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 20, par)).second)
   {
      cerr << "FillMuonId::Fill2() - failed to add data" << endl;
   }
   
   //
   // Do not fill special variables when running production
   //
   if(!fIsStudy) return;

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 21, lower_mean)).second)
   {
      cerr << "FillMuonId::Fill2() - failed to add data" << endl;
   }
   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 22, upper_mean)).second)
   {
      cerr << "FillMuonId::Fill2() - failed to add data" << endl;
   }
   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 23, plane_nstrip)).second)
   {
      cerr << "FillMuonId::Fill4() - failed to add track data" << endl;
   }
   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 24, other_nstrip)).second)
   {
      cerr << "FillMuonId::Fill4() - failed to add track data" << endl;
   }
   if(plane_nstrip > 0.0)
   {
      if(!dmap.insert(DataMap::value_type(fKeyBase + base + 25, other_nstrip/plane_nstrip)).second)
      {
         cerr << "FillMuonId::Fill4() - failed to add track data" << endl;
      }
   }
}

//------------------------------------------------------------------------------------------
void Anp::FillMuonId::Fill3(PIter beg, PIter end, DataMap &dmap, const short base) const
{
   // this function computes mean strip sigcor in last 20% of track
   // and total mean strip sigcor
   // no strips around a track are included

   const unsigned short psize = std::distance(beg, end);
   if(psize < 1)
   {
      cerr << "FillMuonId::Fill3() - track has no strips" << endl;
      return;
   }

   double lower_charge = 0.0, upper_charge = 0.0;
   double lower_nplane = 0.0, upper_nplane = 0.0;

   unsigned short pcount = 0;
   for(PIter pit = beg; pit != end; ++pit)
   {
      const PlaneData &data = pit -> second;

      const double tcharge = std::accumulate(data.track_vec.begin(), data.track_vec.end(), 0.0);
      const double ocharge = std::accumulate(data.other_vec.begin(), data.other_vec.end(), 0.0);

      // !!! hardcoded number !!!
      if((10 * pcount)/psize < 8)
      {
         lower_charge += tcharge;
         lower_charge += ocharge;
         lower_nplane += 1.0;
      }
      else
      {
         upper_charge += tcharge;
         upper_charge += ocharge;
         upper_nplane += 1.0;
      }

      ++pcount;
   }

   if(upper_nplane < 1.0 || lower_nplane < 1.0 || pcount != psize)
   {
      cout << "upper, lower, pcount, psize = "
           << upper_nplane << ", "
           << lower_nplane << ", "
           << pcount << ", "
           << psize << endl;
      cerr << "FillMuonId::Fill3() - failed to split track" << endl;
      return;
   }
   
   const double total_mean = (lower_charge + upper_charge)/double(psize);  
   const double upper_mean = upper_charge/upper_nplane;

   if(!(total_mean > 0.0) && !(total_mean < 0.0))
   {
      cerr << "FillMuonId::Fill3() - mean charge is zero" << endl;
      return;
   }

   const float ratio = upper_mean/total_mean;

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 30, ratio)).second)
   {
      cerr << "FillMuonId::Fill3() - failed to add track data" << endl;
   }

   //
   // do not fill special variables when running production
   //
   if(!fIsStudy) return;

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 31, total_mean)).second)
   {
      cerr << "FillMuonId::Fill3() - failed to add track data" << endl;
   }
   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 32, upper_mean)).second)
   {
      cerr << "FillMuonId::Fill3() - failed to add track data" << endl;
   }
}

//------------------------------------------------------------------------------------------
void Anp::FillMuonId::Fill4(PIter beg, PIter end, DataMap &dmap, const short base) const
{
   // this function computes ratio of track sigcor to sigcor of strips near track
   // use strips within time and strip window around track

   double track_charge = 0.0, other_charge = 0.0;

   for(PIter pit = beg; pit != end; ++pit)
   {
      const PlaneData &data = pit -> second;

      track_charge += std::accumulate(data.track_vec.begin(), data.track_vec.end(), 0.0);
      other_charge += std::accumulate(data.other_vec.begin(), data.other_vec.end(), 0.0);
   }

   const double total_charge = other_charge + track_charge;
   if(!(total_charge > 0.0) && !(total_charge < 0.0))
   {
      cerr << "FillMuonId::Fill4() - total charge is zero" << endl;
      return;
   }
   
   const float par = track_charge/total_charge;

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 40, par)).second)
   {
      cerr << "FillMuonId::Fill4() - failed to add track data" << endl;
   }

   //
   // do not fill special variables when running production
   //
   if(!fIsStudy) return;

   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 41, track_charge)).second)
   {
      cerr << "FillMuonId::Fill4() - failed to add track data" << endl;
   }
   if(!dmap.insert(DataMap::value_type(fKeyBase + base + 42, other_charge)).second)
   {
      cerr << "FillMuonId::Fill4() - failed to add track data" << endl;
   }
}

//------------------------------------------------------------------------------------------
const VldContext Anp::FillMuonId::GetVldc(const Header &header) const
{
   const VldTimeStamp time(header.Sec(), header.NSec());
   if(header.IsData())
   {
      if(header.IsNear())
      {
         return VldContext(Detector::kNear, SimFlag::kData, time);
      }
      else if(header.IsFar())
      {
         return VldContext(Detector::kFar, SimFlag::kData, time);
      }
   }
   else
   {
      if(header.IsNear())
      {
         return VldContext(Detector::kNear, SimFlag::kMC, time);
      }
      else if(header.IsFar())
      {
         return VldContext(Detector::kFar, SimFlag::kMC, time);
      }
   }

   return VldContext(Detector::kUnknown, SimFlag::kUnknown, time);   
}

//------------------------------------------------------------------------------------------
float Anp::FillMuonId::GetLPos(const short plane, const PMap &pmap) const
{
   // 
   // Find longitudinal hit position using nearby planes
   //

   map<short, float> fmap;
   
   for(PMap::const_iterator pit = pmap.begin(); pit != pmap.end(); ++pit)
   {
      const short diff = plane - pit -> first;

      const PlaneData &data = pit -> second;

      // easiest solution: return tpos in the previous plane
      if(diff == -1)
      {
         return data.tpos;
      }

      // skip strips in same plane view
      if(std::abs(diff) % 2 == 0)
      {
         continue;
      }

      if(!fmap.insert(map<short, float>::value_type(diff, data.tpos)).second)
      {
         cerr << "FillMuonId::GetLPos - duplicate plane difference" << endl;
      }
   }

   if(fmap.empty())
   {
      cerr << "FillMuonId::GetLPos - no planes in orthogonal view" << endl;
      return -100.0;
   }

   map<short, float>::iterator dit_pos = fmap.begin();
   map<short, float>::iterator dit_neg = fmap.begin();
   for(map<short, float>::iterator dit = fmap.begin(); dit != fmap.end(); ++dit)
   {
      const short diff = dit -> first;
      if(diff < 0)
      {
         if(abs(dit_neg -> first) > abs(diff))
         {
            dit_neg = dit;
         }
      }
      if(diff > 0)
      {
         if(abs(dit_pos -> first) > abs(diff))
         {
            dit_pos = dit;
         }
      }
   }

   return 0.5*(dit_pos -> second + dit_neg -> second);
}

//------------------------------------------------------------------------------------------
float Anp::FillMuonId::Get(const Strip &strip, const VldContext &vldc,
                           const float lpos, const short option) const
{
   //
   // Calibrate strip pulse height from sigcor to sigmap or from sigmap to mip
   //

   float upos = -1.0, vpos = -1.0;
   if(strip.IsUview())
   {
      upos = strip.TPos();
      vpos = lpos;
   }
   else
   {
      upos = lpos;      
      vpos = strip.TPos();
   }

   UgliGeomHandle ugh(vldc);
   if(!ugh.IsValid())
   {
      cerr << "FillMuonId::Get - invalid UgliGeomHandle: " << vldc << endl;
      return 0.0;
   }
 
   const PlexStripEndId seid(strip.GetEncoded());
   const UgliStripHandle ush = ugh.GetStripHandle(seid);
  
   if(!ush.IsValid())
   {
      cerr << "FillMuonId::Get - invalid UgliStripHandle: " << seid << endl;
      return 0.0;
   }

   const TVector3 guvz(upos, vpos, strip.ZPos());
   const TVector3 lxyz = ush.GlobalToLocal(guvz, false);
   
   float sigmap = 0.0, sigmap_east = -1.0, sigmap_west = -1.0;
   if(vldc.GetDetector() == Detector::kNear)
   {
      sigmap = Calibrator::Instance().GetAttenCorrected(strip.SigCor(), lxyz.x(), seid);
   }
   else if(vldc.GetDetector() == Detector::kFar)
   {
      PlexStripEndId seid_far(seid);
      Calibrator &calibrator = Calibrator::Instance();

      if(strip.SigCorEast() > 0.0)
      {
         seid_far.SetEnd(StripEnd::kEast);
         sigmap_east = calibrator.GetAttenCorrected(strip.SigCorEast(), lxyz.x(), seid_far);
         sigmap += sigmap_west;
      }

      if(strip.SigCorWest() > 0.0)
      {      
         seid_far.SetEnd(StripEnd::kWest);
         sigmap_east = calibrator.GetAttenCorrected(strip.SigCorWest(), lxyz.x(), seid_far);
         sigmap += sigmap_east;
      }      
   }
   else
   {
      cerr << "FillMuonId::Get - unknown detector" << endl;
      return 0.0;
   }

   //
   // Just want sigmap pulse height
   //
   if(option == 1) return sigmap;

   float mip = 0.0;
   if(vldc.GetDetector() == Detector::kNear)
   {
      mip = Calibrator::Instance().GetMIP(sigmap, seid);
   }
   else if(vldc.GetDetector() == Detector::kFar)
   {
      PlexStripEndId seid_far(seid);

      if(sigmap_east > 0.0)
      {
         seid_far.SetEnd(StripEnd::kEast);
         mip += Calibrator::Instance().GetMIP(sigmap_east, seid_far);
      }

      if(sigmap_west > 0.0)
      {      
         seid_far.SetEnd(StripEnd::kWest);
         mip = Calibrator::Instance().GetMIP(sigmap_west, seid_far);
      }
   }

   return mip;
}

---