MeuReco.cxx

Go to the documentation of this file.

// Coded by Jeff Hartnell Jul/2005       
//                                                                    
// Contact: j.j.hartnell@rl.ac.uk

#include <cmath>

#include "TClonesArray.h"
#include "TF1.h"
#include "TGraph.h"

#include "MessageService/MsgService.h"
#include "CandNtupleSR/NtpSRTrack.h"
#include "CandNtupleSR/NtpSRStrip.h"
#include "StandardNtuple/NtpStRecord.h"
#include "DataUtil/PlaneOutline.h"
#include "Plex/PlexHandle.h"
#include "UgliGeometry/UgliGeomHandle.h"

#include "MeuCal/MeuHitInfo.h"
#include "MeuCal/MeuReco.h"
#include "MeuCal/MeuSummary.h"

using std::endl;
using std::cout;
using std::map;
using std::vector;

CVSID("$Id: MeuReco.cxx,v 1.9 2007/05/09 14:08:38 hartnell Exp $");

//......................................................................

MeuReco::MeuReco()
{
  MSG("MeuReco",Msg::kDebug)
    <<"Running MeuReco Constructor..."<<endl;


  MSG("MeuReco",Msg::kDebug)
    <<"Finished MeuReco Constructor"<<endl;
}

//......................................................................

MeuReco::~MeuReco()
{
  MSG("MeuReco",Msg::kDebug)
    <<"Running MeuReco Destructor..."<<endl;
  

  MSG("MeuReco",Msg::kDebug)
    <<"Finished MeuReco Destructor"<<endl;
}

//......................................................................

