FillShortVar.cxx

Go to the documentation of this file.

// C/C++
#include <cmath>
#include <iostream>
#include <numeric>
#include <sstream>
// MINOS
#include "Registry/Registry.h"
#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 "FillShortVar.h"

REGISTER_ANP_OBJECT(AlgSnarl,FillShortVar) 

using namespace std;
//---------------------------------------------------------------------------------------
Anp::FillShortVar::FillShortVar()
  :fMinNPlane(4),
   fMaxNPlane(20),
   fEndNPlane(5),
   fKeyBase(14600),
   fTimeWindow(2.0*(18.83 + 0.1)*1.0e-9),
   fStripWindow(4),
   fMinStripADC(200.0),
   fMinStripMIP(0.3),
   fFillEvent(false),
   fDebug(false),
   fErase(false),
   fqeltrack(false),
   ftruelowE(-10),
   frecolowE(-10),
   fselectsign(0),
   fNTrackAll(0),
   fNTrackAdd(0)
{
}

//------------------------------------------------------------------------------------------
Anp::FillShortVar::~FillShortVar()
{
}
//------------------------------------------------------------------------------------------
bool Anp::FillShortVar::Run(Record &record)
{ 
  map<short,const DataMap> track_dmap;
  DataMap dmap1, dmap2,dmap3;
  DataMap  dmap;
  TrackIterator itrack = record.TrackBegIterator(); 
 
  //
  // Iterate over tracks       
  //
  while(itrack != record.TrackEndIterator())
  { 
    if(IsGoodTrack(itrack,record))
    {
      dmap.clear();
      if(FillShortVar::Study(*itrack, record,dmap1,"SIGCOR"))FillShortVar::Fill(dmap,dmap1,0);
      if(FillShortVar::Study(*itrack, record,dmap2,"SIGMAP"))FillShortVar::Fill(dmap,dmap2,100);
      if(FillShortVar::Study(*itrack, record,dmap3,"MIP"))   FillShortVar::Fill(dmap,dmap3,200);
      
      if(dmap.size()>0)
        {
          for(DataMap::const_iterator dit = dmap.begin(); dit!=dmap.end(); dit++)         
            itrack->Add(dit->first+fKeyBase,dit->second);       
          
          if(!track_dmap.insert(map<short, const DataMap>::value_type(itrack -> TrackIndex(), dmap)).second)
          {
              cerr << "FillShortVar::Run - duplicate track index " << itrack -> TrackIndex() << endl;
          }
          
          ++fNTrackAdd;
        }
      itrack++;
    }
    else{
      if(fDebug) cout<<" FillShortVar:: Run - This track is bad"<<endl;
      if(fErase) itrack = record.Erase(itrack);
      else ++itrack;
    }
    ++fNTrackAll;        
  }
   //
   // Delete events that do not match tracks
   //

  if(!fErase  && (!fFillEvent || track_dmap.empty()))  return true;

   for ( EventIterator ievent = record.EventBegIterator(); ievent !=record.EventEnd();)
   {     
     const TrackVec tvec = Anp::GetTrack(*ievent, record);
     if(fErase &&ievent->GetNTracks()>0 &&(tvec.size()==0||(int)tvec.size()!=ievent->GetNTracks() ))
       {
         ievent=record.Erase(ievent);
         continue;
       }

       if(fFillEvent) 
       {
         TrackIter titer =Anp::LongestTrack(*ievent,record);
         map<short, const DataMap>::const_iterator dit = track_dmap.find(titer -> TrackIndex());

         if(dit == track_dmap.end())
           {
             cerr << "FillShortVar::Run - failed to find track index " << titer -> TrackIndex() << endl;
             ievent++;
             continue;
           }
         const DataMap &dmap = dit -> second;

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

       }
       ievent++;
   }   
   return true;
}

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

   reg.Get("FillShortVarMinNPlane",  fMinNPlane);
   reg.Get("FillShortVarMaxNPlane",  fMaxNPlane);
   reg.Get("FillShortVarEndNPlane",  fEndNPlane);
   reg.Get("FillShortVarKeyBase",    fKeyBase);

   reg.Get("FillShortVarMinStripAdc", fMinStripADC);
   reg.Get("FillShortVarMinStripMip", fMinStripMIP);
   reg.Get("FillShortVarTimeWindow", fTimeWindow);

   //
   // Read bool and string variables
   //
                                                                                          
   Anp::Read(reg, "FillShortVarEvent",    fFillEvent);
   Anp::Read(reg, "FillShortVarDebug",    fDebug);
   Anp::Read(reg, "FillShortVarErase",    fErase);
   Anp::Read(reg, "FillShortVarQELTrack", fqeltrack);

   // Extra Configurations
   reg.Get( "FillShortVarLowTrueE", ftruelowE);
   reg.Get( "FillShortVarLowRecoE", frecolowE);
   reg.Get( "FillShortVarSelectSign", fselectsign);
   if(reg.KeyExists("PrintConfig"))
   {
      cout << "FillShortVar::Config" << endl
           << "   MinNPlane = " << fMinNPlane << endl
           << "   MaxNPlane = " << fMaxNPlane << endl
           << "   EndNPlane = " << fEndNPlane << endl
           << "   KeyBase = " << fKeyBase << endl
           << "   MinStripADC = " << fMinStripADC << endl
           << "   MinStripMIP = " << fMinStripMIP<< endl
           << "   SelectSign = " << fselectsign<< endl
           << "   TimeWindow = " << fTimeWindow << endl
           << "   Debug = " << fDebug << endl
           << "   Erase = " << fErase << endl
           << "   QELTrack  ="<< fqeltrack << endl
           << "   trueE Upper Limit  ="<< ftruelowE << endl
           << "   RecoE Upper Limit  ="<< frecolowE << endl;

   }
}
//-----------------------------------------------------------------------------
void Anp::FillShortVar::End(const DataBlock &)
{
      cout << "   NTrackAll = " << fNTrackAll << endl
           << "   NTrackAdd = " << fNTrackAdd << endl;
}
//-----------------------------------------------------------------------------
bool Anp::FillShortVar::Study(Track &track, const Record &record, DataMap &dmap,string StripMin)
{  
  if(fDebug) cout<<" FillShortVar::Study - START! with StripChargeCut "<<StripMin<<endl;
  //
  // check total number of planes in track given by track container
  //  
  double last_plane=-1; double first_plane=-1; double sum_ph=0.0;   double sum_strip =0.0;
  double trk_ph = 0.0;  double now_ph =0.0;    double now_oph=0.0;  double last_ph=0.0;
  double l_ph = 0.0;    int l_ph_plane =0;     int count_plane = 0;
  int Ibase=20;
  int Sbase=30;

  double minStripPH=FillShortVar::GetStripChargeCut(StripMin);
  if(minStripPH<0) return false;
   
  if((track.GetBasic().NPlaneScint() < fMinNPlane)||
     (track.GetBasic().NPlaneScint() > fMaxNPlane && fMaxNPlane >0)||
     (track.GetBasic().NPlaneScint() < fEndNPlane))
    {
      return false;
    }
  
  dmap[0]=minStripPH;         
  
  //==================== Get Track Map===================
   PlaneMap track_map; 
   std::map<int,float> o_mean;
   if(!GetTrackMap(track,record,track_map,o_mean,StripMin))
   {
     return false;
   }
   //===check number of planes in track with PH cut======   
   const int nplane = static_cast<int>(track_map.size());
   if((nplane < fMinNPlane) ||
      (nplane > fMaxNPlane && fMaxNPlane>0)||
      (nplane < fEndNPlane))
   {
     return false;
   }
   //===Run through track map and find Analysis variables=
   for(PlaneMap::reverse_iterator pit = track_map.rbegin(); pit != track_map.rend(); ++pit)
   {
     ++count_plane;
      const unsigned short plane = pit -> first;
      if(last_plane ==-1) last_plane= plane;
      if(plane< last_plane) first_plane = plane;
      now_ph=0.0;
   
      vector<StripIter> &tvec = pit -> second;
      for(vector<StripIter>::iterator it = tvec.begin(); it!=tvec.end(); it++)
      {
        now_ph+= TrackStripEnergy((*it),track.TrackIndex(),StripMin);   
      }

      vector<StripIter>  ovec = Get(track, plane, record,o_mean,StripMin);

      for(vector<StripIter>::iterator oit = ovec.begin(); oit!=ovec.end(); oit++)
      {       
        double wpos =GetWPos((*oit)->GetPlane(),o_mean);
        double ovec_sc= OtherStripEnergy(*oit,GetVldc(record.GetHeader()),wpos,StripMin);
        if(ovec_sc> minStripPH) now_oph +=ovec_sc;              
      }

      // If you are at the end of the track
      if(count_plane>fEndNPlane) break;
      
      sum_ph+=(now_ph + now_oph);
      trk_ph+=now_ph;      
      sum_strip +=tvec.size();            
      last_ph=now_ph;      
      dmap[count_plane] =sum_ph;
      dmap[count_plane+10]= trk_ph;      
            
      if(now_ph+now_oph>l_ph)
      { 
        l_ph=now_ph + now_oph;
        l_ph_plane=count_plane;
      }           
   }
   
   dmap[Ibase+1]  = trk_ph/sum_ph;
   dmap[Ibase+2]  = sum_ph/sum_strip;
   dmap[Ibase+3]  = l_ph_plane;
   dmap[Ibase+4]  = l_ph;


   //==== Get track scatter quantities===============
   // Straggling
   // Straightness

   float pearu=0.0;    float pearv=0.0;
   float pearmodu=0.0; float pearmodv=0.0;
   
   
   if( FillShortVar::GetScatter(track_map, pearu,pearmodu, 1)&&
       FillShortVar::GetScatter(track_map, pearv,pearmodv, 2))
     {
       if(pearu>1.00) pearu=1.0;
       if(pearv>1.00) pearv=1.0;
       if(pearmodu>1.00) pearmodu=1.0;
       if(pearmodv>1.00) pearmodv=1.0;
       
       float scatu_1 = 0.01/(1.010-pearmodu);
       float scatv_1 = 0.01/(1.010-pearmodv);
       float scatu_2 = 0.01/(1.010-pearu);
       float scatv_2 = 0.01/(1.010-pearv);
       if(scatu_1>1.0 || scatv_1>1.0 ||scatu_2>1.0 || scatv_2>1.0){
         cerr<<"FillShortVar::Study - scattering variables greater than 1.0"<<endl;
       }
       dmap[Sbase+3]  = scatu_1;     
       dmap[Sbase+5]  = scatu_2;   
       dmap[Sbase+4]  = scatv_1;
       dmap[Sbase+6]  = scatv_2;

        dmap[Sbase+7]  = track.QP();
       if(pearmodu >0 || pearmodv>0)dmap[Sbase+8]  = 2*(pearmodu*pearmodv)/(pearmodu*pearmodu+pearmodv*pearmodv);
       else dmap[Sbase+8]  = 0.5;
       dmap[Sbase+9]  = track.QP()/track.ErrorQP();
     }
   return true;
}

