DiFShowerFinder_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       DiFShowerFinder
// Module Type: producer
// File:        DiFShowerFinder_module.cc
//
// Generated at Thu Aug 13 13:48:09 2015 by Nitin Yadav using artmod
// from cetpkgsupport v1_08_06.
// 
//
// Updated to include use at ND 10/30/2019
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "canvas/Utilities/InputTag.h"

#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/CellGeo.h"

#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"

#include "Utilities/func/MathUtil.h"
#include "Utilities/AssociationUtil.h"

#include <vector>
#include <algorithm> 

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

namespace dif
{

class DiFShowerFinder : public art::EDProducer
{
  public:
   explicit DiFShowerFinder(fhicl::ParameterSet const & pset);
   virtual ~DiFShowerFinder();

  //void reconfigure(const fhicl::ParameterSet& pset); //redid hwo parameters are read in
    void beginJob() override;
    void beginRun(art::Run&);
    virtual void produce(art::Event & evt) override;
  
    void findShowerByReco(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane);
    bool inFiducial(double x, double y, double z);
    bool trackStartContained(art::Ptr<rb::Track> &trk);
    bool trackEndContained(art::Ptr<rb::Track> &trk);
    int adjustPlane(art::Ptr<rb::Track> & trk, int p);
    int adjustPlane_end(art::Ptr<rb::Track> & trk, int p);
    bool eparm(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> & trk, int startplane, int endPlane);
    double distancefromtrack( art::Ptr<rb::CellHit> cHit, art::Ptr<rb::Track> &trk);
    void muonstub(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> & trk, int & startplane);

    void endJob() override;


  private:
    std::string fClusterLabel;
    std::string fTrackLabel;

    bool fWriteDiFCluster;
    bool fWriteParentTrackCluster;
    bool fWriteParentSliceCluster;

    bool fShowerByReco;
    bool fTrackForwardGoing; //Assume the track end is the farthest in z instead of lowest in y?
    float fAngleCut;
    float fEnergyCut;
    float fFiducialX;
    float fFiducialY;
    float fFiducialZ;
    float fMinDistToEdge;
    int fPlaneMin;
    int fPlaneMax;
    int fMinTrackLength;//"Optimized" cut on track quality
    int fMaxTrackTail; //"Optimized" cut on Brem background

    fhicl::ParameterSet fPSetND, fPSetFD, fPSetTB;

    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<geo::LiveGeometry> livegeom;
};

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

  DiFShowerFinder::DiFShowerFinder(fhicl::ParameterSet const & pset) :

    
    fClusterLabel      ( pset.get< std::string >("ClusterLabel")),
    fTrackLabel        ( pset.get< std::string >("TrackLabel")),
    fWriteDiFCluster   ( pset.get< bool        >("WriteDiFCluster")),
    fWriteParentTrackCluster ( pset.get< bool  >("WriteParentTrackCluster")),
    fWriteParentSliceCluster ( pset.get< bool  >("WriteParentSliceCluster")),
    fShowerByReco      ( pset.get< bool        >("ShowerByReco")),
    fTrackForwardGoing (false),
    fAngleCut          ( pset.get< float       >("AngleCut")),
    fEnergyCut         ( pset.get< float       >("EnergyCut")),
    fFiducialX         (0),
    fFiducialY         (0),
    fFiducialZ         (0),
    fMinDistToEdge     (0),
    fPlaneMin          (0),
    fPlaneMax          (0),
    fMinTrackLength    (0),
    fMaxTrackTail      (0),
    fPSetND            (pset.get<fhicl::ParameterSet>("nd")),
    fPSetFD            (pset.get<fhicl::ParameterSet>("fd")),
    fPSetTB            (pset.get<fhicl::ParameterSet>("tb"))
  {
    produces<std::vector<rb::Cluster> >();
    if(fWriteDiFCluster) produces< art::Assns<rb::Cluster, rb::Track> >();
  }
//-----------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////

  DiFShowerFinder::~DiFShowerFinder() { }

  void DiFShowerFinder::beginJob()
  {
  }