Bool_t MeuReco::CalcPosOfGaps
(MeuSummary& s,std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //assumes that tracks don't cross the SM boundary

  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"CalcPosOfGaps: Trying to fill gaps"
    <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane
    <<", minPlane2="<<s.MinPlane2<<", maxPlane2="<<s.MaxPlane2
    <<", minPlane3="<<s.MinPlane3<<", maxPlane3="<<s.MaxPlane3
    <<endl;

  Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  Bool_t allGapsFilled=true;

  //get variables for loops
  Int_t pl=s.VtxPlane;
  Int_t tmpEndPlane=s.EndPlane;
  if (posDir) tmpEndPlane+=1;//add 1 so loop uses end plane
  else tmpEndPlane-=1;

  //first loop:
  //recalculate all z's as well as set plane numbers in gaps
  while (pl!=tmpEndPlane){
    
    MeuHitInfo& cp=plInfo[pl];

    MAXMSG("MeuReco",logLevel,1000)
      <<"first: pl="<<cp.Plane<<", cp.z="<<cp.Z
      <<endl;
    
    //recalculate all the z positions! (saw difference in csh)
    VldTimeStamp vldts;
    //VldTimeStamp vldts(ev.UnixTime,0);
    VldContext vc
      (static_cast<Detector::Detector_t>(s.Detector),
       static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
    //get the ugh
    UgliGeomHandle ugh(vc);
    PlexPlaneId plidScint
      (static_cast<Detector::Detector_t>(s.Detector),pl);
    Float_t zScint=ugh.GetScintPlnHandle(plidScint).GetZ0();
    MAXMSG("MeuReco",logLevel,1000)
      <<", zScint="<<zScint<<endl;
    //ugh-csh=0.0141 - what is this difference???
    cp.Z=zScint;
    //set the plane (wont already be set in a gap)
    cp.Plane=pl;

    //now get the view
    PlaneView::PlaneView_t pv=PlaneView::kUnknown;
    pv=ugh.GetScintPlnHandle(plidScint).GetPlaneView();
    cp.View=pv;

    if (posDir) pl++;
    else pl--;
  }

  MAXMSG("MeuReco",logLevel,1000)
    <<"CalcPosOfGaps: second loop"<<endl;

  //second loop:
  //get variables for loops
  pl=s.VtxPlane;
  if (posDir) pl+=2;
  else pl-=2;
  tmpEndPlane=s.EndPlane;
  if (posDir) tmpEndPlane-=1;//use 1 so loop uses to end-2 plane
  else tmpEndPlane+=1;

  //loop over n-4 planes in track looking for gaps to fill
  //only fill the simple 1-plane gaps
  while (pl!=tmpEndPlane){
    MeuHitInfo& cp=plInfo[pl];
    
    if (cp.TPos<=initValue) {    
      //cache the adjacent calibPos
      MeuHitInfo& cpA=plInfo[pl+2];
      MeuHitInfo& cpB=plInfo[pl-2];
      
      if (cpA.TPos<=initValue || cpB.TPos<=initValue) {
        MAXMSG("MeuReco",Msg::kVerbose,1000)
          <<"Gap too big to fill by simple average of neighbours, pl="
          <<cp.Plane<<endl;
        allGapsFilled=false;
      }
      else{
        Float_t avTPos=(cpA.TPos+cpB.TPos)*0.5;
        cp.TPos=avTPos;
      }
    }

    if (posDir) pl++;
    else pl--;
  }

  //have to check if there is a gap just in from track edges
  //since the above loop skips the outside edge 2 planes
  //determine the planes in which a gap would appear
  Int_t plNearVtx=s.VtxPlane-1;
  if (posDir) plNearVtx=s.VtxPlane+1;
  Int_t plNearEnd=s.EndPlane+1;
  if (posDir) plNearEnd=s.EndPlane-1;

  //find out what gaps are present
  Bool_t foundVtxGap=false;
  if (plInfo[plNearVtx].TPos<=initValue){
    foundVtxGap=true;
  }
  Bool_t foundEndGap=false;
  if (plInfo[plNearEnd].TPos<=initValue){
    foundEndGap=true;
  }
  if (foundVtxGap || foundEndGap) allGapsFilled=false;

  return allGapsFilled;
}

//......................................................................

Bool_t MeuReco::CalcPosOfBigGaps
(MeuSummary& s,std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //assumes that tracks don't cross the SM boundary

  MAXMSG("MeuReco",Msg::kDebug,1000)
    <<endl<<"CalcPosOfBigGaps: Trying to fill big gap"<<endl;

  Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  Bool_t allGapsFilled=true;

  //get variables for loops
  Int_t pl=s.VtxPlane;

  //loop and fill more difficult gaps  
  //this loop finds a plane with a gap then looks along the track
  //until it finds a good tpos. It then combines the tpos before
  //the gap with the one after and interpolates between them
  while (pl!=s.EndPlane){
    
    MAXMSG("MeuReco",Msg::kVerbose,10000)
      <<"pl="<<pl<<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 

    MeuHitInfo& cp=plInfo[pl];
    if (cp.TPos<=initValue) {
      //this can happen if the first but one plane is not hit
      MAXMSG("MeuReco",Msg::kDebug,1000)
        <<"plane="<<pl<<" has no tpos"
        <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl;
      allGapsFilled=false;
      if (posDir) pl++;
      else pl--;
      continue;
    }
    
    //work out next plane to check and make sure in track
    Int_t nextPlane=pl;
    if (posDir) nextPlane+=2;
    else nextPlane-=2;
    Bool_t beyondEndPlane=(nextPlane<s.EndPlane);
    if (posDir) beyondEndPlane=(nextPlane>s.EndPlane);
    if (beyondEndPlane) {
      MAXMSG("MeuReco",Msg::kDebug,1000)
        <<"Reached end of track. nextPlane="<<nextPlane
        <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
      break;
    }
    
    MeuHitInfo* cpGap=&plInfo[nextPlane];
    MeuHitInfo* cpAfterGap=0;
    
    //check if next plane is a gap
    if (cpGap->TPos<=initValue){
      
      MAXMSG("MeuReco",Msg::kDebug,1000)
        <<"Found a big gap in plane="<<nextPlane<<endl; 
      
      Bool_t foundGoodTPos=false;

      //now find plane after gap with good tpos
      for (Int_t iteration=0;iteration<100000;iteration++){
        if (posDir) nextPlane+=2;
        else nextPlane-=2;
        
        Bool_t beyondEndPlane=(nextPlane<s.EndPlane);
        if (posDir) beyondEndPlane=(nextPlane>s.EndPlane);
        if (beyondEndPlane) {
          //this can happen if the last but one plane is not hit
          MAXMSG("MeuReco",Msg::kDebug,1000)
            <<"Iterated beyond the end plane in track:"
            <<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
          allGapsFilled=false;
          break;
        }
        
        //cache the calibPos
        cpAfterGap=&plInfo[nextPlane];
        if (cpAfterGap->TPos<=initValue) {
          MAXMSG("MeuReco",Msg::kDebug,1000)
            <<"Found another gap on iteration="<<iteration
            <<", in plane="<<nextPlane<<endl; 
          continue;
        }
        else{
          MAXMSG("MeuReco",Msg::kDebug,1000)
            <<"Found next plane with good tpos, pl="<<nextPlane
            <<", iteration="<<iteration<<endl;
          foundGoodTPos=true;
          break;
        }
      }
      
      if (foundGoodTPos){
        //do interpolation
        //define m=x1_2/y1_2 = x1-x2/y1-y2  
        //(x and y are what is required)
        //define x=tpos, y=z
        MeuHitInfo& cp1=cp;
        MeuHitInfo& cp2=*cpAfterGap;
        Float_t m=(cp1.TPos-cp2.TPos)/(cp1.Z-cp2.Z);
        //Double_t x=(y-y2)*m+x2;// give above formula for m
        Float_t tpos=(cpGap->Z-cp2.Z)*m+cp2.TPos;
        if (cpGap) cpGap->TPos=tpos;
        else {
          MAXMSG("MeuReco",Msg::kWarning,1000)
            <<"Ahhhhh, not found a cp for after gap"<<endl;
          allGapsFilled=false;
        }
      
        MAXMSG("MeuReco",Msg::kDebug,1000)
          <<"Doing interpolation:"<<endl
          <<"Point1:       tpos="<<cp1.TPos<<", z="<<cp1.Z<<endl
          <<"    gap sizes:    ="<<cp1.TPos-cpGap->TPos
          <<", z="<<cp1.Z-cpGap->Z<<endl
          <<"Interp point: tpos="<<cpGap->TPos<<", z="<<cpGap->Z<<endl
          <<"    gap sizes:    ="<<cpGap->TPos-cp2.TPos
          <<",  ="<<cpGap->Z-cp2.Z<<endl
          <<"Point2:       tpos="<<cp2.TPos<<", z="<<cp2.Z<<endl
          <<endl;
        
        //Double_t m=(*ph-*(ph-1))/(*adc-*(adc-1));
        //interpPh=this->Interpolate(reqAdc,*adc,m,*ph);
        //Interpolate(Double_t y,Double_t y2,Double_t m,Double_t x2)
        //interpolation formula
        //Double_t x=(y-y2)*m+x2;
      }
    }
    else {
      MAXMSG("MeuReco",Msg::kDebug,5000)
        <<"Big gap search: no gap in plane="<<nextPlane<<endl;
    }
    
    if (posDir) pl++;
    else pl--;
  }

  MAXMSG("MeuReco",Msg::kDebug,5000)
    <<"CalcPosOfBigGaps: End of big gap search"<<endl<<endl;
  return allGapsFilled;
}

//......................................................................

Bool_t MeuReco::CalcPosOfTrkEndGaps
(MeuSummary& s,std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //assumes that tracks don't cross the SM boundary
  //also assumes that the track is long enough that there are
  //at least two good planes in each view between the track centre
  //and the two ends.
  
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"CalcPosOfTrkEndGaps: Trying to fill trk end gap"
    <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl;

  Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  Bool_t allGapsFilled=true;

  //determine the planes in which a gap would appear
  Int_t plNearVtx=s.VtxPlane-1;
  if (posDir) plNearVtx=s.VtxPlane+1;
  Int_t plNearEnd=s.EndPlane+1;
  if (posDir) plNearEnd=s.EndPlane-1;

  //find out what gaps are present
  Bool_t foundVtxGap=false;
  if (plInfo[plNearVtx].TPos<=initValue){
    foundVtxGap=true;
  }
  Bool_t foundEndGap=false;
  if (plInfo[plNearEnd].TPos<=initValue){
    foundEndGap=true;
  }

  MAXMSG("MeuReco",logLevel,1000)
    <<"foundVtxGap="<<foundVtxGap<<", plNearVtx="<<plNearVtx<<endl
    <<"foundEndGap="<<foundEndGap<<", plNearEnd="<<plNearEnd<<endl;

  //loop twice:
  //first time moving forwards to the end from the middle
  //second time moving backwards to the vtx from the middle
  for (Int_t i=0;i<2;i++){
    Bool_t movingTowardsEnd=true;
    if (i==1) movingTowardsEnd=false;

    //check if a gap exists at each end
    if (movingTowardsEnd){
      if (!foundEndGap) {
        MAXMSG("MeuReco",logLevel,1000)
          <<"No gap at end of track, so not entering loop"<<endl;
        continue;
      }
    }
    else{
      if (!foundVtxGap) {
        MAXMSG("MeuReco",logLevel,1000)
          <<"No gap at vtx of track, so not entering loop"<<endl;
        continue;
      }
    }

    MeuHitInfo* cpA=0;
    MeuHitInfo* cpB=0;
    
    //start in the middle of the event
    Int_t pl=(s.VtxPlane+s.EndPlane)/2;
    //make sure that the pl starts in the right view
    if (!movingTowardsEnd){
      if (pl%2!=plNearVtx%2) pl++;
    }
    else{
      if (pl%2!=plNearEnd%2) pl++;
    }
    
    //calc the planes beyond the track ends so loop stops in right
    //place
    Int_t plBeyondEnd=plNearEnd-2;
    if (posDir) plBeyondEnd=plNearEnd+2;
    Int_t plBeyondVtx=plNearVtx+2;
    if (posDir) plBeyondVtx=plNearVtx-2;
    
    Int_t lastLoopPlane=plBeyondEnd;
    if (!movingTowardsEnd) lastLoopPlane=plBeyondVtx;
    
    MAXMSG("MeuReco",logLevel,1000)
      <<"Starting on pl="<<pl
      <<" and working towards pl="<<lastLoopPlane
      <<", movingTowardsEnd="<<movingTowardsEnd<<endl;
    
    //loop and fill more difficult gaps  
    //move towards end plane from middle first
    while (pl!=lastLoopPlane){
    
      MAXMSG("MeuReco",logLevel,1000)
        <<"CalcPosOfTrkEndGaps: pl="<<pl
        <<", track:"<<s.VtxPlane<<"->"<<s.EndPlane<<endl; 

      MeuHitInfo* cp=&plInfo[pl];
      if (cp->TPos<=initValue){
        MAXMSG("MeuReco",logLevel,1000)
          <<"Found gap at pl="<<pl<<endl;
        allGapsFilled=false;

        //if there are two good tpos then can interpolate
        if (cpA!=0 && cpB!=0){
        
          //do interpolation
          //define m=x1_2/y1_2 = x1-x2/y1-y2  
          //(x and y are what is required)
          //define x=tpos, y=z
          MeuHitInfo& cp1=*cpA;
          MeuHitInfo& cp2=*cpB;
          MeuHitInfo* cpGap=cp;
          Float_t m=(cp1.TPos-cp2.TPos)/(cp1.Z-cp2.Z);
          //Double_t x=(y-y2)*m+x2;// give above formula for m
          Float_t tpos=(cpGap->Z-cp2.Z)*m+cp2.TPos;
          if (cpGap) {
            cpGap->TPos=tpos;
            //record result
            allGapsFilled=true;
          }
          else {
            MAXMSG("MeuReco",Msg::kWarning,1000)
              <<"Ahhhhh, not found a cp for after gap"<<endl;
            allGapsFilled=false;
          }
          
          MAXMSG("MeuReco",logLevel,1000)
            <<"Doing interpolation:"<<endl
            <<"Point1:       tpos="<<cp1.TPos<<", z="<<cp1.Z<<endl
            <<"    gap sizes:    ="<<cp1.TPos-cpGap->TPos
            <<", z="<<cp1.Z-cpGap->Z<<endl
            <<"Interp point: tpos="<<cpGap->TPos
            <<", z="<<cpGap->Z<<endl
            <<"    gap sizes:    ="<<cpGap->TPos-cp2.TPos
            <<",  ="<<cpGap->Z-cp2.Z<<endl
            <<"Point2:       tpos="<<cp2.TPos<<", z="<<cp2.Z<<endl
            <<endl;

          VldTimeStamp vldts;
          //VldTimeStamp vldts(ev.UnixTime,0);
          VldContext vc
            (static_cast<Detector::Detector_t>(s.Detector),
             static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
          //get the ugh
          UgliGeomHandle ugh(vc);
          //Float_t zmin=-1;
          //Float_t zmax=-1;
          Float_t tmin=-1;
          Float_t tmax=-1;
          PlexPlaneId plidScint
            (static_cast<Detector::Detector_t>(s.Detector),
             pl);
          PlaneView::PlaneView_t pv=PlaneView::kUnknown;
          pv=ugh.GetScintPlnHandle(plidScint).GetPlaneView();
          string sPv="";//PlaneView::AsString(pv);
          //ugh.GetZExtent(zmin,zmax);
          ugh.GetTransverseExtent(pv,tmin,tmax);
          MAXMSG("MeuReco",Msg::kVerbose,1000)
            <<"pl="<<pl<<", view="<<sPv<<endl
            //<<"Got ZExtent, min="<<zmin<<", max="<<zmax<<endl
            <<"Got TransverseExtent, min="<<tmin<<", max="<<tmax<<endl;
          
          if (tpos>tmax || tpos<tmin){
            Float_t tposToUse=tmax;
            if (tpos<0) tposToUse=tmin;
            
            MAXMSG("MeuReco",logLevel,1000)
              <<endl<<"Interpolated out of the detector: tpos="<<tpos
              <<", capping the value at "<<tposToUse<<endl
              <<"TransverseExtent: min="<<tmin<<", max="<<tmax
              <<endl;
            
            cpGap->TPos=tposToUse;
            allGapsFilled=true;
          }
        }
      }
      else{//use this as a tpos
        MAXMSG("MeuReco",logLevel,1000)
          <<"Found good tpos in plane="<<pl<<endl;
        cpB=cpA;//overwrite previous cpB with later one
        cpA=cp;
      }
    
      if (movingTowardsEnd) {
        //forwards loop
        if (posDir) pl+=2;//move to next in view
        else pl-=2;
      }
      else{
        //backwards loop
        if (posDir) pl-=2;//move to next in view
        else pl+=2;
      }
    }
  }

  Bool_t warningUnfilledGaps=false;

  if (!allGapsFilled) {
    MAXMSG("MeuReco",Msg::kDebug,1000)
      <<"Failed to fill all gaps"<<endl;
  }
  else{//loop and check all gaps were indeed filled
    //get variables for loops
    Int_t pl=s.VtxPlane;
    Int_t tmpEndPlane=s.EndPlane;
    if (posDir) tmpEndPlane+=1;//add 1 so loop uses end plane
    else tmpEndPlane-=1;
    
    //first loop:
    //recalculate all z's as well as set plane numbers in gaps
    while (pl!=tmpEndPlane){
      MeuHitInfo& cp=plInfo[pl];     
      if (cp.TPos<=initValue){
        MAXMSG("MeuReco",Msg::kWarning,20)
          <<endl<<endl<<"**** Found unfilled gap! ****"<<endl<<endl;
        warningUnfilledGaps=true;
        allGapsFilled=false;
      }
      
      if (posDir) pl++;
      else pl--;
    }
  }

  if (warningUnfilledGaps){
    MAXMSG("MeuReco",Msg::kInfo,20)
      <<"Event="<<s.Event<<", run="<<s.Run<<", fSnarl="<<s.Snarl
      <<endl;
    Int_t pl=s.VtxPlane;
    Int_t tmpEndPlane=s.EndPlane;
    if (posDir) tmpEndPlane+=1;//add 1 so loop uses end plane
    else tmpEndPlane-=1;

    //print out plane gap info
    while (pl!=tmpEndPlane){
      MeuHitInfo& cp=plInfo[pl];     
      MAXMSG("MeuReco",Msg::kInfo,100)
        <<"pl="<<cp.Plane<<", st="<<cp.Strip<<", view="<<cp.View
        <<", sigCor="<<cp.SigCor
        <<", t="<<cp.TPos<<", l="<<cp.LPos<<endl;
      
      if (posDir) pl++;
      else pl--;
    }
  }

  MAXMSG("MeuReco",logLevel,1000)
    <<"CalcPosOfTrkEndGaps: End of trk end gap search"<<endl;
  return allGapsFilled;
}

//......................................................................

void MeuReco::CalcPLCor(MeuSummary& s,
                         std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  Msg::LogLevel_t logLevel2=Msg::kDebug;//for finding +/- planes

  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"CalcPLCor: Get path length correction, track: "
    <<s.VtxPlane<<"->"<<s.EndPlane<<endl;
  
  Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  static Int_t event=0;
  event++;//for getting MC PLCor

  //get variables for loops
  Int_t pl=s.VtxPlane;
  Int_t tmpEndPlane=s.EndPlane;
  if (posDir) tmpEndPlane++;//make sure last plane is used
  else tmpEndPlane--;
  
  //loop and determine PLCor in each plane
  while (pl!=tmpEndPlane){    
    MAXMSG("MeuReco",Msg::kVerbose,1000)
      <<"pl="<<pl<<endl;

    //cache the current cp
    MeuHitInfo& cp=plInfo[pl];

    //first step is to determine the planes to use in determining PLCor

    //work out the plane to use that is towards the track end
    Int_t nextPlane=pl;
    if (posDir) nextPlane+=2;
    else nextPlane-=2;

    Bool_t beyondEndPlane=(nextPlane<s.EndPlane);
    if (posDir) beyondEndPlane=(nextPlane>s.EndPlane);
    if (beyondEndPlane) {
      MAXMSG("MeuReco",logLevel2,1000)
        <<"pl+2="<<nextPlane<<", gone beyond the end plane in track: "
        <<s.VtxPlane<<"->"<<s.EndPlane<<endl; 

      //try looking only 1 plane further
      nextPlane=pl;
      if (posDir) nextPlane+=1;
      else nextPlane-=1;
      beyondEndPlane=(nextPlane<s.EndPlane);
      if (posDir) beyondEndPlane=(nextPlane>s.EndPlane);

      if (beyondEndPlane) {
        MAXMSG("MeuReco",logLevel2,1000)
          <<"pl+1="<<nextPlane<<", gone beyond the end plane in track: "
          <<s.VtxPlane<<"->"<<s.EndPlane<<endl;
        nextPlane=pl;//must be at one track end so just take pl
      }
      else{
        MAXMSG("MeuReco",logLevel2,1000)
          <<"pl+1="<<nextPlane<<", is ok for track: "
          <<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
      }
    }

    //work out the plane to use that is towards the track vtx
    Int_t prevPlane=pl;
    if (posDir) prevPlane-=2;
    else prevPlane+=2;

    Bool_t beyondVtxPlane=(prevPlane>s.VtxPlane);
    if (posDir) beyondVtxPlane=(prevPlane<s.VtxPlane);
    if (beyondVtxPlane) {
      MAXMSG("MeuReco",logLevel2,1000)
        <<"pl+2="<<prevPlane<<", gone beyond the end plane in track: "
        <<s.VtxPlane<<"->"<<s.EndPlane<<endl; 

      //try looking only 1 plane further
      prevPlane=pl;
      if (posDir) prevPlane-=1;
      else prevPlane+=1;
      beyondVtxPlane=(prevPlane>s.VtxPlane);
      if (posDir) beyondVtxPlane=(prevPlane<s.VtxPlane);

      if (beyondVtxPlane) {
        MAXMSG("MeuReco",logLevel2,1000)
          <<"pl+1="<<prevPlane<<", gone beyond the end plane in track: "
          <<s.VtxPlane<<"->"<<s.EndPlane<<endl;
        prevPlane=pl;//must be at one track end so just take pl
      }
      else{
        MAXMSG("MeuReco",logLevel2,1000)
          <<"pl+1="<<prevPlane<<", is ok for track: "
          <<s.VtxPlane<<"->"<<s.EndPlane<<endl; 
      }
    }

    MAXMSG("MeuReco",logLevel2,1000)
      <<"For pl="<<pl<<" using planes: "<<prevPlane<<"->"<<nextPlane
      <<endl;

    //second step is to create a TGraph of 3-5 planes and then 
    //fit to a straight line

    TGraph gx_z(3);//min of 3 points
    TGraph gy_z(3);
    Int_t tmpPl=prevPlane;
    Int_t tmpPlEnd=nextPlane-1;
    if (posDir) tmpPlEnd=nextPlane+1;
    Int_t i=0;
    //loop over planes to fit
    while (tmpPl!=tmpPlEnd){
      MeuHitInfo& tmpcp=plInfo[tmpPl];

      Float_t x=tmpcp.X;
      Float_t y=tmpcp.Y;
      Float_t z=tmpcp.Z;
      gx_z.SetPoint(i,z,x);
      gy_z.SetPoint(i,z,y);
      MAXMSG("MeuReco",Msg::kVerbose,1000)
        <<"Adding points to graphs: x="<<x<<", y="<<y<<", z="<<z<<endl;
      //count the number of points added
      i++;
      //increment the plane
      if (posDir) tmpPl++;
      else tmpPl--;
    }

    //calc the MC plcor if this is MC
    Double_t PLCorMC=-1;
    Double_t pathLength=-1;//=this->TrueMainPL(pl,event);
    if (pathLength>=0){//-1 if no MC
      static Double_t lastPL=1.3*Munits::cm;
      //make a best guess if a zero is returned
      if (pathLength<=0) pathLength=lastPL;
      lastPL=pathLength;
      
      PLCorMC=pathLength/(0.95*Munits::cm);
    }

    Float_t PLCor=initValue;
    if (i>=3){
      gx_z.Fit("pol1","q");
      gy_z.Fit("pol1","q");
      Float_t mx_z=gx_z.GetFunction("pol1")->GetParameter(1);
      Float_t my_z=gy_z.GetFunction("pol1")->GetParameter(1);

      //calc plcor; don't need to worry about 0.95 since
      //it's already per unit z
      PLCor=sqrt(pow(mx_z,2)+pow(my_z,2)+1);
      MAXMSG("MeuReco",logLevel,1000)
        <<"pl="<<pl<<", reconstructed PLCor="<<PLCor
        <<", MC PLCor="<<PLCorMC
        <<endl;
      //set the value in cp
      cp.PLCor=PLCor;
    }
    else {
      MSG("MeuReco",Msg::kWarning)
        <<"Less than 3 points, can't do local grad. fit"<<endl;
    }
    
    if (posDir) pl++;
    else pl--;
  }
}

//......................................................................

void MeuReco::CalcMEU(MeuSummary& s,
                       std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"CalcMEU: Get the MEU, track: "
    <<s.VtxPlane<<"->"<<s.EndPlane<<endl;

  //Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  //static Int_t event=0;
  //event++;//for getting MC PLCor

  Float_t meuSigLinOnly=0;
  Float_t meuSigDrf=0;
  Float_t meuAdc=0;
  Float_t meuPe=0;
  Float_t meuSigLin=0;
  Float_t meuSigCor=0;
  Float_t meuSigMap=0;
  Float_t meuNumDigits=0;
  Float_t meuNumStrips=0;
  Float_t meuEnDep=0;
  Float_t meuSigMap_MeV=0;
  Float_t meuSigCor_MeV=0;
  Float_t sigCorSum=0;
  Float_t sigMapSum=0;
  Float_t enDepSum=0;
  Float_t materialFromTrkEnd=0;
  Float_t materialInWindow=0;
  Int_t winCounter=0;
  Float_t avPLCor=0;
  Int_t winStopSidePl=-1;
  Int_t winVtxSidePl=-1;
  //const Float_t material01Pl=(2.5+1)*Munits::cm;
  const Float_t material01Pl=(5.94)*Munits::cm;//including air gap
  const Float_t material16Pl=16*material01Pl;//distance to trk end
  const Float_t material14Pl=14*material01Pl;//length of window
  MAXMSG("MeuReco",Msg::kInfo,1)
    <<endl<<"Using 14 pls material="<<material14Pl
    <<", 16 plns="<<material16Pl
    <<", 1 plns="<<material01Pl<<endl;

  //calculate energy deposition in the window
  if (s.EndPlane>=0){
    
    //get variables for loops
    //NOTE: start at the end of the track NOT the vtx
    Int_t pl=s.EndPlane;
    Int_t tmpVtxPlane=s.VtxPlane;
    if (posDir) tmpVtxPlane--;//make sure last plane is used
    else tmpVtxPlane++;
    
    //loop over the track starting at the end
    while (pl!=tmpVtxPlane){
      
      //cache the current cp
      MeuHitInfo& cp=plInfo[pl];
      Float_t PLCor=cp.PLCor;
      Float_t plMaterial=PLCor*material01Pl;
      
      //sum up energy in window
      //require that the material from trk end is >16pl already 
      //or when you include the current plane
      //ALSO check if already passed window
      if ((materialFromTrkEnd>material16Pl || 
           materialFromTrkEnd+plMaterial>material16Pl) && 
          materialInWindow<material14Pl){
        
        //record the first plane included in window
        if (winStopSidePl<0) {
          winStopSidePl=pl;
        }
        winVtxSidePl=pl;//the last plane will be the last time called
        
        //add the energy in this plane to the total
        materialInWindow+=plMaterial;

        //sum the path length corrections
        avPLCor+=PLCor;

        //sum the energy deposition in the window
        meuSigLinOnly+=cp.SigLinOnly/PLCor;
        meuSigDrf+=cp.SigDrf/PLCor;
        meuAdc+=cp.Adc/PLCor;
        meuPe+=cp.Pe/PLCor;
        meuSigLin+=cp.SigLin/PLCor;
        meuSigCor+=cp.SigCor/PLCor;
        meuSigMap+=cp.SigMap/PLCor;
        meuNumDigits+=cp.NumDigits;
        meuNumStrips+=cp.NumStrips;

        if (cp.MCEnDep>0) meuEnDep+=cp.MCEnDep/PLCor;

        sigCorSum+=cp.SigCor;
        sigMapSum+=cp.SigMap;
        enDepSum+=cp.MCEnDep;
        winCounter++;
      }

      materialFromTrkEnd+=plMaterial;

      MAXMSG("MeuReco",logLevel,1000)
        <<"plMat="<<plMaterial
        <<", matTrk="<<materialFromTrkEnd<<" (of "<<material16Pl<<")"
        <<", matWin="<<materialInWindow<<" (of "<<material14Pl<<")"
        <<endl
        <<"meuSigMap="<<meuSigMap<<", PLCor="<<PLCor<<endl;
      
      //increment the plane
      if (posDir) pl--;//working backwards from end to vtx
      else pl++;
    }
  }
  else{
    MAXMSG("MeuReco",Msg::kWarning,1000)
      <<"Can't CalcMEU: fEndPlane="<<s.EndPlane
      <<", vtxPl="<<s.VtxPlane<<endl;
  }

  //store the material traversed - do this for all events
  s.TotalMatTraversed=materialFromTrkEnd;
  
  //if a window was completed
  if (materialInWindow>material14Pl){

    if (winCounter>0){
      meuSigLinOnly/=winCounter;
      meuSigDrf/=winCounter;
      meuAdc/=winCounter;
      meuPe/=winCounter;
      meuSigLin/=winCounter;
      meuSigCor/=winCounter;
      meuSigMap/=winCounter;
      meuNumDigits/=winCounter;
      meuNumStrips/=winCounter;
      meuEnDep/=winCounter;

      if (enDepSum>0){
        meuSigMap_MeV=sigMapSum/(enDepSum/Munits::MeV);
        meuSigCor_MeV=sigCorSum/(enDepSum/Munits::MeV);
      }

      avPLCor/=winCounter;
      
      //store the meu values
      s.WinSigLinOnly=meuSigLinOnly;
      s.WinSigDrf=meuSigDrf;
      s.WinAdc=meuAdc;
      s.WinPe=meuPe;
      s.WinSigLin=meuSigLin;
      s.WinSigCor=meuSigCor;
      s.WinSigMap=meuSigMap;
      s.WinAvNumDigits=meuNumDigits;
      s.WinAvNumStrips=meuNumStrips;
      s.MCWinEnDep=meuEnDep;
      
      s.MCWinSigCor_MeV=meuSigCor_MeV;
      s.MCWinSigMap_MeV=meuSigMap_MeV;
      
      //store others
      s.WinAvPLCor=avPLCor;
      s.WinAvCosThetaZ=1./s.WinAvPLCor;
      if (!posDir) s.WinAvCosThetaZ*=-1;//make it negative
      s.WinStopSidePl=winStopSidePl;
      s.WinVtxSidePl=winVtxSidePl;
    }
  }
  //else they just stay as the defaults
}

//......................................................................

void MeuReco::PrintMeuHitInfo
(const std::map<Int_t,MeuHitInfo>& plInfo) const
{
  Msg::LogLevel_t logLevel=Msg::kInfo;
    
  for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.begin();
       it!=plInfo.end();++it){
        
    const MeuHitInfo& cp=it->second;
        
    MAXMSG("MeuReco",logLevel,1000)
      <<"Print: pl="<<cp.Plane<<", st="<<cp.Strip
      <<", view="<<cp.View
      <<", LPos="<<cp.LPos
      <<", TPos="<<cp.TPos
      <<endl;
  }
}

//......................................................................

void MeuReco::PrintMeuHitInfoFull
(const std::map<Int_t,MeuHitInfo>& plInfo) const
{
  Msg::LogLevel_t logLevel=Msg::kInfo;
    
  MSG("MeuReco",logLevel)<<"Printing full MeuHitInfo:"<<endl;

  for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.begin();
       it!=plInfo.end();++it){
        
    const MeuHitInfo& cp=it->second;
        
    MSG("MeuReco",logLevel)
      <<"pl="<<cp.Plane<<", st="<<cp.Strip
      <<", view="<<cp.View
      <<", T="<<cp.TPos<<", L="<<cp.LPos
      <<", PLCor="<<cp.PLCor
      <<", xyz="<<cp.X<<","<<cp.Y<<","<<cp.Z
      <<endl
      <<"ADC: "<<cp.Adc<<" ("<<cp.Adc1<<","<<cp.Adc2<<")"
      <<", SigLin: "<<cp.SigLin
      <<" ("<<cp.SigLin1<<","<<cp.SigLin2<<")"
      <<endl
      <<"PE: "<<cp.Pe<<" ("<<cp.Pe1<<","<<cp.Pe2<<")"
      <<endl
      <<"SigCor: "<<cp.SigCor
      <<" ("<<cp.SigCor1<<","<<cp.SigCor2<<")"
      <<"SigCorTrk: "<<-1
      <<" ("<<cp.SigCorTrk1<<","<<cp.SigCorTrk2<<")"
      <<endl
      <<"SigMap: "<<cp.SigMap
      <<" ("<<cp.SigMap<<","<<cp.SigMap<<")"
      <<", MCSigMap: "<<cp.MCSigMap
      <<endl
      <<"MCEnDep="<<cp.MCEnDep
      <<", MCSuppEnDep="<<cp.MCSuppEnDep
      <<endl;
  }
}