//-----------------------------------------------------------------------------
std::vector<Anp::StripIter> Anp::FillShortVar::Get(const Track &track,  const unsigned short plane, 
                                                   const Record &record, std::map<int,float> other_map,string StripMin) const
{
  //
  // Collect all non-track strips in a given plane within window after some given cuts   
  //  
  vector<StripIter> svec;
  int min_strip=-1;  int max_strip=-1;

  double minStripPH = FillShortVar::GetStripChargeCut(StripMin);
  if(minStripPH<0)  return svec;
  
   for(StripIter istrip = record.StripBeg(); istrip != record.StripEnd(); ++istrip)
   {
     if(istrip -> GetPlane() != plane)  continue;
     
     if(istrip -> MatchTrk(track.TrackIndex()))
     {
       if(min_strip<0 || min_strip > istrip->GetStrip()) min_strip = istrip->GetStrip();
       if(max_strip<0 || max_strip < istrip->GetStrip()) max_strip = istrip->GetStrip();
       continue;
     }
     if(fTimeWindow > 0.0)
     {
       if(istrip->Time() < track.GetBasic().MinTime() - fTimeWindow ||
          istrip->Time() > track.GetBasic().MaxTime() + fTimeWindow)
         continue;     
     }
      svec.push_back(istrip);
   }
   std::vector<Anp::StripIter>::iterator si = svec.begin() ;
   while(si!=svec.end()){
     if((*si)->GetStrip() > max_strip+fStripWindow || (*si)->GetStrip()< min_strip-fStripWindow)  si=svec.erase(si);
     else si++;
   }
   return svec;
}
//-------------------------------------------------------------------------
bool Anp::FillShortVar::GetTrackMap(Track& track, const Record& record,PlaneMap& track_map, 
                                    std::map<int,float>& other_map,string StripMin) const
{
  //
  // collect all track strips in a plane and insert plane, vector into the 'plane map'
  //  find the average position of all track (and near?) strips
  //
  double minStrip= GetStripChargeCut(StripMin);
  int  Uview=0;
  int  Vview=0;
  PlaneMap o_map;
  for(StripIter istrip = record.StripBeg(); istrip != record.StripEnd(); ++istrip)
  {

    // seperate strips in record not from track
    bool istrack = istrip -> MatchTrk(track.TrackIndex());       
    if(!istrack) continue;
    // Ignore  non calorimeter strips/events
    if(record.GetHeader().IsNear() && istrip -> GetPlane() > 120 && istrack) return false;     
    
    if(TrackStripEnergy(istrip,track.TrackIndex(),StripMin)>minStrip){  
      vector<StripIter> &svec = track_map[istrip -> GetPlane()]; 
      svec.push_back(istrip);                                           
      if(istrip->IsUview()) Uview++;
      if(istrip->IsVview()) Vview++;
    }          
  }
                                                                                                                        
  for(PlaneMap::iterator track_it =track_map.begin(); track_it!=track_map.end(); track_it++)   
    {
      float strip_mean =0.0;
      const std::vector<StripIter>  tsvec = track_it->second;
      for(std::vector<StripIter>::const_iterator it = tsvec.begin(); it!= tsvec.end(); it++){
        strip_mean += (*it) -> TPos();  
      }
      strip_mean /=(float) tsvec.size();
      other_map [track_it->first]= strip_mean;
    }

  if(Uview<=0 || Vview<=0){
    if(fDebug) cout<<" No track hits in  view  (hits in U, hits in V) = ("<<Uview<<","<<Vview<<")"<<endl;
    return false;
  }

  return true;
}
//-------------------------------------------------------------------------
bool Anp::FillShortVar::GetScatter(PlaneMap track_map, float& pear,
                                   float& dpear, int keyview) const 
{

  double TPos_mean =0.0000;  double ZPos_mean =0.0000;   double TPos2_mean =0.0000;  double ZPos2_mean =0.000000;
  double TPosZPos_mean =0.0000;

  int count=0; 

  for(PlaneMap::iterator pit = track_map.begin(); pit != track_map.end(); ++pit){
   vector<StripIter> &vec = pit->second;
    if((vec[0]->IsUview() &&!(keyview==1)) ||(vec[0]->IsVview() &&!(keyview==2))) continue; 

    for(vector<StripIter>::iterator ivec = vec.begin(); ivec!=vec.end() ; ivec++)
    {
      double tpos=(double)((*ivec))->TPos();                             
      double zpos=(double)((*ivec))->ZPos();
      TPos_mean  +=tpos;
      ZPos_mean  +=zpos;
      TPos2_mean +=pow(tpos,2);
      ZPos2_mean +=pow(zpos,2);
      TPosZPos_mean +=tpos*zpos;
      count++;
    }
  }
  if(count <=0){
    cerr<<"count is = "<<count<<endl;
    return false;
  }

  TPos2_mean/=count;
  ZPos2_mean/=count;
  TPos_mean/=count;
  ZPos_mean/=count;
  TPosZPos_mean/=count;

  if(TPos2_mean <pow(TPos_mean,2.0) || ZPos2_mean < pow(ZPos_mean,2.0))
    {
      cerr<<"FillShortVar::GetScatter count is zero or TPos,ZPos wrong "<<count<<" "<<endl;
      cerr<<"TPOS \t"<<TPos2_mean<<" "<<pow(TPos_mean,2)<<" "<<TPos2_mean-pow(TPos_mean,2)<<endl;
      cerr<<"ZPOS \t"<<ZPos2_mean<<" "<<pow(ZPos_mean,2)<<" "<<ZPos2_mean-pow(ZPos_mean,2)<<endl;
      return false;
    }

  if(count<3) {
    pear=1.0;
    dpear=1.0;
    return true;
  }

  double sdev_t=0.0;  double sdev_z=0.0;
  sdev_t = sqrt(TPos2_mean - pow(TPos_mean,2));
  sdev_z = sqrt(ZPos2_mean - pow(ZPos_mean,2));

  if(sdev_t*sdev_z>0){
    pear = abs(TPosZPos_mean - (TPos_mean*ZPos_mean))/(sdev_t*sdev_z);
  }
  else pear =0;
  if(pear>1.05){
    cerr<<"FillShortVar::GetScatter - something went terribly wrong "<<count<<" "<<pear<<endl;
    cerr<<"tpos:  "<<TPos_mean<<endl;
    cerr<<"tpos2: "<<TPos2_mean<<endl;
    cerr<<"zpos:  "<<ZPos_mean<<endl;
    cerr<<"zpos2: "<<ZPos2_mean<<endl;
    cerr<<"tzpos: "<<TPosZPos_mean<<endl;
    cerr<<"sdevt: "<<sdev_t<<endl;
    cerr<<"sdevz: "<<sdev_z<<endl;
    cerr<<"numer: "<< abs(TPosZPos_mean - (TPos_mean*ZPos_mean))<<endl;
    cerr<<"denom: "<<(sdev_t*sdev_z)<<endl;
    pear=1.0;
    dpear=1.0;
    return false;
  }

  if(sdev_t>0 && sdev_z>0) dpear  = sqrt(pow(2*pear*sdev_t*sdev_z,2) + pow((sdev_t*sdev_t)-(sdev_z*sdev_z),2))/(sdev_t*sdev_t+sdev_z*sdev_z);
  else dpear=0;

  return true;
}

                                                                                                          
float Anp::FillShortVar::GetWPos(const short plane, std::map< int,float> pmap) const
{
  //
  // Find transverse hit position using nearby planes
  //
  std::map<short, float> fmap;
  for(std::map<int,float>::iterator pit = pmap.begin(); pit != pmap.end(); ++pit)
    {     
      const short diff = plane - pit -> first;
      // easiest solution: return tpos in the previous plane      
      if(diff == -1)
        {         
          return  pit->second;
        }

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

      if(!fmap.insert(map<short, float>::value_type(diff,pit->second)).second)
        {
          cerr << "FillShortVar::GetWPos - duplicate plane difference" << endl;
        }
    }

  if(fmap.empty())
    {
      cerr << "FillShortVar::GetWPos - 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;
            }
        }
    }

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

