ReMIdDedxFD_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////           
// \file    ReMIdDedxFD_module.cc
// \brief   A module for creating trees for FD cosmics. Based on the code for ND ReMIdDedx_module.cc
//          ReMId info is not available for tracks outside the spill window. The pieces of code
//          to extract the dE and dx info  are adapted from RecoMuon_module.cc
// \version $Id: ReMIdDedxFD_module.cc,v 2017/02/07  sgermani Exp $  
// \author  Stefano Germani - s.germani@ucl.ac.uk                                                   
//////////////////////////////////////////////////////////////////// 
#include <string>

// ROOT includes
#include "TVector3.h"
#include "TTree.h"
#include "TH1D.h"
#include "TH2D.h"
#include "TMath.h"

// Framework includes
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "Utilities/AssociationUtil.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Principal/SubRun.h"

// NOvA includes
//....................................
#include "Calibrator/Calibrator.h"
#include "PhotonTransport/FiberBrightness.h"
//....................................
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "RunHistory/service/RunHistoryService.h"
#include "ReMId/ReMId.h"
#include "NumuEnergy/NumuE.h"
#include "QEEventFinder/QePId.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "GeometryObjects/CellGeo.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/Track.h"
//................................
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"

#include "MCCheater/SimHit.h" 

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

#include "Simulation/Simulation.h"
#include "SummaryData/EventQuality.h"
#include "SummaryData/POTSum.h"
#include "SummaryData/SpillData.h"
#include "MCCheater/BackTracker.h"



namespace remid {

  // A module to analyze ReMId objects for pid training
  class ReMIdDedxFD : public art::EDAnalyzer {
  public:
    explicit ReMIdDedxFD(fhicl::ParameterSet const& pset);
    ~ReMIdDedxFD();

    void analyze(art::Event const& evt);

    void beginJob();
    void beginRun(const art::Run& run);
    void endJob();
    //virtual void beginSubRun(art::SubRun& sr);

  private:
    std::string fHitModuleLabel;
    std::string fSliceLabel;         ///<Where to find reconstructed slices
    std::string fTrackLabel;         ///<Where to find recondtructed tracks
    std::string fPidLabel;           ///<Where to find the pids
    std::string fGeneratorLabel;
    std::string fNumiBeamLabel;
    std::string fSpillQualLabel;
    std::string fNumuEnergyLabel;
    std::string fQepidLabel;

    TH1D*  fPOT;   ///<Histogram of POT passing spill cuts
    TH1D*  fSlicesPassingSelectionCuts; 
    TH1D*  fTrueSlicesPassingSelectionCuts; 


    double fTotalPOT;
    bool                     fFD;                      ///< Is detector FD?

    // variables used in filling the kNN training trees
    float fTrackLength;
    //float fProtonTrackLength;
    //float fDedxSep;
    //float fScatSep;
    //float fnuE;
    //float fvisE;
    //float frecoE;
    //float fMeasFrac;

    TTree* DedxSample;  ///< Tree for holding dE/dx variables needed for making signal sample histograms 

    // variables used in filling the sample histogram trees
    float fDedx;
    int fnDedx;
    int fnDedxUsed;
    float fDistDedx;
    //float fCosScat;
    //float fDistLastScat;
    //float fScatVar;
    //int fnScat;
    //float fDistScat;
    //int fnPlane;
    int fpdg;
    //int fintType;
    bool fVert;
    bool fUsed;
    float fzStart;
    float fzEnd;
    bool fPassOffTrackHits;
    bool fPassOffTrackEnergy;
    bool fPassHadEOnMuonTrack;

    // variables added to trees for Mini Prod mismatch studies (S. Germani) ....
    int fnCellsPerPlane;
    float fmyDedx;
    float fDe;
    int   fView;
    int   fPlane;
    int   fCell;
    float fPE;
    float fADC;
    float fPECorr;
    float fMIP;
    float fDx;
    float fxDead;
    float fxStart;
    float fyStart;
    float fxEnd;
    float fyEnd;
    float fxHit;
    float fyHit;
    float fzHit;
    float ftHit;
    float fxDir;
    float fyDir;
    float fzDir;
    float fmcDe;
    float fmcDeBirks;
    float fmcDx;
   