//......................................................................

void MeuReco::PrintMeuSummary(const MeuSummary& s) const
{
  Msg::LogLevel_t logLevel=Msg::kInfo;
    
  MSG("MeuReco",logLevel)
    <<"Printing full MeuHitInfo:"<<endl
    <<"fCount"<<s.Count<<endl
    <<"fTemperature"<<s.Temperature<<endl
        
    <<"fEvent"<<s.Event<<endl
    <<"fRun"<<s.Run<<endl
    <<"fSubRun"<<s.SubRun<<endl
    <<"fTimeSec"<<s.TimeSec<<endl
    <<"fSnarl"<<s.Snarl<<endl
    <<"fNStrip"<<s.NStrip<<endl
        
    <<"fGoodContainment"<<s.GoodContainment<<endl
    <<"fGoodWindow"<<s.GoodWindow<<endl
        
    <<"fSM1"<<s.SM1<<endl
    <<"fSM2"<<s.SM2<<endl
    <<"fSMBoth"<<s.SMBoth<<endl
        
    <<"fEntSM1Front"<<s.EntSM1Front<<endl
    <<"fEntSM1Back"<<s.EntSM1Back<<endl
    <<"fEntSM2Front"<<s.EntSM2Front<<endl
    <<"fEntSM2Back"<<s.EntSM2Back<<endl
    <<"fEntSMEnd"<<s.EntSMEnd<<endl
        
    <<"fExitSM1Front"<<s.ExitSM1Front<<endl
    <<"fExitSM1Back"<<s.ExitSM1Back<<endl
    <<"fExitSM2Front"<<s.ExitSM2Front<<endl
    <<"fExitSM2Back"<<s.ExitSM2Back<<endl
    <<"fExitSMEnd"<<s.ExitSMEnd<<endl
        
    <<"fPC"<<s.PC<<endl
    <<"fFC"<<s.FC<<endl
    <<"fAwayFromCoil"<<s.AwayFromCoil<<endl
        
    <<"fMCHighEn"<<s.MCHighEn<<endl
    <<"fMCLowEn"<<s.MCLowEn<<endl
    <<"fMCParticleId"<<s.MCParticleId<<endl
    <<"fRFid"<<s.RFid<<endl
    <<"fRFidCoil"<<s.RFidCoil<<endl
    <<"fRVtx"<<s.RVtx<<endl
    <<"fREnd"<<s.REnd<<endl
    <<"fVtxPlane"<<s.VtxPlane<<endl
    <<"fEndPlane"<<s.EndPlane<<endl
    <<"fMCVtxPlane"<<s.MCVtxPlane<<endl
    <<"fMCEndPlane"<<s.MCEndPlane<<endl;
}