float Anp::FillShortVar::TrackStripEnergy(Anp::StripIter si, short index,string StripMin) const
{
  if(StripMin.find("SIGCOR")!=string::npos && si->SigCor()>0) return si->SigCor();
  Strip::TrackInfoIter info_it = si->GetInfo(index);

  if(info_it == si-> InfoEndIter())
  {
    cout<<" FillShortVar::TrackStripEnergy -  Can't find info for track index  "<<index<<endl;
    return 0;
  }

  const Strip::TrackInfo &info = info_it -> second;
  double strip_charge = -1.0;

  if(StripMin.find("SIGMAP")  !=string::npos)strip_charge =(info.sigmap_east+info.sigmap_west);
  else if(StripMin.find("MIP")!=string::npos)strip_charge =(info.mip_east   +info.mip_west   );
  else cerr<<"FillShortVar::TrackStripEnergy - unable to find unit "<<StripMin<<endl;
  if(strip_charge >0) return strip_charge;
  
  cerr<<" bad strip charge "<<strip_charge<<" for "<<StripMin<<endl;
  return 0;
}

//-----------------------------------------------------------                                                                
float Anp::FillShortVar::OtherStripEnergy(StripIter sp, const VldContext &vldc,
                                          const float lpos, string StripMin) const
{
  // This is stolen from rustem's FillMuonID code.
  // Calibrate strip pulse height from sigcor to sigmap or from sigmap to mip 
  
  //
  // RETURN SIGCOR
  //
  if(StripMin.find("SIGCOR")!=string::npos) return sp->SigCor();

  const Strip& strip = *sp;
  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);
  Calibrator &calibrator = Calibrator::Instance();
  calibrator.ReInitialise(vldc);
 float sigmap = 0.0, sigmap_east = -1.0, sigmap_west = -1.0;
 if(vldc.GetDetector() == Detector::kNear)
   {
     if(fDebug) cout<<"FillShortVar::OtherStripEnergy  -near detector "<<lxyz.x()<<" "<<strip.SigCor()<<" "<<seid<<endl;
     sigmap = calibrator.GetAttenCorrected(strip.SigCor(), lxyz.x(), seid);
     if(fDebug) cout<<" here !"<<endl;
     sigmap = calibrator.GetAttenCorrectedTpos(strip.SigCor(), lpos, seid);
     if(fDebug) cout<<" after calibrator "<<endl;
   }
 else if(vldc.GetDetector() == Detector::kFar)
   {
     PlexStripEndId seid_far(seid);
     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  
 //
 //
 // RETURN SIGMAP
 //
 if(StripMin.find("SIGMAP")!=string::npos){
   if(fDebug) cout<<" FillShortVar::OtherStripEnergy - return SIGMAP "<<sigmap<<endl;
   return sigmap;
 }
 if(StripMin.find("MIP")!=string::npos){
      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
     //
     if(fDebug) cout<<" FillShortVar::OtherStripEnergy - return MIP "<<mip<<endl;
     return mip;
   }
   
   cerr<<" FillShortVar::OtherStripEnergy - Invalid type "<<StripMin<<endl;   
   return 0;
}