    float fTrueE;
    int   fFiberBrightness;

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

    art::ServiceHandle<geo::Geometry> fgeom;
    art::ServiceHandle<nova::dbi::RunHistoryService> frh;
    art::ServiceHandle<geo::LiveGeometry>  flivegeom;
    art::ServiceHandle<cheat::BackTracker>           bt;

    art::ServiceHandle<photrans::FiberBrightness>    fb;

    //novadaq::cnv::DetId fDetID;  ///< Detector ID in nova daq convention namespace typedef


    //unsigned int fdibfirst = 0;
    //unsigned int fdiblast = 0;

    bool ContainedEvent(art::Ptr<rb::Cluster> slice, art::Ptr<rb::Track> trk);


  };

}


namespace remid{


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

  ReMIdDedxFD::ReMIdDedxFD(fhicl::ParameterSet const& pset):
    EDAnalyzer(pset)
  {
    fHitModuleLabel = (pset.get< std::string >("HitModuleLabel"));
    fSliceLabel = (pset.get< std::string >("SliceLabel"));
    fTrackLabel = (pset.get< std::string >("TrackLabel"));
    //    fPidLabel   = (pset.get< std::string >("PidLabel"  ));
    fGeneratorLabel = pset.get< std::string >("GeneratorInput");
    fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
    //    fQepidLabel = pset.get< std::string >("QepidInput");
  }

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

  ReMIdDedxFD::~ReMIdDedxFD()
  {
  }

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

  void ReMIdDedxFD::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    
    DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
    DedxSample->Branch("dedx", &fDedx, "dedx/F");
    DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
    DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
    DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
    DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
    DedxSample->Branch("vert", &fVert, "vert/O");
    DedxSample->Branch("used", &fUsed, "used/O");
    DedxSample->Branch("zStart",&fzStart,"zStart/F");
    DedxSample->Branch("zEnd",&fzEnd,"zEnd/F");
    DedxSample->Branch("passOffTrackHits",&fPassOffTrackHits,"passOffTrackHits/O");
    DedxSample->Branch("passOffTrackEnergy",&fPassOffTrackEnergy,"passOffTrackEnergy/O");
    DedxSample->Branch("passHadEOnMuonTrack",&fPassHadEOnMuonTrack,"passHadEOnMuonTrack/O");

    DedxSample->Branch("ncellsPerPlane",   &fnCellsPerPlane,    "ncellsPerPlane/I");    
    DedxSample->Branch("mydedx",   &fmyDedx,    "mydedx/F");    
    DedxSample->Branch("de",       &fDe,        "de/F");    
    DedxSample->Branch("view",     &fView,      "view/I");    
    DedxSample->Branch("cell",     &fCell,      "cell/I");
    DedxSample->Branch("plane",    &fPlane,     "plane/I");
    DedxSample->Branch("adc",      &fADC,       "adc/F");    
    DedxSample->Branch("pe",       &fPE,        "pe/F");    
    DedxSample->Branch("pecorr",   &fPECorr,    "pecorr/F");    
    DedxSample->Branch("mip",      &fMIP,        "mip/F");    
    DedxSample->Branch("dx",       &fDx,        "dx/F");    
    DedxSample->Branch("deadx",    &fxDead,     "deadx/F");    
    DedxSample->Branch("xStart",   &fxStart,    "xStart/F");
    DedxSample->Branch("yStart",   &fyStart,    "yStart/F");
    DedxSample->Branch("xEnd",     &fxEnd,      "xEnd/F");  
    DedxSample->Branch("yEnd",     &fyEnd,      "yEnd/F");
    DedxSample->Branch("xHit",     &fxHit,      "xHit/F");  
    DedxSample->Branch("yHit",     &fyHit,      "yHit/F");
    DedxSample->Branch("zHit",     &fzHit,      "zHit/F");
    DedxSample->Branch("tHit",     &ftHit,      "tHit/F");
    DedxSample->Branch("xDir",     &fxDir,      "xDir/F");
    DedxSample->Branch("yDir",     &fyDir,      "yDir/F");
    DedxSample->Branch("zDir",     &fzDir,      "zDir/F");
    DedxSample->Branch("mcde",     &fmcDe,      "mcde/F");    
    DedxSample->Branch("mcdeBirks",&fmcDeBirks, "mcdeBirks/F");    
    DedxSample->Branch("mcdx",     &fmcDx,      "mcdx/F");    