//......................................................................

Bool_t MeuReco::CalcLPos(MeuSummary& s,
                          std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  
  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"CalcLPos: Get position along strip, track: "
    <<s.VtxPlane<<"->"<<s.EndPlane<<endl;

  Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  
  //get variables for loops
  Int_t pl=s.VtxPlane;
  if (posDir) pl+=1;//start one plane in from vtx
  else pl-=1;
  Int_t tmpEndPlane=s.EndPlane;

  //loop over n-4 planes in track looking for gaps to fill
  //only fill the simple 1-plane gaps
  while (pl!=tmpEndPlane){//does not actually use tmpEndPlane
    MeuHitInfo& cp=plInfo[pl];
    
    if (cp.LPos<=initValue) {    
      //cache the adjacent calibPos
      MeuHitInfo& cpA=plInfo[pl+1];
      MeuHitInfo& cpB=plInfo[pl-1];
      
      if (cpA.TPos<=initValue || cpB.TPos<=initValue) {
        MAXMSG("MeuReco",Msg::kWarning,1000)
          <<"Found a gap, pl="<<cp.Plane
          <<", pl+1 tposA="<<cpA.TPos
          <<", pl-1 tposB="<<cpB.TPos
          <<endl;
      }
      else{
        Float_t avTPos=(cpA.TPos+cpB.TPos)*0.5;
        cp.LPos=avTPos;//note that setting LPos=TPos
        MAXMSG("MeuReco",logLevel,1000)
          <<"Calc LPos: pl="<<cp.Plane<<", st="<<cp.Strip
          <<", tposA="<<cpA.TPos
          <<", tposB="<<cpB.TPos
          <<", LPos="<<cp.LPos<<endl;
      }
    }
    else{
      MAXMSG("MeuReco",Msg::kWarning,1000)
        <<"LPos value already found="<<cp.LPos
        <<", pl="<<cp.Plane<<endl;
    }

    if (posDir) pl++;
    else pl--;
  }

  //now calc the ends
  //determine the planes to use to calc the lpos of the last plane
  Int_t plNearVtxA=s.VtxPlane-1;
  Int_t plNearVtxB=s.VtxPlane-3;
  if (posDir) {
    plNearVtxA=s.VtxPlane+1;
    plNearVtxB=s.VtxPlane+3;
  }
  Int_t plNearEndA=s.EndPlane+1;
  Int_t plNearEndB=s.EndPlane+3;
  if (posDir) {
    plNearEndA=s.EndPlane-1;
    plNearEndB=s.EndPlane-3;
  }
  
  Bool_t interpError=false;

  MeuHitInfo& cpEnd=plInfo[s.EndPlane];
  MeuHitInfo& cpEndA=plInfo[plNearEndA];
  MeuHitInfo& cpEndB=plInfo[plNearEndB];
  MeuHitInfo& cpVtx=plInfo[s.VtxPlane];
  MeuHitInfo& cpVtxA=plInfo[plNearVtxA];
  MeuHitInfo& cpVtxB=plInfo[plNearVtxB];

  //calc the lpos for the vtx and end planes
  for (Int_t i=0;i<2;i++){

    MeuHitInfo* cpGap=&cpVtx;
    MeuHitInfo* cp1=&cpVtxA;
    MeuHitInfo* cp2=&cpVtxB;
    if (i==1) {
      cpGap=&cpEnd;
      cp1=&cpEndA;
      cp2=&cpEndB;
    }

    Int_t pl=cpGap->Plane;

    //do interpolation
    //define m=x1_2/y1_2 = x1-x2/y1-y2  
    //(x and y are what is required)
    //define x=tpos, y=z
    Float_t m=(cp1->TPos-cp2->TPos)/(cp1->Z-cp2->Z);
    //Double_t x=(y-y2)*m+x2;// give above formula for m
    Float_t tpos=(cpGap->Z-cp2->Z)*m+cp2->TPos;
    if (cpGap) {
      //record result
      cpGap->LPos=tpos;//note that setting LPos=TPos
    }
    else {
      MAXMSG("MeuReco",Msg::kWarning,1000)
        <<"Ahhhhh, not found a cp for after gap"<<endl;
      interpError=true;
    }
    
    MAXMSG("MeuReco",logLevel,1000)
      <<"Doing interpolation:"<<endl
      <<"Point1:       tpos="<<cp1->TPos<<", z="<<cp1->Z<<endl
      <<"    gap sizes:    ="<<cp1->TPos-cpGap->LPos
      <<", z="<<cp1->Z-cpGap->Z<<endl
      <<"Interp point: tpos="<<cpGap->LPos
      <<", z="<<cpGap->Z<<endl
      <<"    gap sizes:    ="<<cpGap->LPos-cp2->TPos
      <<",  ="<<cpGap->Z-cp2->Z<<endl
      <<"Point2:       tpos="<<cp2->TPos<<", z="<<cp2->Z<<endl
      <<endl;
    
    VldTimeStamp vldts;
    //VldTimeStamp vldts(ev.UnixTime,0);
    VldContext vc
      (static_cast<Detector::Detector_t>(s.Detector),
       static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
    //get the ugh
    UgliGeomHandle ugh(vc);
    //Float_t zmin=-1;
    //Float_t zmax=-1;
    Float_t tmin=-1;
    Float_t tmax=-1;
    PlexPlaneId plidScint
      (static_cast<Detector::Detector_t>(s.Detector),pl);
    PlaneView::PlaneView_t pv=PlaneView::kUnknown;
    pv=ugh.GetScintPlnHandle(plidScint).GetPlaneView();
    string sPv="";//PlaneView::AsString(pv);
    
    if (cpGap->Strip<0) cout<<"Ahhhhh, bad strip="
                             <<cpGap->Strip<<endl; 
    PlexStripEndId seid
      (static_cast<Detector::Detector_t>(s.Detector),
       pl,cpGap->Strip,StripEnd::kWest);

    if (!seid.IsValid()){
      interpError=true;
      MSG("MeuReco",Msg::kInfo)
        <<"Strip not valid: pl="<<pl<<", st="<<cpGap->Strip
        <<", view="<<sPv<<endl;
      seid.Print();
      break;
    }

    UgliStripHandle ush=ugh.GetStripHandle(seid);
    Float_t  halflen=ush.GetHalfLength();
    TVector3 localCenter(0,0,0);//centre of the "logical" strip
    TVector3 localEastEnd(-halflen,0,0);
    TVector3 localWestEnd(+halflen,0,0);

    // calculate the positions in global XYZ space
    TVector3 xyzE = ush.LocalToGlobal(localEastEnd);
    TVector3 xyzW = ush.LocalToGlobal(localWestEnd);
    // calculate the positions in UVZ space
    TVector3 uvzE = ugh.xyz2uvz(xyzE);
    TVector3 uvzW = ugh.xyz2uvz(xyzW);
    //uvzE.Print();
    //uvzW.Print();

    Float_t lposE=-999;
    Float_t lposW=-999;
    if (pv==PlaneView::kU){
      lposE=uvzE.Y();//in u-view strip lies along v=y in vector
      lposW=uvzW.Y();//therefore to get strip length need y
    }
    else if (pv==PlaneView::kV){
      lposE=uvzE.X();//in v-view strip lies along u=x in vector
      lposW=uvzW.X();//therefore to get strip length need x
    }
    else cout<<"ahhhhhhh"<<endl;

    //work out min and max
    Float_t lmax=lposE;
    Float_t lmin=lposW;
    if (lposE<lposW){
      lmin=lposE;
      lmax=lposW;
    }

    //ugh.GetZExtent(zmin,zmax);
    ugh.GetTransverseExtent(pv,tmin,tmax);
    MAXMSG("MeuReco",Msg::kVerbose,1000)
      <<"pl="<<pl<<", st="<<cpGap->Strip<<", view="<<sPv<<endl
      //<<"Got ZExtent, min="<<zmin<<", max="<<zmax<<endl
      <<"Got TransverseExtent,   min="<<tmin<<", max="<<tmax<<endl
      <<"Got LongitudinalExtent, min="<<lmin<<", max="<<lmax<<endl;
      
    //check that the interpolation did not go crazy
    if (tpos>lmax || tpos<lmin){
      Float_t lposToUse=lmax;
      if (tpos<0) lposToUse=lmin;
      
      Float_t newLPos=tpos;
      
      //check difference
      //20 cm shift in 6 cm is not too steep
      if (fabs(lposToUse-tpos)<20*Munits::cm) {
        //just cap the tpos value at slightly below min/max
        newLPos=lposToUse;

        //now use the new tpos value
        if (cpGap) {
          //record result
          cpGap->LPos=newLPos;//note that setting LPos=TPos
        }
        else {
          MAXMSG("MeuReco",Msg::kWarning,1000)
            <<"Ahhhhh, not found a cp for after gap"<<endl;
          interpError=true;
        }
      }
      else{
        MAXMSG("MeuReco",logLevel,1000)
          <<endl<<"Interpolated out of the detector significantly"
          <<" for trk end plane"<<endl
          <<"LPos="<<tpos<<", not capping the value at "
          <<lposToUse<<endl
          <<"LongitudinalExtent: min="<<lmin<<", max="<<lmax<<endl;
        //don't cap it if it's this bad
        //newLPos=lposToUse;
        interpError=true;
      }
    }
  } 
  
  //last final check
  if (interpError){
    MAXMSG("MeuReco",logLevel,1000)
      <<"Couldn't fill LPos of track end planes"<<endl;
  }

  return !interpError;//return true if all ok
}

//......................................................................

void MeuReco::CalcStripDists
(const MeuSummary& s,std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  
  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"CalcStripDists: Get position along strip, track: "
    <<s.VtxPlane<<"->"<<s.EndPlane<<endl;

  VldTimeStamp vldts;
  //VldTimeStamp vldts(ev.UnixTime,0);
  VldContext vc
    (static_cast<Detector::Detector_t>(s.Detector),
     static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
  //get the ugh
  UgliGeomHandle ugh(vc);
  
  typedef map<Int_t,MeuHitInfo>::iterator cpIt;
  for (cpIt it=plInfo.begin();it!=plInfo.end();++it){

    MeuHitInfo& cp=it->second;
    Int_t pl=cp.Plane;
    
    PlexPlaneId plidScint
      (static_cast<Detector::Detector_t>(s.Detector),pl);
    PlaneView::PlaneView_t pv=PlaneView::kUnknown;
    pv=ugh.GetScintPlnHandle(plidScint).GetPlaneView();
    string sPv="";//PlaneView::AsString(pv);

    //work out the strip from u,v coord.
    Int_t strip=ugh.GetScintPlnHandle(plidScint).
      GetClosestStrip(cp.TPos,cp.LPos).GetSEId().GetStrip();
    //UgliStripHandle GetClosestStrip (tpos,orthCoord=9999.)
    if (cp.Strip>=0 && cp.Strip!=strip) {
      PlexStripEndId seidW
        (static_cast<Detector::Detector_t>(s.Detector),
         pl,strip,StripEnd::kWest);
      PlexStripEndId seidWOrig
        (static_cast<Detector::Detector_t>(s.Detector),
         pl,cp.Strip,StripEnd::kWest);
      
      MAXMSG("MeuReco",Msg::kInfo,5)
        <<"Weird strip: orig="<<cp.Strip
        <<", ugli="<<strip
        <<", ugli2="<<ugh.GetScintPlnHandle(plidScint).
        GetClosestStrip(cp.TPos).GetSEId().GetStrip()
        <<", pl="<<pl<<endl
        <<"tposOrig="<<ugh.GetStripHandle(seidWOrig).GetTPos()
        <<", tposNew="<<ugh.GetStripHandle(seidW).GetTPos()
        <<", tposNtp="<<cp.TPos<<endl;
      strip=cp.Strip;//use the original for now
    }
      
    //set the best guess strip
    cp.StripBestGuess=strip;

    PlexStripEndId seidW
      (static_cast<Detector::Detector_t>(s.Detector),
       pl,strip,StripEnd::kWest);

    MAXMSG("MeuReco",Msg::kDebug,100)
      <<"pl="<<pl<<", st="<<strip<<endl;
    
    if (!seidW.IsValid()){//only west because of ND
      MSG("MeuReco",Msg::kError)
        <<"Strip not valid: pl="<<pl<<", st="<<strip
        <<", view="<<sPv<<endl;
      seidW.Print();
      break;
    }
    
    //set some of the strip's variables
    UgliStripHandle ushW=ugh.GetStripHandle(seidW);
    cp.StripHalfLength=ushW.GetHalfLength();
    cp.StripClearFibre2=ushW.ClearFiber(StripEnd::kWest);
    cp.StripPigtail2=ushW.WlsPigtail(StripEnd::kWest);

    if (s.Detector!=Detector::kNear){
      PlexStripEndId seidE
        (static_cast<Detector::Detector_t>(s.Detector),
         pl,strip,StripEnd::kEast);
      UgliStripHandle ushE=ugh.GetStripHandle(seidE);
      cp.StripClearFibre1=ushE.ClearFiber(StripEnd::kEast);
      cp.StripPigtail1=ushE.WlsPigtail(StripEnd::kEast);
    }

    //get the distance to the strip centre
    TVector3 xyz(cp.X,cp.Y,cp.Z);
    //static Int_t i=0;
    //++i;
    //if (i<200) xyz.Print();
    TVector3 localPos=ushW.GlobalToLocal(xyz);
    //if (i<200) localPos.Print();
    cp.DistToStripCentre=localPos.X();
    
    //TVector3 localCenter(0,0,0);//centre of the "logical" strip
    //TVector3 localEastEnd(-halflen,0,0);
    //TVector3 localWestEnd(+halflen,0,0);

    // calculate the positions in global XYZ space
    //TVector3 xyzE = ush.LocalToGlobal(localEastEnd);
    //TVector3 xyzW = ush.LocalToGlobal(localWestEnd);
    // calculate the positions in UVZ space
    //TVector3 uvzE = ugh.xyz2uvz(xyzE);
    //TVector3 uvzW = ugh.xyz2uvz(xyzW);
    //uvzE.Print();
    //uvzW.Print();
  } 
}

//......................................................................

void MeuReco::ConvertToXYZ(MeuSummary& s,
                            std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  
  MAXMSG("MeuReco",logLevel,1000)
    <<endl<<"ConvertToXYZ: Turn TPos, LPos and z -> xyz"<<endl;

  //Float_t initValue=-999;
  Bool_t posDir=(s.EndPlane>s.VtxPlane);
  
  //get variables for loops
  Int_t pl=s.VtxPlane;
  Int_t tmpEndPlane=s.EndPlane;
  if (posDir) tmpEndPlane++;//make sure last plane is used
  else tmpEndPlane--;
  
  VldTimeStamp vldts;
  //VldTimeStamp vldts(ev.UnixTime,0);
  VldContext vc
    (static_cast<Detector::Detector_t>(s.Detector),
     static_cast<SimFlag::ESimFlag>(s.SimFlag),vldts);
  //get the ugh
  UgliGeomHandle ugh(vc);
  TVector3 uvz;

  //loop and calc xyz using tpos, lpos and z in each plane
  while (pl!=tmpEndPlane){
    MAXMSG("MeuReco",Msg::kVerbose,1000)
      <<"pl="<<pl<<endl;
    
    //cache the current cp
    MeuHitInfo& cp=plInfo[pl];

    //set z first
    uvz.SetZ(cp.Z);
    
    //now decide what tpos and lpos measure depending on the view
    if (cp.View==PlaneView::kU){
      uvz.SetX(cp.TPos);//in u-view tpos measures u
      uvz.SetY(cp.LPos);//in u-view lpos measures v
    }
    else if (cp.View==PlaneView::kV){
      uvz.SetY(cp.TPos);//in v-view tpos measures v
      uvz.SetX(cp.LPos);//in v-view lpos measures u
    }
    else cout<<"ahhhhhhh"<<endl;

    //now calc xyz
    TVector3 xyz=ugh.uvz2xyz(uvz);
    
    //now set the cp x and y
    cp.X=xyz.X();
    cp.Y=xyz.Y();

    MAXMSG("MeuReco",logLevel,1000)
      <<"pl="<<pl<<", cp: X="<<cp.X<<", Y="<<cp.Y<<", Z="<<cp.Z
      <<endl;
      
    if (posDir) pl++;
    else pl--;
  } 
}

//......................................................................

Bool_t MeuReco::Reconstruct(MeuSummary& s,
                             std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
   
  MSG("MeuReco",logLevel)
    <<"Running Reconstruct..."<<endl;

  //steps:
  //determine lpos for all planes
  //determine the x,y,z 

  //this function should only accept tracks which finish
  //in similar places in two views
    
  Bool_t allGapsFilled=this->CalcPosOfGaps(s,plInfo);
  if (!allGapsFilled) {
    MSG("MeuReco",logLevel)
      <<"Running CalcPosOfBigGaps..."<<endl;
    allGapsFilled=this->CalcPosOfBigGaps(s,plInfo);
  }
  if (!allGapsFilled) {
    MSG("MeuReco",logLevel)
      <<"Running CalcPosOfTrkEndGaps..."<<endl;
    allGapsFilled=this->CalcPosOfTrkEndGaps(s,plInfo);
  }
  if (!allGapsFilled) {
    MAXMSG("MeuReco",Msg::kDebug,100)
      <<"Giving up on this event, gaps not filled"<<endl;
    //this->PrintMeuHitInfo(plInfo);
    return false;//have to give up now...
  }
    
  MSG("MeuReco",logLevel)
    <<"Running CalcLPos..."<<endl;
    
  Bool_t goodLPos=this->CalcLPos(s,plInfo);
  if (!goodLPos) {
    MAXMSG("MeuReco",Msg::kDebug,100)
      <<"Giving up on this event, LPos is bad"<<endl;
    return false;//have to give up now...
  }

  MSG("MeuReco",logLevel)
    <<"Running ConvertToXYZ..."<<endl;
  //calcuate xyz from tpos and lpos
  this->ConvertToXYZ(s,plInfo);
    
  MSG("MeuReco",logLevel)
    <<"Running CalcPLCor..."<<endl;
  //do the fit to determine the PLCor
  this->CalcPLCor(s,plInfo);//needs xyz  

  return true;
}

//......................................................................

void MeuReco::CalcWindow(MeuSummary& s,
                          std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
    
  MSG("MeuReco",logLevel)
    <<"Running CalcWindow..."<<endl;
  //MC
  //this->TrueFirstLastPl(s.MCVtxPlane,s.MCEndPlane,
  //            s.MCTrueWinVtxSidePl,s.MCTrueWinStopSidePl);

  //this->ExtractTrueLPos(plInfo,s);

  //this->ExtractTrueEnDep(plInfo,s);

  MSG("MeuReco",logLevel)
    <<"Running CalcMEU..."<<endl;
  //actually get the MEU value!!!
  this->CalcMEU(s,plInfo);

  if (s.WinSigMap>0){
    static Float_t avMeuSigCor=0;
    static Float_t avMeuSigMap=0;
    static Int_t counter=0;
    avMeuSigCor+=s.WinSigCor;
    avMeuSigMap+=s.WinSigMap;
    counter++;
      
    MAXMSG("MeuReco",Msg::kDebug,1000)
      <<"Completed the window:"<<endl
      <<"  s.WinSigMap="<<s.WinSigMap
      <<", runing av="<<avMeuSigMap/counter<<endl
      <<"  s.WinSigCor="<<s.WinSigCor
      <<", runing av="<<avMeuSigCor/counter
      <<" (N="<<counter<<")"<<endl;
  }
}

//......................................................................

void MeuReco::CalcVtxSpecialND
(MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //this function allows muons to start in the spectrometer
  //but it converts the vtx plane to be 120 to get round 
  //the wide active plane spacing in the ND

  if (s.Detector!=Detector::kNear) return;

  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;

  if (s.VtxPlane<=120 && s.EndPlane<=120) {
    //do nothing
  }
  else if (s.VtxPlane>120 && s.EndPlane<=118) {
    for (Int_t pl=120;pl>0;pl--){
      if (plInfo.find(pl)!=plInfo.end()){
        const MeuHitInfo& cp=plInfo.find(pl)->second;
        MSG("MeuReco",logLevel)
          <<"MeuReco::CalcVtxSpecialND pl="<<pl
          <<", st="<<cp.Strip<<endl;
        if (cp.Strip!=-1) {
          s.VtxPlane=pl;
          MSG("MeuReco",logLevel)
            <<"Breaking out of loop, s.VtxPlane="
            <<s.VtxPlane<<endl;
          break;
        }
      }
      else {
        MSG("MeuReco",logLevel)
          <<"CalcVtxSpecialND: No pl="<<pl<<endl;
      }
    }
  }
  else if (s.EndPlane>120) {
    MSG("MeuReco",Msg::kError)
      <<"End plane is in Spectrometer"<<endl;
  }
  else MSG("MeuReco",Msg::kError)<<"Spectrometer issues!!!"<<endl;

  //last final check
  if (s.VtxPlane>120) {
    MSG("MeuReco",Msg::kError)
      <<"Ahhh, last part of special missed vtx plane >120, pl="
      <<s.VtxPlane<<endl;
    this->PrintMeuHitInfo(plInfo);
  }
}

//......................................................................

void MeuReco::CalcVtx(MeuSummary& s,
                       const std::map<Int_t,MeuHitInfo>& plInfo) const
{
  //vtx plane is the one with the greatest Y for downward events
  //for cosmics
  
  map<Int_t,MeuHitInfo>::const_iterator itBegin=plInfo.begin();
  map<Int_t,MeuHitInfo>::const_iterator itEnd=--plInfo.end();
  
  //loop and find the first tracked plane
  while (itBegin->second.Strip<0 || itEnd->second.Strip<0){
    
    MAXMSG("MeuReco",Msg::kDebug,500)
      <<"Beg/End plane not on track:"<<endl
      <<"Beg: pl="<<itBegin->second.Plane
      <<", st="<<itBegin->second.Strip<<endl
      <<", End: pl="<<itEnd->second.Plane
      <<", st="<<itEnd->second.Strip
      <<endl;
  
    //iterate the required iterator
    if (itBegin->second.Strip<0) ++itBegin;
    if (itEnd->second.Strip<0) --itEnd;

    MAXMSG("MeuReco",Msg::kDebug,500)
      <<"trying next plane in on track..."<<endl
      <<"Beg: pl="<<itBegin->second.Plane
      <<", st="<<itBegin->second.Strip<<endl
      <<", End: pl="<<itEnd->second.Plane
      <<", st="<<itEnd->second.Strip
      <<endl;

    //check that haven't iterated off the end
    if (itBegin==plInfo.end()) {
      MAXMSG("MeuReco",Msg::kWarning,500)
        <<"Failed to find beg/end strip"<<endl;
      break;
    }

  }

  if (s.TrigSrc>60000) {//spill mode
    MAXMSG("MeuReco",Msg::kInfo,1)
      <<"CalcVtx in spill mode"<<endl;
    s.VtxPlane=itBegin->second.Plane;
    s.EndPlane=itEnd->second.Plane;
    
    const MeuHitInfo& cvtx=itBegin->second;
    const MeuHitInfo& cend=itEnd->second;
    if (s.Detector==Detector::kNear){
      s.RVtx=sqrt(pow(cvtx.X-1.5,2)+pow(cvtx.Y,2));
      s.REnd=sqrt(pow(cend.X-1.5,2)+pow(cend.Y,2));
    }
    else if (s.Detector==Detector::kFar){
      s.RVtx=sqrt(pow(cvtx.X,2)+pow(cvtx.Y,2));
      s.REnd=sqrt(pow(cend.X,2)+pow(cend.Y,2));
    }
    else cout<<"No detector type recognised"<<endl;
  }
  else {//dynode triggers i.e. cosmics (probably)
    MAXMSG("MeuReco",Msg::kInfo,1)
      <<"CalcVtx in cosmics mode"<<endl;
    //now actually work out which end is the vtx
    if (itBegin->second.Y>itEnd->second.Y) {
      s.VtxPlane=itBegin->second.Plane;
      s.EndPlane=itEnd->second.Plane;
        
      const MeuHitInfo& cvtx=itBegin->second;
      const MeuHitInfo& cend=itEnd->second;
      if (s.Detector==Detector::kNear) {
        s.RVtx=sqrt(pow(cvtx.X-1.5,2)+pow(cvtx.Y,2));
        s.REnd=sqrt(pow(cend.X-1.5,2)+pow(cend.Y,2));
      }
      else if (s.Detector==Detector::kFar) {
        s.RVtx=sqrt(pow(cvtx.X,2)+pow(cvtx.Y,2));
        s.REnd=sqrt(pow(cend.X,2)+pow(cend.Y,2));
      }
      else cout<<"No detector type recognised"<<endl;
    }
    else {
      s.VtxPlane=itEnd->second.Plane;
      s.EndPlane=itBegin->second.Plane;
        
      const MeuHitInfo& cvtx=itEnd->second;
      const MeuHitInfo& cend=itBegin->second;
      if (s.Detector==Detector::kNear) {
        s.RVtx=sqrt(pow(cvtx.X-1.5,2)+pow(cvtx.Y,2));
        s.REnd=sqrt(pow(cend.X-1.5,2)+pow(cend.Y,2));
      }
      else if (s.Detector==Detector::kFar) {
        s.RVtx=sqrt(pow(cvtx.X,2)+pow(cvtx.Y,2));
        s.REnd=sqrt(pow(cend.X,2)+pow(cend.Y,2));
      }
      else cout<<"No detector type recognised"<<endl;
    }
  }
}

//......................................................................

Float_t MeuReco::CalcSmallestDeepDistToEdge(const MeuHitInfo& cp) const
{
  //Deep containment is defined as being within the overlap of 
  //all different plane coverages

  Float_t distUPI=0;
  Float_t distVPI=0;
  Float_t distUFI=0;
  Float_t distVFI=0;
  Float_t xedge=0;
  Float_t yedge=0;

  PlaneOutline po;

  //calc for U-view partial
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kU,
                    PlaneCoverage::kNearPartial,
                    distUPI,xedge,yedge);
  Bool_t isInsideUPI=po.IsInside(cp.X,cp.Y,PlaneView::kU,
                                 PlaneCoverage::kNearPartial);

  //calc for U-view full
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kU,
                    PlaneCoverage::kNearFull,
                    distUFI,xedge,yedge);
  Bool_t isInsideUFI=po.IsInside(cp.X,cp.Y,PlaneView::kU,
                                 PlaneCoverage::kNearFull);

  //calc for V-view partial
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kV,
                    PlaneCoverage::kNearPartial,
                    distVPI,xedge,yedge);
  Bool_t isInsideVPI=po.IsInside(cp.X,cp.Y,PlaneView::kV,
                                 PlaneCoverage::kNearPartial);

  //calc for V-view full
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kV,
                    PlaneCoverage::kNearFull,
                    distVFI,xedge,yedge);
  Bool_t isInsideVFI=po.IsInside(cp.X,cp.Y,PlaneView::kV,
                                 PlaneCoverage::kNearFull);

  //check that hit has "deep" containment and find 
  //the closest distance to the edge of the combined overlap region
  Float_t smallestDist=999;
  if (isInsideUPI && isInsideUFI && isInsideVPI && isInsideVFI){
    if (smallestDist>distUPI) smallestDist=distUPI;
    if (smallestDist>distUFI) smallestDist=distUFI;
    if (smallestDist>distVPI) smallestDist=distVPI;
    if (smallestDist>distVFI) smallestDist=distVFI;
  }

  //if not deeply contained then find the closest distance
  //to the combined overlap region
  //this is actually the furthest distance!
  if (smallestDist==999) {
    smallestDist=0;
    if (smallestDist<distUPI && !isInsideUPI) smallestDist=distUPI;
    if (smallestDist<distUFI && !isInsideUFI) smallestDist=distUFI;
    if (smallestDist<distVPI && !isInsideVPI) smallestDist=distVPI;
    if (smallestDist<distVFI && !isInsideVFI) smallestDist=distVFI;

    //make the distance negative for hits outside deep containment
    smallestDist*=-1;
  }

  if (smallestDist==999) cout<<"Ahhhhhhh"<<endl;

  MAXMSG("MeuReco",Msg::kDebug,500)
    <<"smallestDist="<<smallestDist
    <<", UPI="<<distUPI<<" (in="<<isInsideUPI<<")"
    <<", UFI="<<distUFI<<" (in="<<isInsideUFI<<")"
    <<", VPI="<<distVPI<<" (in="<<isInsideVPI<<")"
    <<", VFI="<<distVFI<<" (in="<<isInsideVFI<<")"<<endl;
    
  return smallestDist;
}