  void DiFShowerFinder::beginRun(art::Run&)
  {
    fhicl::ParameterSet pset;
    switch(geom->DetId()){
    case novadaq::cnv::kNEARDET:
      pset = fPSetND;
      break;
    case novadaq::cnv::kFARDET:
      pset = fPSetFD;
      break;
    case novadaq::cnv::kTESTBEAM:
      pset = fPSetTB;
      break;
    default:
      assert(0 && "Unknown detector");
    }

  fTrackForwardGoing = pset.get< bool        >("TrackForwardGoing");
  fFiducialX         = pset.get< float       >("FiducialX");
  fFiducialY         = pset.get< float       >("FiducialY");
  fFiducialZ         = pset.get< float       >("FiducialZ");
  fMinDistToEdge     = pset.get< float       >("MinDistToEdge");
  fPlaneMin          = pset.get< int         >("PlaneMin");
  fPlaneMax          = pset.get< int         >("PlaneMax");
  fMinTrackLength    = pset.get< int         >("MinTrackLength");
  fMaxTrackTail      = pset.get< int         >("MaxTrackTail");

  if((fWriteDiFCluster || fWriteParentTrackCluster) && fWriteParentSliceCluster){
    mf::LogWarning("DiFShowerFinder") << "DifShowerFinder: WARNING! your configured options will include some hits in more than one output cluster. Recommend WriteParentSliceCluster is not used in combination with the other options" << std::endl
                                      << pset.to_string();
  }

  if(!fWriteDiFCluster && !fWriteParentTrackCluster && !fWriteParentSliceCluster){
    mf::LogWarning("DiFShowerFinder") << "DifShowerFinder: WARNING! you haven't configured any output. You probably want one or more of the Write*Cluster settings." << std::endl
                                      << pset.to_string();
  }

  }


void DiFShowerFinder::produce(art::Event & evt)
{ 
  art::Handle< std::vector<rb::Cluster> > slicecol; 
  evt.getByLabel(fClusterLabel, slicecol);

  art::FindManyP<rb::Track> fmt(slicecol, evt, fTrackLabel);

  std::unique_ptr< art::Assns<rb::Cluster, rb::Track>   > trackAssassin(new art::Assns<rb::Cluster, rb::Track>);
  std::unique_ptr< std::vector<rb::Cluster> > showercol(new std::vector<rb::Cluster>);
  // Downstream modules love their noise slices
  showercol->emplace_back();
  showercol->back().SetNoise(true);


  // Loop over slices 
  for( unsigned int iSlice =0; iSlice < slicecol->size(); ++iSlice){
    art::Ptr<rb::Cluster> slc(slicecol, iSlice);

    if(slc->NCell()<=0) continue;

    bool goodSlice = true;

    if (slc->MeanTNS () <25000 || slc->MeanTNS () > 475000) {goodSlice =false;}

    std::vector <art::Ptr<rb::Track> > tracks = fmt.at(iSlice);

    int startPlane = -1;
    int endPlane = -1;

    double disToBadX=-1;
    double disToBadY=-1;
    double disToBadZ=-1;
    bool isShower = false; 
    bool isCandTrack  = false; 

    if(!slc->IsNoise() && tracks.size()==1 && tracks[0]->Is3D() && goodSlice){

      bool startContained = false;
      bool endContained =false; 

      //Decide if this track is a candidat track?
      if(int(tracks[0]->MaxPlane() ) >fPlaneMin && int(tracks[0]->MinPlane() ) <fPlaneMax){ 

        startContained = trackStartContained(tracks[0]);
        endContained = trackEndContained(tracks[0]);

        //part of "optimized" cuts
        int nPlaneTrack = 0;

        if(int(tracks[0]->MaxPlane())>fPlaneMax) nPlaneTrack = fPlaneMax - int(tracks[0]->MinPlane());
        else if(int(tracks[0]->MinPlane())<fPlaneMin) nPlaneTrack = int(tracks[0]->MaxPlane()) - fPlaneMin;
        else nPlaneTrack = tracks[0]->MaxPlane() - tracks[0]->MinPlane();

        disToBadX=livegeom->DistToEdgeX(tracks[0]->Stop());
        disToBadY=livegeom->DistToEdgeY(tracks[0]->Stop());
        disToBadZ=livegeom->DistToEdgeZ(tracks[0]->Stop()); 


        if(startContained && !endContained) isCandTrack =false;
        else{
          
          if(tracks[0]->Dir().Z()<fAngleCut) isCandTrack = false;
        else {

          if(nPlaneTrack<fMinTrackLength) isCandTrack =false; //"optimized" cut
          else {

              if(abs(disToBadX)<fMinDistToEdge || abs(disToBadY)<fMinDistToEdge  || abs(disToBadZ)<fMinDistToEdge) isCandTrack = false;
               else {
                 isCandTrack = true;
               }
          }
        }
          
        
        } 
    } // candidate track ends   


      if(isCandTrack){ 

        if(fShowerByReco && tracks[0]->Dir().Z()>0){

          findShowerByReco(slc, tracks[0], startPlane, endPlane);
          if(startPlane>0 && endPlane>0) {

             int parameter = std::abs(int(tracks[0]->MaxPlane()) - endPlane);

            if(parameter<fMaxTrackTail){ //"optimized" cut
              if(!eparm(slc,tracks[0],startPlane, endPlane)){

                isShower=true;
                // endPlane = tracks[0]->MaxPlane();
                
              } //eparm
             } //param        
          } // startplane and endplane >0

        } // fShowerByReco

    } // isCandTrack
  } //good slice ends end if track size is 1

  if(isShower){

    std::cout<<"DiF Shower Found! Run: "<<evt.run()<<" Subrun: "<<evt.subRun()<<" Event: "<<evt.event()<<"\n";

    if(fWriteDiFCluster){
      // removing the muon stub
      muonstub(slc, tracks[0],startPlane );

      rb::Cluster shower;
   
      for(unsigned int iCell =0; iCell<slc->NCell(); ++iCell){
        art::Ptr<rb::CellHit> chit = slc->Cell(iCell);

        if(chit->Plane()<startPlane || chit->Plane()>endPlane) continue;

        if(distancefromtrack(chit, tracks[0])<=80){
          shower.Add(chit);
        }
      }//looping cells   
      
      //non 3D clusters aren't useful and break things down the line
      if(shower.Is3D()){
        showercol->push_back(shower);
        if(fWriteDiFCluster) util::CreateAssn( *this, evt, *showercol, tracks[0], *trackAssassin );
      }
    }

    if(fWriteParentTrackCluster){
      showercol->push_back(*tracks[0]);
    }

    if(fWriteParentSliceCluster){
      showercol->push_back(*slc);
    }
  } //isShower

} // looping over slice ends 

  evt.put(std::move(showercol));
  if(fWriteDiFCluster) evt.put(std::move(trackAssassin));
}//end produce

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

bool DiFShowerFinder::inFiducial(double x, double y, double z){
  bool isContained = false;

  double xyzMax[3], xyzMin[3];
  geom->Plane(fPlaneMax)->Cell(0)->GetCenter(xyzMax);
  geom->Plane(fPlaneMin)->Cell(0)->GetCenter(xyzMin);
  double zMax = xyzMax[2] + geom->Plane(fPlaneMax)->Cell(0)->HalfD();
  double zMin = xyzMin[2] - geom->Plane(fPlaneMin)->Cell(0)->HalfD();

  if(x < (geom->DetHalfWidth() - fFiducialX) && x > (-geom->DetHalfWidth() + fFiducialX) &&
      y < (geom->DetHalfHeight() - fFiducialY) && y > (-geom->DetHalfHeight() + fFiducialY) &&
      z < (zMax - fFiducialZ) && z > zMin + fFiducialZ
    )

    isContained = true;
  return isContained;

}
//-----------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////

bool DiFShowerFinder::trackStartContained(art::Ptr<rb::Track> &trk){

  return inFiducial(trk->Start().X(), trk->Start().Y(), trk->Start().Z());
}

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

bool DiFShowerFinder::trackEndContained(art::Ptr<rb::Track> &trk){

  double ymin = 10000;
  double zmax = 0;

  bool isContained = false;
  int nTrajPt = trk->NTrajectoryPoints();
 
  for(int iTrajPt=0; iTrajPt<nTrajPt; iTrajPt++){
    double xTrajpt =  trk->TrajectoryPoint(iTrajPt).X();
    double yTrajpt =  trk->TrajectoryPoint(iTrajPt).Y();
    double zTrajpt =  trk->TrajectoryPoint(iTrajPt).Z();

    //assume forward track if specified
    if(fTrackForwardGoing == true){
      if(zTrajpt>zmax){
        isContained = inFiducial(xTrajpt, yTrajpt, zTrajpt);
        zmax = zTrajpt;
      }
    }
    
    //assume downward track otherwise
    else{
      if(yTrajpt<ymin){
        isContained = inFiducial(xTrajpt, yTrajpt, zTrajpt);
        ymin = yTrajpt;
      }     
    }

  }

  return isContained;
}
//-----------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////

int DiFShowerFinder::adjustPlane(art::Ptr<rb::Track> & trk, int p){

 int pn=p;

  int minPlane = trk->MinPlane();
  int maxPlane = trk->MaxPlane();
  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
  std::map<unsigned int, double> ePlane;

  ePlane.clear();

  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
    const rb::CellHit* cHit = trk->Cell(iCell).get();
    const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
    if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
    if( rHit.IsCalibrated() ){
      ePlane[cHit->Plane()] += rHit.GeV();
    }
  }//end loop over cells