    DedxSample->Branch("pdg",     &fpdg,        "pdg/I");    
    DedxSample->Branch("trueE",   &fTrueE,      "trueE/F");    
    DedxSample->Branch("fiberBrightness", &fFiberBrightness, "fiberBrightness/I");

    fPOT = tfs->make<TH1D>("POT",        ";;Total POT",1, 0.0, 1.0);
    fSlicesPassingSelectionCuts = tfs->make<TH1D>("SlicesPassingSelectionCuts",        ";;Slices Passing Selection Cuts",1, 0.0, 1.0);
    fTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("TrueSlicesPassingSelectionCuts",        ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);

    fTotalPOT = 0.0;

  }

  //......................................................................
  
  void ReMIdDedxFD::beginRun(const art::Run &r)

  { 

    fFD = false;
    
    novadaq::cnv::DetId detID = fgeom->DetId();
    if (detID == novadaq::cnv::kFARDET) fFD = true;

    return;
  }
  

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

  void ReMIdDedxFD::analyze(art::Event const& evt)
  {
    //Only written for FD
    if (!fFD) return;
    
    // get a geometry handle   ************************************** 
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<calib::Calibrator> cal;
    //***************************************************************

    //Getting event hits
    art::Handle< std::vector<rb::CellHit> > hithdl;
    evt.getByLabel(fHitModuleLabel, hithdl);

    art::PtrVector<rb::CellHit> allhits;
    for(size_t h = 0; h < hithdl->size(); ++h){
      art::Ptr<rb::CellHit> hit(hithdl, h);
      allhits.push_back(hit);
    }

    // get the slices out of the event
    art::Handle<std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fSliceLabel,slicecol);

    art::PtrVector<rb::Cluster> slicelist;
    
    for(unsigned int i = 0; i<slicecol->size();++i){
      art::Ptr<rb::Cluster>slice(slicecol,i);
      slicelist.push_back(slice);
    }

    // get  assosciations between slices and tracks
    art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);

    // get  assosciations between slices and numu energy
    art::FindManyP<numue::NumuE> fndmnyNumuEnergy(slicecol, evt, fNumuEnergyLabel);

    // loop over the slices
    for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
    
      // If slice is noise, skip it
      if(slicelist[iSlice]->IsNoise()){ continue; }

      // get all of the tracks associated to this slice
      std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
      // if no tracks move on
      if(tracks.size() == 0){ continue; }

      //Require precisely one track
      if(tracks.size() != 1){ continue; }
      

      //Require Best Track is 3D
      bool is3D = true;

      int bestidx    = 0; 
      if(!tracks[bestidx]->Is3D())    is3D = false;
      if (!is3D){ continue; }

      //Default all pass cuts to false
      fPassOffTrackHits = false;
      fPassOffTrackEnergy = false;
      fPassHadEOnMuonTrack = false;

      //Determine if tracks come from outside and stop
      art::ServiceHandle<geo::LiveGeometry> livegeo;

      /*
      double distStartToFront = livegeo->DistToFront( tracks[bestidx]->Start() );
      double distStartToTop   = livegeo->DistToTop(   tracks[bestidx]->Start() );
      double distStartToSide  = livegeo->DistToEdgeX( tracks[bestidx]->Start() );


      //physical starts
      if (distStartToFront < 0) continue;
      if (distStartToTop   < 0) continue;
      if (distStartToSide  < 0) continue;
      */

      // endpoint isn't contained
      double distEndToEdge = std::min(livegeo->DistToEdgeXY(tracks[bestidx]->Stop()), livegeo->DistToEdgeZ(tracks[bestidx]->Stop()));
      bool isStopper = true;
      if (distEndToEdge < 50) isStopper = false;
      
      // Process only contained tracks
      if (!isStopper) continue;
      
     //Filling muon track length properties
     fzStart = tracks[bestidx]->Start().Z();
     fzEnd = tracks[bestidx]->Stop().Z();
     fTrackLength = tracks[bestidx]->TotalLength();
     
     fxStart = tracks[bestidx]->Start().X();
     fyStart = tracks[bestidx]->Start().Y();
     fxEnd = tracks[bestidx]->Stop().X();
     fyEnd = tracks[bestidx]->Stop().Y();

     fxDir = tracks[bestidx]->Dir().X();
     fyDir = tracks[bestidx]->Dir().Y();
     fzDir = tracks[bestidx]->Dir().Z();
    
     double trklen = tracks[bestidx]->TotalLength();

     unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();

     // vector holds effective active track length from the end of the track
     std::vector<double> lengths(ntraj); 

     // loop over trajectory points starting from the end
     for(int itrj = ntraj-1; itrj >= 0; --itrj){
       TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
       if(itrj == (int)ntraj-1)  lengths[itrj] = 0.0;
       else {
         TVector3 lasttrajpoint = tracks[bestidx]->TrajectoryPoint(itrj+1);
         lengths[itrj] = (trajpoint-lasttrajpoint).Mag() + lengths[itrj+1];
       } // else
     } // for itrj
     
     unsigned int minPlane = tracks[bestidx]->MinPlane();
     unsigned int maxPlane = tracks[bestidx]->MaxPlane();
     int maxXPlane = -1;
     int maxYPlane = -1;
     double hitCosx[tracks[bestidx]->NXCell()];
     double hitCosy[tracks[bestidx]->NYCell()];

     for(unsigned int ihit = 0; ihit<tracks[bestidx]->NXCell(); ++ihit){
       if(tracks[bestidx]->XCell(ihit)->Plane() > maxXPlane){ maxXPlane = tracks[bestidx]->XCell(ihit)->Plane(); }
       // get direction with this x hit
       double xyzplane[3];
       geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
       double dz = geom->Plane(tracks[bestidx]->XCell(ihit)->Plane())->Cell(0)->HalfD();

       double minz = xyzplane[2]-dz;
       double maxz = xyzplane[2]+dz;
       unsigned int currtrj = 0;
       // find the matching trajectory point for this plane;
       for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
         if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
       } // for itaraj
       TVector3 direction = tracks[bestidx]->Dir();
       if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
         direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
       }
       direction = direction.Unit();
       hitCosx[ihit] = sqrt(1.0/(direction.X()*direction.X()+direction.Z()*direction.Z()))*direction.Z();
     } // for ihit XCell

     for(unsigned int ihit = 0; ihit<tracks[bestidx]->NYCell(); ++ihit){
       if(tracks[bestidx]->YCell(ihit)->Plane() > maxYPlane){ maxYPlane = tracks[bestidx]->YCell(ihit)->Plane(); }
       // get direction with this y hit
       double xyzplane[3];
       geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->GetCenter(xyzplane);
       double dz = geom->Plane(tracks[bestidx]->YCell(ihit)->Plane())->Cell(0)->HalfD();
       //       double dz = HalfD;
       double minz = xyzplane[2]-dz;
       double maxz = xyzplane[2]+dz;
       unsigned int currtrj = 0;
       // find the matching trajectory point for this plane;
       for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
         if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){currtrj = itraj;}
       }
       TVector3 direction = tracks[bestidx]->Dir();
       if(currtrj>0 && currtrj<tracks[bestidx]->NTrajectoryPoints()-2){
         direction = tracks[bestidx]->TrajectoryPoint(currtrj) - tracks[bestidx]->TrajectoryPoint(currtrj-1);
       }
       direction = direction.Unit();
       hitCosy[ihit] = sqrt(1.0/(direction.Y()*direction.Y()+direction.Z()*direction.Z()))*direction.Z();
     } // for ihit YCell

     ///

     // calculate dedx info for this track
     std::vector<double> dedxs;
     std::vector<double> des;
     std::vector<double> pecorrs;
     std::vector<double> mips;
     std::vector<double> dxs;
     std::vector<double> deads;
     std::vector<double> xdedxs;
     std::vector<unsigned int> dedxpls;

     std::vector<double> pes;
     std::vector<double> adcs;
     std::vector<double> xhits;
     std::vector<double> yhits;
     std::vector<double> zhits;
     std::vector<double> thits;
     std::vector<int> views;
     std::vector<int> cells;
     std::vector<int> planes;
     
     std::vector<double> mcdes;
     std::vector<double> mcdebs;
     std::vector<double> mcdxs;
     
     for(unsigned int iPlane = minPlane;iPlane<=maxPlane;++iPlane){

         double planeGeV    = 0;
         double planePECorr = 0;
         double planeMIP    = 0;

         int    planeView   = -1;
         int    plane       = -1;
         int    cell        = -1;
         double planeX      = 0;
         double planeY      = 0;
         double planeZ      = 0;
         double planeT      = 0;
         double planeADC    = 0;
         double planePE     = 0;
         
         double mcPlaneE    = 0;
         double mcPlaneEBirks = 0;
         double mcPathLength  = 0;
         double totcosinview = 0;
         double avgCosOtherView = 1.0;
         bool useOtherView = true;
         double nhits = 0;
         double xyzplane[3];
         geom->Plane(iPlane)->Cell(0)->GetCenter(xyzplane);
         double dz = geom->Plane(iPlane)->Cell(0)->HalfD();
         //      double dz = HalfD;
         double minz = xyzplane[2]-dz;
         double maxz = xyzplane[2]+dz;
         int matchpt = 0;
         int minCell = 10000;
         int maxCell = 0;
         double nhitsother = 0;
         int view = 0;
         int fbIndex = -1;
         fFiberBrightness = -1;
         // find the matching trajectory point for this plane;
         for(unsigned int itraj = 0;itraj<tracks[bestidx]->NTrajectoryPoints();++itraj){
             if(tracks[bestidx]->TrajectoryPoint(itraj).Z()>minz && tracks[bestidx]->TrajectoryPoint(itraj).Z()<maxz){matchpt = itraj;}
         }

         trklen = lengths[0];
         
         // find the track length up to this trajectory point
         double xlength = trklen - lengths[matchpt];
         

         if(minPlane == maxPlane){xlength = trklen/2;}
         if(tracks[bestidx]->Stop().Y()>tracks[bestidx]->Start().Y()){xlength = (trklen-xlength);}

         if(geom->Plane(iPlane)->View() == geo::kX){
           planeView = 0;
           for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
             if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane && tracks[bestidx]->XCell(iHit)->View() != geo::kY){
               const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
               planeADC+=chit->ADC();
               planePE +=chit->PE();
               cell  = chit->Cell();
               plane = chit->Plane();
               rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->XCell(iHit));
               if(minPlane == maxPlane){
                 rhit = cal->MakeRecoHit(*(tracks[bestidx]->XCell(iHit)), 0.25);
               }
               if(rhit.IsCalibrated()){
                 planeGeV+=rhit.GeV();
                 planePECorr+=rhit.PECorr();
                 planeMIP+=rhit.MIP();
                 planeX += rhit.X(); 
                 planeY += rhit.Y(); 
                 planeZ += rhit.Z(); 
                 planeT += rhit.T(); 
                 totcosinview+=hitCosx[iHit];
                 nhits+=1.0;
                 if(tracks[bestidx]->XCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->XCell(iHit)->Cell();}
                 if(tracks[bestidx]->XCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->XCell(iHit)->Cell();}

                 fbIndex = fb->getBrightnessIndex(iPlane,  tracks[bestidx]->XCell(iHit)->Cell() );

                 if(bt->HaveTruthInfo()){
                   const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->XCell(iHit);
                   const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
                   for(int j = 0; j < (int)flsHits.size(); j++){
                     mcPlaneE      += flsHits[j].GetEdep();
                     mcPlaneEBirks += flsHits[j].GetEdepBirks();
                     mcPathLength  += flsHits[j].GetTotalPathLength();

                   }//End loop over flshits in hit
                 } // if BackTracker

               }
             }
           } // for iHit
         
           // try to get the avgcos from the other view
           double nhitsbefore = 0;
           double nhitsafter = 0;
           double totCosBefore = 0;
           double totCosAfter = 0;
           for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
             if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane-1){
               ++nhitsbefore;
               totCosBefore+=hitCosy[iHit];
             }
             else if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane+1){
               ++nhitsafter;
               totCosAfter+=hitCosy[iHit];
             }
           }
           if(nhitsbefore > 0 && nhitsafter > 0){
             avgCosOtherView = ((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter));
                 nhitsother = (nhitsbefore+nhitsafter)/2;
           }
           else if(nhitsbefore > 0){ 
             avgCosOtherView = totCosBefore/nhitsbefore;
             nhitsother = nhitsbefore;
           }
           else if(nhitsafter > 0){ 
             avgCosOtherView = totCosAfter/nhitsafter;
             nhitsother = nhitsafter;
           }
           else{ useOtherView = false; }
           
         } else if(geom->Plane(iPlane)->View() == geo::kY){
           planeView = 1;
           view = 1;
           for(unsigned int iHit = 0; iHit<tracks[bestidx]->NYCell();++iHit){
             if(tracks[bestidx]->YCell(iHit)->Plane() == iPlane && tracks[bestidx]->YCell(iHit)->View() != geo::kX 
                && tracks[bestidx]->YCell(iHit)->View() != geo::kXorY){
               const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
               planeADC+=chit->ADC();
               planePE +=chit->PE();
               cell  = chit->Cell();
               plane = chit->Plane();

               rb::RecoHit rhit = tracks[bestidx]->RecoHit(tracks[bestidx]->YCell(iHit));
               if(minPlane == maxPlane){
                 rhit = cal->MakeRecoHit(*(tracks[bestidx]->YCell(iHit)), 0.25);
               }
               if(rhit.IsCalibrated()){
                 planeGeV+=rhit.GeV();
                 planePECorr+=rhit.PECorr();
                 planeMIP+=rhit.MIP();
                 planeX += rhit.X(); 
                 planeY += rhit.Y(); 
                 planeZ += rhit.Z(); 
                 planeT += rhit.T(); 
                 totcosinview+=hitCosy[iHit];
                 nhits+=1.0;
                 if(tracks[bestidx]->YCell(iHit)->Cell()<minCell){minCell = tracks[bestidx]->YCell(iHit)->Cell();}
                 if(tracks[bestidx]->YCell(iHit)->Cell()>maxCell){maxCell = tracks[bestidx]->YCell(iHit)->Cell();}

                 fbIndex = fb->getBrightnessIndex(iPlane,  tracks[bestidx]->YCell(iHit)->Cell() );

                 if(bt->HaveTruthInfo()){
                   const art::Ptr<rb::CellHit>& chit = tracks[bestidx]->YCell(iHit);
                   const std::vector< sim::FLSHit > flsHits = bt->HitToFLSHit(chit);
                   for(int j = 0; j < (int)flsHits.size(); j++){
                     mcPlaneE      += flsHits[j].GetEdep();
                     mcPlaneEBirks += flsHits[j].GetEdepBirks();
                     mcPathLength  += flsHits[j].GetTotalPathLength();

                   }//End loop over flshits in hit
                 } // if BackTracker

               } // if isCalibrated
             } //            
           } // for iHit           
     
         // try to get the avgcos from the other view
         double nhitsbefore = 0;
         double nhitsafter = 0;
         double totCosBefore = 0;
         double totCosAfter = 0;
         for(unsigned int iHit = 0; iHit<tracks[bestidx]->NXCell();++iHit){
           if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane-1){
             ++nhitsbefore;
             totCosBefore+=hitCosx[iHit];
           }
           else if(tracks[bestidx]->XCell(iHit)->Plane() == iPlane+1){
             ++nhitsafter;
             totCosAfter+=hitCosx[iHit];
           }
         }
         if(nhitsbefore > 0 && nhitsafter > 0){
           avgCosOtherView = TMath::Abs(((totCosAfter+totCosBefore)/(nhitsbefore+nhitsafter)));
           nhitsother = (nhitsbefore+nhitsafter)/2;
         }
         else if(nhitsbefore > 0){
           avgCosOtherView = TMath::Abs(totCosBefore/nhitsbefore);
           nhitsother = nhitsbefore;
         }
         else if(nhitsafter > 0){
           avgCosOtherView = TMath::Abs(totCosAfter/nhitsafter);
           nhitsother = nhitsbefore;
         }
         else{ useOtherView = false; }

         } // if kY 
         fnCellsPerPlane = nhits;
         if(nhits >0){

           if(nhits == 1) {
             fFiberBrightness = fbIndex;
           }

           double avgCos = TMath::Abs(totcosinview/nhits);
           double xyzp[3];
           geom->Plane(iPlane)->Cell(minCell)->GetCenter(xyzp);
           double xyzp1[3];
           geom->Plane(iPlane)->Cell(maxCell)->GetCenter(xyzp1);
           double planewidth = 2*geom->Plane(iPlane)->Cell(1)->HalfD();
           //     double planewidth = 2*HalfD;
           // calculate the amount of dead material to subtract off
           // first add up live material 
           double liveLength = 0;
           for(int i = minCell; i<maxCell+1;++i){
             // get the cell half width
             const geo::PlaneGeo* p = geom->Plane(iPlane);
             const geo::CellGeo* c = p->Cell(i);
             if(i == minCell || i == maxCell){ liveLength+=c->HalfW(); }
             else{ liveLength+=(2.0*c->HalfW() ); }
           }
           // dead material is total length - live length only if min and max cell are not the same
           double deadLength = 0;
           if(minCell != maxCell){ deadLength = (TMath::Abs(xyzp1[view]-xyzp[view])-liveLength); }
           double planeheight = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
           double planePathLength = TMath::Sqrt(planewidth*planewidth + planeheight*planeheight); 
           //      std::cout << "  ++++ " << iPlane << " " <<  planePathLength << " " << planewidth << " " << planeheight << std::endl; 
           // in the case of completely vertical track  be a bit careful
           if(avgCos == 0 || avgCos!=avgCos){ planePathLength = liveLength; }      
           if(useOtherView){
               double deltaz = planewidth;
               double deltav = planewidth*TMath::Tan(TMath::ACos(avgCos))-deadLength;
               if(avgCos == 0 || avgCos!=avgCos){
                   deltav = liveLength;
                   // if the track looks completely vertical in both views,
                   // than no part of the track length should be in z direction
                   if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){ deltaz = 0; }
               }

               // in opposite view dead length 
               double deltaov = (planewidth-deadLength)*TMath::Tan(TMath::ACos(avgCosOtherView));
               // in the case that the other view looks completely vertical scale up
               // the length between cells in this view by the number of hits in the other view
               // this case should really never happen
               if(avgCosOtherView == 0 || avgCosOtherView != avgCosOtherView){
                   deltaov = nhitsother*(TMath::Abs(xyzp1[view]-xyzp[view]))/nhits;
                   }
               planePathLength = TMath::Sqrt(deltaz*deltaz+deltav*deltav+deltaov*deltaov);
           }
           // also the case of completely vertical track
           if(minPlane == maxPlane){ planePathLength = liveLength; }
          
           xdedxs.push_back(xlength);

           assert(planeGeV/(planePathLength)>=0);
           dedxs.push_back(planeGeV/(planePathLength));
           des.push_back(planeGeV);
           pecorrs.push_back(planePECorr);
           mips.push_back(planeMIP);
           dxs.push_back(planePathLength);
           deads.push_back(deadLength);
           xhits.push_back(planeX/nhits);
           yhits.push_back(planeY/nhits);
           zhits.push_back(planeZ/nhits);
           thits.push_back(planeT/nhits);
           views.push_back(planeView);
           planes.push_back(plane);
           cells.push_back(cell);
           pes.push_back(planePE);
           adcs.push_back(planeADC);

           mcdes.push_back(mcPlaneE);
           mcdebs.push_back(mcPlaneEBirks);
           mcdxs.push_back(mcPathLength);
           
           //      std::cout << Form(" Dedx  %d    %2.6f  %2.6f  %2.6f", iPlane ,  planeGeV/(planePathLength), planeGeV, planePathLength ) <<  std::endl;
           } //nhits > 0         
     } // iPlane
     //.....
       
     //Looking at numu energy things
     numue::NumuE numuE;
     if (!fndmnyNumuEnergy.isValid()){ continue; }
     std::vector<art::Ptr<numue::NumuE> > energyVec = fndmnyNumuEnergy.at(iSlice);
     if (energyVec.empty()){ continue; }
     numuE = *energyVec[0];

     //Require decent energy object; not too much hadronic energy in muon catcher
     if (numuE.E() == -5){ continue; }
     //Not much hadronic energy contamination on muon track
     if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }

     //Slice energy mostly on one tracks
     double calEDiff = (slicelist[iSlice]->CalorimetricEnergy())-(tracks[bestidx]->CalorimetricEnergy());
     if (calEDiff < 0.1){ fPassOffTrackEnergy = true; }

     //Slice hits mostly on one track
     int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
     if (hitDiff < 4){ fPassOffTrackHits = true; }
 
     fpdg = 0;
    
     if (bt->HaveTruthInfo()){
         const std::vector<const sim::Particle*> test_mu = bt->HitsToParticle(tracks[bestidx]->AllCells());
         if (!test_mu.empty()) {
           fpdg   = test_mu[0]->PdgCode();
           fTrueE = test_mu[0]->E();
         }
     } // if HaveTruthInfo   

     if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){ 
       // get all the dedx variables filled
       fnDedx = des.size() -1;      // pid->NDedx()-1;
       fnDedxUsed = des.size(); // pid->NMeasurementPlanes();
       // loop over all of the dE/dx measurements
       //for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
       for( int idedx = 0; idedx<fnDedx; ++idedx){

         fDe     = des[idedx];
         fPECorr = pecorrs[idedx];
         fMIP    = mips[idedx];
         fDx     = dxs[idedx];
         fxDead  = deads[idedx];
         fmyDedx = fDe/fDx;

         fxHit   = xhits[idedx];
         fyHit   = yhits[idedx];
         fzHit   = zhits[idedx];
         ftHit   = thits[idedx];
         fView   = views[idedx];
         fPlane  = planes[idedx];
         fCell   = cells[idedx];
         fADC    = adcs[idedx];
         fPE     = pes[idedx];

         fmcDe      = mcdes[idedx];
         fmcDeBirks = mcdebs[idedx];
         fmcDx      = mcdxs[idedx];

         // ............
         // get the distance from the end of the track
         // fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[bestidx]->TotalLength();
         fDistDedx = (1.0 -xdedxs[idedx]/trklen)*tracks[bestidx]->TotalLength();
         // fill the dedx sample tree
         DedxSample->Fill();
       }//Loop over DedxSample

     }//If passes contaiment  (now  temporarely always true)
     
    } // end loop over slices
        
  }


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

  void ReMIdDedxFD::endJob()
  {
    fPOT->Fill(.5, fTotalPOT);
  }


  //.......................................................................
  bool ReMIdDedxFD::ContainedEvent(art::Ptr<rb::Cluster> slice, art::Ptr<rb::Track> trk)
  {
    
    bool contained = false;
    
    // check that this is a contained event
    
    // 2) how close is the track to the edge of the detector
    // in planes
    //SL - NOPE, need to ask slice, not track
    //    unsigned int firstpl = slice->MinPlane();
    //    unsigned int lastpl = slice->MaxPlane();
    // in projected cells, just define directions here actually get projected cell below
    TVector3 dirbwk = -trk->Dir();
    TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();

    // 3) hits in slice
    int sliceNHit = slice->NCell();

    // 4) slice mean time
    //    double slcTime = slice->MeanTNS();

    // 5) hits on track
    int trackNHit = trk->NCell();

    //std::cout<<"In contained event loop!"<<std::endl;
    // cut the slice number of cells from edge
    if((sliceNHit > 20)&&(trackNHit > 5)){
      int fwdcell = flivegeom->ProjectedLiveCellsToEdge(trk->Stop(), dirfwd);
      //      int bwkcell = flivegeom->ProjectedLiveCellsToEdge(trk->Start(),dirbwk);

      // also make sure the track doesn't start too close to the end of the active region
      if(fwdcell > 4)   contained = true; 
    } // end if track ends projections contained
        
    return contained;
    
  }//  ContainedEvent
 

}


namespace remid{

  DEFINE_ART_MODULE(ReMIdDedxFD)

}