//......................................................................

Float_t MeuReco::CheckDeepDistToEdge(const MeuHitInfo& cp) const
{
  //Deep containment is defined as being within the overlap of 
  //all different plane coverages

  Float_t distUPI=0;
  Float_t distVPI=0;
  Float_t distUFI=0;
  Float_t distVFI=0;
  Float_t xedge=0;
  Float_t yedge=0;

  PlaneOutline po;

  //calc for U-view partial
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kU,
                    PlaneCoverage::kNearPartial,
                    distUPI,xedge,yedge);
  Bool_t isInsideUPI=po.IsInside(cp.X,cp.Y,PlaneView::kU,
                                 PlaneCoverage::kNearPartial);

  //calc for U-view full
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kU,
                    PlaneCoverage::kNearFull,
                    distUFI,xedge,yedge);
  Bool_t isInsideUFI=po.IsInside(cp.X,cp.Y,PlaneView::kU,
                                 PlaneCoverage::kNearFull);

  //calc for V-view partial
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kV,
                    PlaneCoverage::kNearPartial,
                    distVPI,xedge,yedge);
  Bool_t isInsideVPI=po.IsInside(cp.X,cp.Y,PlaneView::kV,
                                 PlaneCoverage::kNearPartial);

  //calc for V-view full
  po.DistanceToEdge(cp.X,cp.Y,PlaneView::kV,
                    PlaneCoverage::kNearFull,
                    distVFI,xedge,yedge);
  Bool_t isInsideVFI=po.IsInside(cp.X,cp.Y,PlaneView::kV,
                                 PlaneCoverage::kNearFull);

  //check that hit has "deep" containment
  Float_t smallestDist=999;
  if (isInsideUPI && isInsideUFI && isInsideVPI && isInsideVFI){
    if (smallestDist>distUPI) smallestDist=distUPI;
    if (smallestDist>distUFI) smallestDist=distUFI;
    if (smallestDist>distVPI) smallestDist=distVPI;
    if (smallestDist>distVFI) smallestDist=distVFI;
  }

  //set to be right on the edge of the detector
  if (smallestDist==999) smallestDist=0;

  MAXMSG("MeuReco",Msg::kDebug,500)
    <<"smallestDist="<<smallestDist
    <<", UPI="<<distUPI
    <<", UFI="<<distUFI
    <<", VPI="<<distVPI
    <<", VFI="<<distVFI<<endl;
    
  return smallestDist;
}