  for(int i=pn; i >=minPlane; i--){
    double mipDeDx  = 0.00157;
    double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
    pn=i;
    double avg = (ePlane[i-2] + ePlane[i] + ePlane[i-1])/3;//the average of 3 planes
    //double avg = (ePlane[i-3] + ePlane[i-2] + ePlane[i-1]+ ePlane[i])/4;//the average of 3 planes
    //if(ePlane[i]<0.5*mipPlane) break;
    //if(avg<=1.0*mipPlane)
    if(avg<=0.9*mipPlane)

      break;
  }

  return pn;
}
//-----------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////

int DiFShowerFinder::adjustPlane_end(art::Ptr<rb::Track> & trk, int p){
  int pln=p;

  int minPlane = trk->MinPlane();
  int maxPlane = trk->MaxPlane();
  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
  std::map<unsigned int, double> ePlane;

  ePlane.clear();

  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
    const rb::CellHit* cHit = trk->Cell(iCell).get();
    const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
    if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
    if( rHit.IsCalibrated() ){
      ePlane[cHit->Plane()] += rHit.GeV();
    }
  }//end loop over cells

  for(int i=pln; i <=maxPlane-2; i++){
    double mipDeDx  = 0.00157;
    double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
    pln=i;
    double avg = (ePlane[i+2] + ePlane[i] + ePlane[i+1])/3;//the average of 3 planes
    //double avg = (ePlane[i-3] + ePlane[i-2] + ePlane[i-1]+ ePlane[i])/4;//the average of 3 planes
    //if(ePlane[i]<0.5*mipPlane) break;
    if(avg<=0.9*mipPlane)

      break;
  }

  if(pln+3 <= int(trk->MaxPlane())) return pln+3;
  else return pln;
}

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