//-----------------------------------------------------------------------------------------
float Anp::FillShortVar::GetStripChargeCut(string StripMin) const 
{
  //
  // Given a string, find the corresponding minimum (MIP, SIGMAP or SIGCOR) in ADC or MIP
  //
  double minStripPH=-1.0;
  if(StripMin.find("MIP")!=string::npos){
    minStripPH = fMinStripMIP;
  }
  else if(StripMin.find("SIGMAP")!=string::npos){
    minStripPH= fMinStripADC;
  }
  else if(StripMin.find("SIGCOR")!=string::npos){
    minStripPH= fMinStripADC;
  }
  else {
    cerr<<" FillShortVar::GetStripChargeCut - Invalid Strip Charge Unit "<<StripMin<<endl;
  }

  return minStripPH;
}
//---------------------------------------------------------------------------- 
void Anp::FillShortVar::Fill(DataMap& dmapout, DataMap dmapin, int base){
  for(DataMap::iterator dmapit= dmapin.begin(); dmapit!=dmapin.end(); dmapit++){
    if(!dmapout.insert(DataMap::value_type(dmapit->first+ base, dmapit->second)).second){
      cerr<<" The value of "<<dmapit->first+base<<", "<< dmapit->first<<" + "<<base<<" = "<<dmapit->second
          <<" has been taken already by  "<<dmapout[dmapit->first]<<endl;
    }
  }
}
//------------------------------------------------------------------------------
std::vector<Anp::StripIter> Anp::FillShortVar::FirstPassOtherStrip(const std::vector<StripIter> instrip,const std::vector<StripIter> track_strip, double minT, double maxT) const
{
  // Do I really need this function?
  // Rustem uses (in FillMuonID) the track strips to find the orthogonal strip location for the 'other' strips (to get sigmaps)
  // keep but don't need
  std::vector<StripIter> ret_strip;

  double min_strip=-1;  double max_strip=-1;
  for( std::vector<StripIter>::const_iterator istrip = track_strip.begin(); istrip!=track_strip.end(); istrip++)
  {
    if(min_strip<0 || min_strip > (*istrip)->GetStrip()) min_strip = (*istrip)->GetStrip();
    if(max_strip<0 || max_strip < (*istrip)->GetStrip()) max_strip = (*istrip)->GetStrip();
    //add trck strips too
    ret_strip.push_back(*istrip);
  }


  for( std::vector<StripIter>::const_iterator oiter = instrip.begin(); oiter!=instrip.end();oiter++){
    // Decent SigCor
    if((*oiter)->SigCor()<fMinStripADC){
      continue;
    }
    // Within Time Window
    if(fTimeWindow > 0.0)
      {
        if((*oiter)->Time() < minT - fTimeWindow || (*oiter)->Time() > maxT + fTimeWindow)
         {
           continue;
         }
      }
    // Within Strip Window
    if((*oiter)->GetStrip() > max_strip+fStripWindow || (*oiter)->GetStrip()< min_strip-fStripWindow){
      continue;
    }

    ret_strip.push_back(*oiter);
  }
  return ret_strip;
}
//--------------------------------------------------------------------------
                                                                                                
bool Anp::FillShortVar::IsGoodTrack(TrackIter itrack, Record& record) {
  //
  // Returns whether a track is 'good' or not (study)
  //
  bool goodtrack=true;
  EventIter itevent= record.FindEvent(*itrack);

  if(     (ftruelowE>=0)   && ftruelowE<record.FindTruth(*itevent)->Energy())           goodtrack= false;
  else if((frecolowE>=0)   && (itevent!=record.EventEnd())&&frecolowE<itevent->Gev())   goodtrack= false;
  else if((fqeltrack)      && (itevent!=record.EventEnd())&&!(itevent->GetNTracks()==1))goodtrack= false;
  else if((fselectsign!=0) && itrack->QP()/fselectsign<0)                               goodtrack= false;
  else if((fMaxNPlane>0)   && itrack->GetBasic().NPlaneScint()>fMaxNPlane)              goodtrack= false;

  return goodtrack;
}

//----------------------------------------------------------------------------   

const VldContext Anp::FillShortVar::GetVldc(const Header &header) const
{
  // 
  // This is stolen from Rustem's FillMuonID code  
  //
  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);
}

---