//......................................................................

void MeuReco::CalcFCPC(MeuSummary& s,
                        const std::map<Int_t,MeuHitInfo>& plInfo) const
{

  if (s.Detector==Detector::kNear) {   
    //set these just to be explicit
    s.SM1=true;
    s.SM2=false;
    s.SMBoth=false;

    //check if muon entered SM end
    s.EntSM1Front=s.VtxPlane<=5;//first 5 planes
    s.EntSM1Back=s.VtxPlane>=116;//120,119,18,17,16
    s.EntSMEnd=s.EntSM1Front || s.EntSM1Back;
    
    //check if muon exited SM end
    s.ExitSM1Front=s.EndPlane<=5;//first 5 planes
    s.ExitSM1Back=s.EndPlane>=116;//120,119,18,17,16
    s.ExitSMEnd=s.ExitSM1Front || s.ExitSM1Back;
    
    //check containment
    //vtx plane
    if (plInfo.find(s.VtxPlane)==plInfo.end()) {
      MSG("MeuReco",Msg::kFatal)<<"Ahhhhh! No vtx"<<endl;
    }
    const MeuHitInfo& cp=plInfo.find(s.VtxPlane)->second;
    s.VtxDistToEdge=this->CheckDeepDistToEdge(cp);
    
    //end plane
    if (plInfo.find(s.EndPlane)==plInfo.end()) {
      MSG("MeuReco",Msg::kError)<<"Ahhhhh2! No end"<<endl;
    }
    const MeuHitInfo& cpEnd=plInfo.find(s.EndPlane)->second;
    s.EndDistToEdge=this->CheckDeepDistToEdge(cpEnd);
    
    //decide if FC or PC
    if (s.TrigSrc>60000) {//spill mode
      MAXMSG("MeuReco",Msg::kInfo,1)
        <<"CalcFCPC: Spill mode"<<endl;
      s.PC=!s.ExitSMEnd && s.EndDistToEdge>s.DistToEdgeFid;
    }
    else {//cosmics mode
      MAXMSG("MeuReco",Msg::kInfo,1)
        <<"CalcFCPC: Cosmics mode"<<endl;
      s.PC=!s.ExitSMEnd && s.EndDistToEdge>s.DistToEdgeFid && 
        (s.VtxDistToEdge<s.DistToEdgeFid || s.EntSMEnd);
    }

    if (!s.PC) {
      s.FC=!s.ExitSMEnd && !s.EntSMEnd && 
        s.VtxDistToEdge>s.DistToEdgeFid && 
        s.EndDistToEdge>s.DistToEdgeFid;
    }
    else s.FC=false;//don't want something both FC and PC
    
    //s.PC=!s.ExitSMEnd && s.REnd<s.RFid && 
    //(s.RVtx>s.RFid || s.EntSMEnd);
    //s.FC=!s.ExitSMEnd && !s.EntSMEnd && 
    //s.RVtx<s.RFid && s.REnd<s.RFid;
    //if (s.PC){
    //this->PrintContainment(s);
    //}
  }
  else if (s.Detector==Detector::kFar) {

    Float_t maxPlane=s.MaxPlane2;
    if (s.MaxPlane3>s.MaxPlane2) maxPlane=s.MaxPlane3;
    Float_t minPlane=s.MinPlane2;
    if (s.MinPlane3<s.MinPlane2) minPlane=s.MinPlane3;

    //work out which SM or if in both
    s.SM1=maxPlane<=248; 
    s.SM2=minPlane>=250;
    s.SMBoth=minPlane<=248 && maxPlane>=250;
    
    //record the min and max planes in each view
    s.MinPlane2=s.MinPlane2;  
    s.MinPlane3=s.MinPlane3;
    s.MaxPlane2=s.MaxPlane2;  
    s.MaxPlane3=s.MaxPlane3;
    
    //check if muon entered SM end
    s.EntSM1Front=s.VtxPlane<=5;//first 5 planes
    s.EntSM1Back=s.VtxPlane>=244 && s.VtxPlane<=248;//248,7,6,5,4
    s.EntSM2Front=s.VtxPlane>=250 && s.VtxPlane<=254;//250,1,2,3,4
    s.EntSM2Back=s.VtxPlane>=481;//485,4,3,2,1
    s.EntSMEnd=s.EntSM1Front || s.EntSM1Back || s.EntSM2Front || 
      s.EntSM2Back;
    
    //check if muon exited SM end
    s.ExitSM1Front=s.EndPlane<=5;//first 5 planes
    s.ExitSM1Back=s.EndPlane>=244 && s.EndPlane<=248;//248,7,6,5,4
    s.ExitSM2Front=s.EndPlane>=250 && s.EndPlane<=254;//250,1,2,3,4
    s.ExitSM2Back=s.EndPlane>=481;//485,4,3,2,1
    s.ExitSMEnd=s.ExitSM1Front || s.ExitSM1Back ||
      s.ExitSM2Front || s.ExitSM2Back;
    
    //work out if FC or PC
    s.PC=!s.ExitSMEnd && s.REnd<s.RFid &&
      (s.RVtx>s.RFid || s.EntSMEnd);
    s.FC=!s.ExitSMEnd && !s.EntSMEnd && s.RVtx<s.RFid &&
      s.REnd<s.RFid;

    //this->PrintContainment(s);
  }
  else {
    MAXMSG("MeuReco",Msg::kFatal,500)
      <<"CalcFCPC: Detector not implemented = "<<s.Detector<<endl;
  }

  
  Bool_t print=false;
  if (print){
    if (s.PC){
      MAXMSG("MeuReco",Msg::kInfo,500)
        <<"PC"
        <<", vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
        <<", dVtx="<<s.VtxDistToEdge
        <<", dEnd="<<s.EndDistToEdge<<endl;
      for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.begin();
           it!=plInfo.end();++it){
        const MeuHitInfo& c=it->second;
        MAXMSG("MeuReco",Msg::kInfo,500)
          <<"pl="<<c.Plane<<", y="<<c.Y<<endl;
      }
    }
    else if (s.FC){
      MAXMSG("MeuReco",Msg::kInfo,500)
        <<"FC!!!!!"
        <<", vtx Pl="<<s.VtxPlane<<", end Pl="<<s.EndPlane
        <<", dVtx="<<s.VtxDistToEdge
        <<", dEnd="<<s.EndDistToEdge<<endl;
      for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.begin();
           it!=plInfo.end();++it){
        const MeuHitInfo& c=it->second;
        MAXMSG("MeuReco",Msg::kInfo,500)
          <<"pl="<<c.Plane<<", y="<<c.Y<<endl;
      }
    }
    else {
      MAXMSG("MeuReco",Msg::kInfo,500)
        <<"NOT PC or FC"<<endl;
    }
  }
}