bool DiFShowerFinder::eparm(art::Ptr<rb::Cluster> & slc,art::Ptr<rb::Track> & trk, int startPlane , int endPlane){

  int minPlane = trk->MinPlane();
  int maxPlane = trk->MaxPlane();
  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;
  std::map<unsigned int, double> ePlane;
  // std::map<unsigned int, double> ePlanec;

  ePlane.clear();
  //ePlanec.clear();

  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
    const rb::CellHit* cHit = trk->Cell(iCell).get();
    const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));
    if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;
    if( rHit.IsCalibrated() ){
      ePlane[cHit->Plane()] += rHit.GeV();
    }
  }//end loop over cells

  bool eparm=false;
  //int eparmcount=0;
 
  for(int i=endPlane; i >endPlane-7; i--){
    double mipDeDx  = 0.00157;
    double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
    
    for(int k=0;k<=3;k++){
      if(k==0) eparm=(ePlane[i-k] >0.5*mipPlane && ePlane[i-k]<=1.5*mipPlane);
      else eparm*= (ePlane[i-k] >0.5*mipPlane && ePlane[i-k]<=1.5*mipPlane);
      if(!eparm) break;
    }
    if(eparm) break;
  }

  return eparm;
}
//-----------------------------------------------------------------------------
/////////////////////////////////////////////////////////////////////////////////