//......................................................................

void MeuReco::PrintContainment(const MeuSummary& s) const
{
  MSG("MeuReco",Msg::kInfo)
    <<"Containment summary:"<<endl;
  MSG("MeuReco",Msg::kInfo)
    <<"  Track: "<<s.VtxPlane<<" -> "<<s.EndPlane<<endl;
  MSG("MeuReco",Msg::kInfo)
    <<"  EntSMEnd="<<s.EntSMEnd
    <<" (SM1: F="<<s.EntSM1Front<<", B="<<s.EntSM1Back
    <<" SM2: F="<<s.EntSM2Front<<", B="<<s.EntSM2Back<<")"<<endl;
  MSG("MeuReco",Msg::kInfo)
    <<"  ExitSMEnd="<<s.ExitSMEnd
    <<" (SM1: F="<<s.ExitSM1Front<<", B="<<s.ExitSM1Back
    <<" SM2: F="<<s.ExitSM2Front<<", B="<<s.ExitSM2Back<<")"<<endl;
  MSG("MeuReco",Msg::kInfo)
    <<"  Radius: fREnd="<<s.REnd<<", fRVtx="<<s.RVtx
    <<", (fRFid="<<s.RFid<<")"<<endl;
  MSG("MeuReco",Msg::kInfo)
    <<"  VtxDistToEdge="<<s.VtxDistToEdge
    <<", EndDistToEdge="<<s.EndDistToEdge
    <<", PC="<<s.PC<<", FC="<<s.FC<<endl;
}

//......................................................................

Bool_t MeuReco::CheckWinContainment
(const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo) const
{
  Bool_t contained=true;
  Float_t containDisttoEdge=0.1;//fiducial distance from edge for win

  if (s.Detector==Detector::kNear) {
    //perform check that every plane in window is in fiducial vol
    Int_t lowPl=s.WinVtxSidePl;
    Int_t highPl=s.WinStopSidePl;
    if (lowPl>highPl){
      lowPl=s.WinStopSidePl;
      highPl=s.WinVtxSidePl;
    }
    
    //check that planes exist
    if (plInfo.find(lowPl)==plInfo.end() ||
        plInfo.find(highPl)==plInfo.end()) {
      MAXMSG("MeuReco",Msg::kFatal,200)
        <<"Can't find track window edge planes, returning !contained"
        <<endl;
      contained=false;
      return contained;
    }
    
    //loop over planes in window
    for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.find(lowPl);
         it!=++plInfo.find(highPl);++it){//stop at 1 past highPl
      
      const MeuHitInfo& cp=it->second;
      
      Float_t distToEdge=this->CheckDeepDistToEdge(cp);
      
      if(distToEdge<containDisttoEdge) {
        contained=false;
        MAXMSG("MeuReco",Msg::kDebug,500)
          <<"pl="<<cp.Plane<<", st="<<cp.Strip
          <<", d="<<distToEdge
          <<", con="<<contained<<endl;
        break;
      }
    }
  }
  else if (s.Detector==Detector::kFar) {
    //do nothing
  }
  else {
    cout<<"Detector type not supported"<<endl;
  }
  
  return contained;
}

//......................................................................