void DiFShowerFinder::findShowerByReco(art::Ptr<rb::Cluster> & slc, art::Ptr<rb::Track> &trk, int &startPlane, int &endPlane){

  int minPlane = trk->MinPlane();
  int maxPlane = trk->MaxPlane();
  if(minPlane<fPlaneMin) minPlane = fPlaneMin;
  if(maxPlane>fPlaneMax) maxPlane = fPlaneMax;

  std::map<unsigned int, double> ePlanec;
  std::map<unsigned int, double> ePlane;
  std::map<unsigned int, double> ePlaneAve;

  ePlane.clear();
  ePlanec.clear();

for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
    const rb::CellHit* cHit = slc->Cell(iCell).get();
    const rb::RecoHit rHit = slc->RecoHit(slc->Cell(iCell));

    if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;

    if( rHit.IsCalibrated() ){
      ePlanec[cHit->Plane()] += rHit.GeV();
      //  fHEdep->Fill(rHit.GeV());
      //  fHADC->Fill(cHit->ADC(0));
      //  fHPE->Fill(cHit->PE());
      //  fHPECorr->Fill(rHit.PECorr());
    }
  }//end loop over cells



  for(unsigned int iCell=0; iCell<trk->NCell(); ++iCell){
    const rb::CellHit* cHit = trk->Cell(iCell).get();
    const rb::RecoHit rHit = trk->RecoHit(trk->Cell(iCell));

    if(cHit->Plane()<fPlaneMin || cHit->Plane()>fPlaneMax) continue;

    if( rHit.IsCalibrated() ){  
      ePlane[cHit->Plane()] += rHit.GeV();
      //  fHEdep->Fill(rHit.GeV());
      //  fHADC->Fill(cHit->ADC(0));
      //  fHPE->Fill(cHit->PE());
      //  fHPECorr->Fill(rHit.PECorr());
    }
  }//end loop over cells

  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
        //opt stuff
     //if(iPlane == minPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane] + ePlane[iPlane+1]);
     //else if(iPlane == maxPlane) ePlaneAve[iPlane] = 0.5*(ePlane[iPlane-1] + ePlane[iPlane]);
    // else ePlaneAve[iPlane] = (ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1])/3;//the average of 3 planes

     //ePlaneAve[iPlane] = ePlane[iPlane];


  if(iPlane == minPlane) ePlaneAve[iPlane] = (ePlane[iPlane] + ePlane[iPlane+1]  + ePlane[iPlane+2])/3;
    else if(iPlane == minPlane+1) ePlaneAve[iPlane] = ( ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1]+ePlane[iPlane+2] )/4;
    else if(iPlane == maxPlane) ePlaneAve[iPlane] = (ePlane[iPlane-2] + ePlane[iPlane-1] + ePlane[iPlane])/3;
    else if(iPlane == maxPlane-1) ePlaneAve[iPlane] = (ePlane[iPlane-2]+ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1])/4;
    else ePlaneAve[iPlane] = (ePlane[iPlane-2] + ePlane[iPlane-1] + ePlane[iPlane] + ePlane[iPlane+1]  + ePlane[iPlane+2])/5;//the average of 5 planes

  }

  bool isShower = false;
  bool showerStart = false;
  bool showerEnd = false;

  const int nPlaneStart =5;// it is 2 
  const int nPlaneEnd = 2;//opt stuff change from 2 to 1

  double eshower = 0;
  for(int iPlane=minPlane; iPlane<=maxPlane; iPlane++){
    double mipDeDx  = 0.00157;
    double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
    if(!isShower){//if it's not already in a shower, try to find one   
      if(iPlane<=(maxPlane-nPlaneStart))
        for(int k=0; k<nPlaneStart; ++k){
          // if(k==0) showerStart= (ePlane[iPlane+k]>3.5*mipPlane);
          // else showerStart *= (ePlane[iPlane+k]>3.5*mipPlane);
          if(k==0) showerStart= (ePlane[iPlane+k]>2.5*mipPlane);//opt stuff 3.0 to 1.5
          else showerStart *= (ePlane[iPlane+k]>2.5*mipPlane);  //opt stuff
          if(!showerStart) break;
        }

      else showerStart = false;


      if(showerStart){
        isShower = true; 
        startPlane=adjustPlane(trk, iPlane ) ;
      //    startPlane = iPlane;
      }//end if shower start   
    }//end if not shower  

    else{
      for(int k=0; k<nPlaneEnd; k++)
        if((iPlane+k)<=maxPlane){

          if(k==0) showerEnd = (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.01*mipPlane); //opt stuff change from 1.5 to 2.5
          else showerEnd *= (ePlaneAve[iPlane+k]<1.5*mipPlane && ePlaneAve[iPlane+k]>0.01*mipPlane);

          if(!showerEnd) break;
        }
      if(showerEnd){
        isShower = false;
        endPlane= adjustPlane_end(trk,iPlane);
//        endPlane = iPlane;
      }

    }//end else

  }//end loop over track planes


  if(startPlane>0 && endPlane>0 && startPlane>fPlaneMin && endPlane<fPlaneMax){
    for(int iPlane = startPlane; iPlane <= endPlane; iPlane++) {
      //double mipDeDx  = 0.00157;
      //double mipPlane = 2*mipDeDx*geom->Plane(iPlane)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
      eshower += ePlanec[iPlane] ;//- mipPlane;
    }

    // fHEShower->Fill(eshower);
  }

  if(eshower < fEnergyCut) startPlane = -1;
  return;
}//end findShowerByReco

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

double DiFShowerFinder::distancefromtrack( art::Ptr<rb::CellHit> cHit, art::Ptr<rb::Track> &trk)
{
  TVector3 start  = trk->Start();
  TVector3 dir    = trk->Dir();

  TVector3 xyz;
  geom->Plane(cHit->Plane())->Cell(cHit->Cell())->GetCenter(xyz);
  const TVector3 dxyz = (cHit->View() == geo::kX) ? TVector3(0, 1, 0): TVector3(1, 0, 0);

  double di;
  return sqrt(geo::ClosestApproach(start, start+dir, xyz, xyz+dxyz, &di));

} // end of DistanceToCore 

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

void DiFShowerFinder::muonstub(art::Ptr<rb::Cluster> & slc, art::Ptr<rb::Track> & trk,int  & startPlane ){

  int sp=startPlane;
  //int minPlane = slc->MinPlane();
  int maxPlane = slc->MaxPlane();

  int spp=sp+7;
  if(spp>maxPlane) return ;
  std::map<unsigned int, double> ePlane;

  ePlane.clear();

  for(unsigned int iCell=0; iCell<slc->NCell(); ++iCell){
    const rb::CellHit* cHit = slc->Cell(iCell).get();
    const rb::RecoHit rHit = slc->RecoHit(slc->Cell(iCell));
    if( rHit.IsCalibrated() ){
      ePlane[cHit->Plane()] += rHit.GeV();
    }
  }//end loop over cells

  bool eparm=false;
  //int eparmcount=0;


  for(int i=sp; i <sp+7; i++){
    double mipDeDx  = 0.00157;
    double mipPlane = 2*mipDeDx*geom->Plane(i)->Cell(0)->HalfD()/fabs(trk->Dir().Z());
  
    for(int k=0;k<=4;k++){
      if(i+k>maxPlane) continue;
      if(k==0) eparm=(ePlane[i+k] >0.2*mipPlane && ePlane[i+k]<=2.5*mipPlane);
      else eparm*= (ePlane[i+k] >0.2*mipPlane && ePlane[i+k]<=2.5*mipPlane);
      if(!eparm) break;
    }

    if(eparm){ startPlane=i+4;}
  }
}

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

void DiFShowerFinder::endJob() { }

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

DEFINE_ART_MODULE(DiFShowerFinder)

} // end namespace