Bool_t MeuReco::CheckTrackGaps
(const std::map<Int_t,MeuHitInfo>& plInfo) const
{
  Bool_t gapsOk=true;
  Int_t lastTrkPl=-1;
  Int_t gap=0;

  //strips TPos:
  //plane 6 goes -2.64 -> 1.32 m
  //plane 11 goes -1.32 -> 2.64 m
    
  //plane 4,8,10 goes -2.40 -> 0.24 m (V-view)
  //plane 5,7,9 goes -0.24 -> 2.40 m (U-view)
    
  //2.64 - 2.40 = 24 cm = 5.85 strips
  //the FI planes "stick out" by ~6 strips

  //Tobi's code snippet:
  // in the near detector, a further check is needed:
  // partial U planes have strips 0-63
  // partial V planes have strips 4-67
  //if (det==static_cast<Detector::Detector_t>(s.Detector)) {
  //  if (((pl-1)%5) && (pl%2)    && st>63) continue;
  //  if (((pl-1)%5) && (pl%2)==0 && st<4 ) continue;
  //}

  //define the fiducial volume
  Float_t minU=-0.14*Munits::m;
  Float_t maxU=2.3*Munits::m;
  Float_t minV=-maxU;
  Float_t maxV=-minU;
  MAXMSG("MeuReco",Msg::kInfo,50)
    <<"minU="<<minU<<", maxU="<<maxU
    <<", minV="<<minV<<", maxV="<<maxV
    <<endl;

  for (map<Int_t,MeuHitInfo>::const_iterator it=plInfo.begin();
       it!=plInfo.end();++it){
    const MeuHitInfo& cp=it->second;
    Int_t pl=cp.Plane;
    if (pl>120) continue;

    Float_t u=cp.TPos;
    Float_t v=cp.LPos;
    if (cp.View==PlaneView::kV){
      v=cp.TPos;//in the V-view tpos measures V
      u=cp.LPos;
    }
        
    Bool_t fid=false;
    if (u>minU && u<maxU && v>minV && v<maxV){
      fid=true;
    }
        
    if (cp.Strip!=-1) {
      gap=0;
      lastTrkPl=pl;
    }
    else gap++;

    MAXMSG("MeuReco",Msg::kInfo,500)
      <<"pl="<<cp.Plane<<", st="<<cp.Strip
      <<", sigcor="<<cp.SigCor
      <<", gap="<<gap
      <<", fid="<<fid<<", u="<<u<<", v="<<v<<endl;

    if (gap>=3) gapsOk=false;

  }

  return gapsOk;
}

//......................................................................

Bool_t MeuReco::CalcTrace
(const MeuSummary& s,const std::map<Int_t,MeuHitInfo>& plInfo,
 Float_t* trace,Float_t* traceZ) const
{
  if (s.Detector!=Detector::kFar) {
    MAXMSG("MeuReco",Msg::kInfo,1)
      <<"Not calculating trace for detector="<<s.Detector
      <<", only implemented for FD"<<endl;
    return true;
  }

  //variable to turn on all the useful messages if required
  Msg::LogLevel_t logLevel=Msg::kDebug;
  
  MAXMSG("MeuReco",logLevel,50)
    <<endl<<"CalcTrace: Get position along strip, track: "
    <<s.VtxPlane<<"->"<<s.EndPlane<<endl;

  Bool_t posDir=(s.EndPlane>s.VtxPlane);

  Int_t pl=s.EndPlane;
  if (posDir) pl-=5;//step back 5 if pos direction
  else pl+=5;

  const MeuHitInfo& cpEnd=plInfo.find(s.EndPlane)->second;
  const MeuHitInfo* tmpcp=&(plInfo.find(pl)->second);
  
  if (tmpcp->Strip<0){
    MAXMSG("MeuReco",Msg::kDebug,100)
      <<"Can't find strip for this plane="<<pl
      <<", strip="<<tmpcp->Strip
      <<", trk: "<<s.VtxPlane<<"->"<<s.EndPlane<<endl;

    if (posDir) pl--;//step back another plane if pos direction
    else pl++;
    tmpcp=&(plInfo.find(pl)->second);

    if (tmpcp->Strip<0){
      MAXMSG("MeuReco",Msg::kDebug,1000)
        <<"Can't find strip (again) for this plane="<<pl
        <<", strip="<<tmpcp->Strip
        <<", trk: "<<s.VtxPlane<<"->"<<s.EndPlane<<endl;

      if (posDir) pl--;//step back another plane if pos direction
      else pl++;
      tmpcp=&(plInfo.find(pl)->second);
      
      if (tmpcp->Strip<0){
        MAXMSG("MeuReco",logLevel,100)
          <<"Can't find strip (again2) for this plane="<<pl
          <<", strip="<<tmpcp->Strip
          <<", trk: "<<s.VtxPlane<<"->"<<s.EndPlane<<endl
          <<"Given up looking... and returned false"<<endl;
        return false;
      }
    }
    MAXMSG("MeuReco",Msg::kDebug,10)
      <<"Found pl="<<pl<<", strip="<<tmpcp->Strip<<endl;
  }
  
  //get a reference to the calibpos pl5
  const MeuHitInfo& cp=*tmpcp;

  static vector<TVector3> va;
  static vector<TVector3> vb;
  static vector<TVector3> vc;

  //initialise all the vectors if they are not already
  if (va.size()==0){
    //get the points on the planes    

    //get the points on the 8 detector sides (of the octagon)
    va.push_back(TVector3( 0, 4   ,0));//top 0 or 12 o'clock
    va.push_back(TVector3(-4, 1.66,0));
    va.push_back(TVector3(-4, 0   ,0));//3 o'clock
    va.push_back(TVector3(-4,-1.66,0));
    va.push_back(TVector3( 0,-4   ,0));//bottom or 6 o'clock
    va.push_back(TVector3( 4,-1.66,0));
    va.push_back(TVector3( 4, 0   ,0));//9 o'clock
    va.push_back(TVector3( 4, 1.66,0));

    //get the points on the 4 super-module ends
    va.push_back(TVector3( 0, 0   , 0));
    va.push_back(TVector3( 0, 0   ,14.772));
    va.push_back(TVector3( 0, 0   ,15.957));
    va.push_back(TVector3( 0, 0   ,29.948));

    //get the 1st vector

    //get the vectors on the 8 detector sides (of the octagon)
    //going round clockwise
    vb.push_back(TVector3(-1, 0, 0));
    vb.push_back(TVector3(-1,-1, 0));//45 degrees
    vb.push_back(TVector3( 0,-1, 0));
    vb.push_back(TVector3( 1,-1, 0));//45 degrees

    vb.push_back(TVector3( 1, 0, 0));
    vb.push_back(TVector3( 1, 1, 0));//45 degrees
    vb.push_back(TVector3( 0, 1, 0));
    vb.push_back(TVector3(-1, 1, 0));//45 degrees

    //get the vectors of the planes on the 4 super-module ends
    for (Int_t i=0;i<4;i++){
      vb.push_back(TVector3( 0, 1, 0));
    }

    //get the 2nd vector
    
    //get the vectors along the 8 detector sides (of the octagon)
    for (Int_t i=0;i<8;i++){
      vc.push_back(TVector3( 0, 0, 1));
    }
    //get the vectors of the planes on the 4 super-module ends
    for (Int_t i=0;i<4;i++){
      vc.push_back(TVector3( 1, 0, 0));
    }
  }

  MAXMSG("MeuReco",logLevel,1)
    <<"Sizes: va="<<va.size()
    <<", vb="<<vb.size()
    <<", vc="<<vc.size()<<endl;

  const TVector3 d(cpEnd.X,cpEnd.Y,cpEnd.Z);
  Float_t dX=cpEnd.X-cp.X;
  Float_t dY=cpEnd.Y-cp.Y;

  //make sanity checks on the directions
  //to ensure distance to all planes can be found
  if (dX==0){//will never hit some edges
    dX+=0.001;//add on 1mm so it just hits
    MAXMSG("MeuReco",Msg::kDebug,1000)
      <<"Tweaking dX so !=0"<<endl;
  }
  if (dY==0){//will never hit some edges
    dY-=0.001;//subtract 1mm so it just hits
    MAXMSG("MeuReco",Msg::kDebug,1000)
      <<"Tweaking dY so !=0"<<endl;
  }
  //the denominator below is zero if dX and dY are equal
  //so just tweak them a bit if so.
  if (fabs(dX)==fabs(dY)){
    dX+=0.001;//add on 1mm just to make them not the same
    MAXMSG("MeuReco",Msg::kDebug,1000)
      <<"Tweaking dX so !=dY"<<endl;
  }
  const TVector3 e(dX,dY,cpEnd.Z-cp.Z);//end-pl5

  
  vector<Float_t> dist(va.size(),-999);//hold the dist. to intersection 
  vector<TVector3> vintersect(va.size());//hold the intersection points
  Float_t closestDist=999999;
  Int_t closesti=-1;
  Float_t closestdZ=999;

  MAXMSG("MeuReco",logLevel,50)
    <<"pl="<<pl<<", strip="<<cp.Strip
    <<", trk: "<<s.VtxPlane<<"->"<<s.EndPlane<<endl
    <<"pl5.X,Y,Z="<<cp.X<<" , "<<cp.Y<<" , "<<cp.Z<<endl
    <<"d.X,Y,Z="<<d.X()<<" , "<<d.Y()<<" , "<<d.Z()<<endl
    <<"e.X,Y,Z="<<e.X()<<" , "<<e.Y()<<" , "<<e.Z()
    <<endl;
  
  for (UInt_t i=0;i<va.size();i++){
    //get references to the TVector3s
    TVector3& a=va[i];
    TVector3& b=vb[i];
    TVector3& c=vc[i];
    TVector3& intersect=vintersect[i];

    //work out d-a
    TVector3 d_a=d-a;
    
    Float_t denominator=b.Cross(c).Dot(-e);
    MAXMSG("MeuReco",logLevel,500)
      <<"Denominator="<<denominator<<endl;

    Float_t distance=999999;
    if (denominator!=0){
      Float_t numerator=b.Cross(c).Dot(d_a);
      Float_t lambda=numerator/denominator;

      intersect=d+(lambda*e);
      Float_t dZ=fabs(intersect.Z()-d.Z());

      distance=lambda*e.Mag();

      MAXMSG("MeuReco",logLevel,500)
        <<"i="<<i<<", numerator="<<numerator<<", Mag="<<e.Mag()
        <<", closest distance="<<distance<<endl
        <<"     in.X,Y,Z="<<intersect.X()<<" , "<<intersect.Y()
        <<" , "<<intersect.Z()
        <<", R="<<sqrt(pow(intersect.X(),2)+pow(intersect.Y(),2))
        <<", dZ="<<dZ
        <<" (pl="<<dZ/(5.94*Munits::cm)<<")"
        <<endl;

      //store the distance
      dist[i]=distance;

      //get closest distance, has to be positive
      if (distance<closestDist && distance>=0) {
        closesti=i;
        closestdZ=dZ;
        closestDist=distance;
      }
    }
    else{
      MAXMSG("MeuReco",logLevel,200)
        <<endl<<endl<<"Bad denominator, pl="<<pl<<", strip="<<cp.Strip
        <<", trk: "<<s.VtxPlane<<"->"<<s.EndPlane<<endl
        <<"pl5.X,Y,Z="<<cp.X<<" , "<<cp.Y<<" , "<<cp.Z<<endl
        <<"d.X,Y,Z="<<d.X()<<" , "<<d.Y()<<" , "<<d.Z()<<endl
        <<"e.X,Y,Z="<<e.X()<<" , "<<e.Y()<<" , "<<e.Z()
        <<endl;
      cout<<"CalcTrace: denominator=0!!!"<<endl<<endl;
    }
  }

  MAXMSG("MeuReco",logLevel,100)
    <<"ClosestDist="<<closestDist<<endl;

  if (closestDist<0.5 || (closestdZ<0.5 && closesti<8)){
    MAXMSG("MeuReco",logLevel,100)
      <<"ClosestDist="<<closestDist<<", dZ="<<closestdZ
      <<", i="<<closesti<<endl;
  }

  if (trace!=0 && traceZ!=0) {
    *trace=closestDist;
    if (closesti<8){
      *traceZ=closestdZ;
    }
    else *traceZ=999;//if it exits the SM end then don't calc traceZ
    //already have a 5 plane cut, which is good enough
  }

  return true;//return true if all ok
}

//......................................................................

---