ReMIdDedxRock_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////           
// \file    ReMIdDedxRock_module.cc
// \brief   A module for creating trees which can be used to study dE/dx for Rock Muons
//          Adapted from ReMIdDedx_module.cc
// \version $Id: ReMIdDedxRock_module.cc,v 2017/06/20  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_root_io/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 ReMIdDedxRock : public art::EDAnalyzer {
  public:
    explicit ReMIdDedxRock(fhicl::ParameterSet const& pset);
    ~ReMIdDedxRock();

    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*  fEventsPassingSpillCuts; 
    TH1D*  fSlicesPassingSelectionCuts; 
    TH1D*  fTrueSlicesPassingSelectionCuts; 
    TH1D*  fSuperTrueSlicesPassingSelectionCuts; 
    TH2D*  fTrueProtonEFromTrackLength;
    TH2D*  fSuperTrueProtonEFromTrackLength;

    double fTotalPOT;
    bool                     fND;                      ///< Is detector ND?

    // 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 
    //TTree* ProtonDedxSample;  ///< 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 fProtonDedx;
    //int fProtonnDedx;
    //int fProtonnDedxUsed;
    //float fProtonDistDedx;   
    //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 fProtonVert;
    //bool fProtonUsed;
    //float fProtonzStart;
    //float fProtonzEnd;
    bool fPassOffTrackHits;
    bool fPassOffTrackEnergy;
    bool fPassHadEOnMuonTrack;
    bool fPassRockMuonCut;
    bool fPassVertX;
    bool fPassVertY;
    bool fPassVertZ;
    bool fPassProtonRemid;
    bool fPassDedxRatio;
    bool fPassEDiffZTest;
    bool fTrueNumuCCQE;
    bool fTrueNumuCCWithProton;
    bool fSuperTrueNumuCCWithProton;

    // variables added to trees for detailed dE/dx studies
    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 fxEnd;
    float fyStart;
    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;

    std::vector <int> fMultiCell;
    std::vector <float> fMultiPE;
    std::vector <float> fMultiADC;
    std::vector <float> fMultiPECorr;

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

    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{


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

  ReMIdDedxRock::ReMIdDedxRock(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");
    fNumiBeamLabel  = pset.get< std::string >("NuMIBeamInput");
    fSpillQualLabel = pset.get< std::string >("SpillQualInput");
    fNumuEnergyLabel = pset.get< std::string >("NumuEnergyInput");
    fQepidLabel = pset.get< std::string >("QepidInput");
  }

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

  ReMIdDedxRock::~ReMIdDedxRock()
  {

  }

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

  void ReMIdDedxRock::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("passRockMuonCut",&fPassRockMuonCut,"passRockMuonCut/O");
    DedxSample->Branch("passVertX",&fPassVertX,"passVertX/O");
    DedxSample->Branch("passVertY",&fPassVertY,"passVertY/O");
    DedxSample->Branch("passVertZ",&fPassVertZ,"passVertZ/O");
    DedxSample->Branch("passProtonRemid",&fPassProtonRemid,"passProtonRemid/O");
    DedxSample->Branch("passDedxRatio",&fPassDedxRatio,"passDedxRatio/O");
    DedxSample->Branch("passEDiffZTest",&fPassEDiffZTest,"passEDiffZTestX/O");
    DedxSample->Branch("trueNumuCCQE",&fTrueNumuCCQE,"trueNumuCCQE/O");

    // Branches added wrt standard ReMIdDedx module for detailed dE/dx studies
    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("xEnd",     &fxEnd,      "xEnd/F");  
    DedxSample->Branch("yStart",   &fyStart,    "yStart/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");

    DedxSample->Branch("multiCell",     &fMultiCell );
    DedxSample->Branch("multiADC",      &fMultiADC  );
    DedxSample->Branch("multiPE",       &fMultiPE   );
    DedxSample->Branch("multiPECorr",   &fMultiPECorr);


    fPOT = tfs->make<TH1D>("POT",        ";;Total POT",1, 0.0, 1.0);
    fEventsPassingSpillCuts = tfs->make<TH1D>("EventsPassingSpillCuts",        ";;Events Passing Spill Cuts",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);
    fSuperTrueSlicesPassingSelectionCuts = tfs->make<TH1D>("SuperTrueSlicesPassingSelectionCuts",        ";;True Slices Passing Selection Cuts",1, 0.0, 1.0);
    fTrueProtonEFromTrackLength = tfs->make<TH2D>("TrueProtonEFromTrackLength",        ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);
    fSuperTrueProtonEFromTrackLength = tfs->make<TH2D>("SuperTrueProtonEFromTrackLength",        ";Track Length (cm);True Proton Energy (GeV)",100, 0.0, 150.0,100,0.0,2.0);

    fTotalPOT = 0.0;



  }

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

  { 
    fND = false;

    novadaq::cnv::DetId detID = fgeom->DetId();
    if (detID == novadaq::cnv::kNEARDET) fND = true;

    return;
  }
  

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

  void ReMIdDedxRock::analyze(art::Event const& evt)
  {
    //Only written for ND

    if (!fND) return;
    
    // get a geometry handle   ************************************** 
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<calib::Calibrator> cal;
    //***************************************************************

    //Does this event pass spill cuts?
    art::Handle<sumdata::SpillData> spillPot;
    if(bt->HaveTruthInfo()){
      evt.getByLabel(fGeneratorLabel, spillPot);
    }
    else{
      evt.getByLabel(fNumiBeamLabel, spillPot);
    }

    if (spillPot.failedToGet()){
      mf::LogError("Susie") << "Spill Data not found for real data event";
      return;
    }

    //SL -- only cut on beam quality for data
    if (!bt->HaveTruthInfo()){
      if (std::abs(spillPot->deltaspilltimensec) > 0.5e9) return;
      if ((spillPot->spillpot)*1e12 < 2e12) return;
      if ((spillPot->hornI < -202)||(spillPot->hornI > -198)) return;
      if ((spillPot->posx < -2.00)||(spillPot->posx > 2.00)) return;
      if ((spillPot->posy < -2.00)||(spillPot->posy > 2.00)) return;
      if ((spillPot->widthx < 0.57)||(spillPot->widthx > 1.58)) return;
      if ((spillPot->widthy < 0.57)||(spillPot->widthy > 1.58)) return;
    }

    //Data quality spill cuts

    // Get the EventQuality information
    art::Handle<sumdata::EventQuality> spillQual;
    evt.getByLabel(fSpillQualLabel, spillQual);

    if(spillQual.failedToGet()){
      mf::LogError("CAFMaker") << "Spill EventQuality not found for real data event";
      return;
    }

    if (spillQual->fracdcm3hits > 0.45) return;
    if (spillQual->nmissingdcms > 0) return;


    //Filling spill pot after it passes the cuts
    float spillpot = spillPot->spillpot;
    if (!bt->HaveTruthInfo()) spillpot *= 1e12;
    fTotalPOT += spillpot;

    fEventsPassingSpillCuts->Fill(0.5);

    //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);
    // get  assosciations between slices and qepid
    // ecm art::FindManyP<qeef::QePId> fndmnyQepid(slicecol, evt, fQepidLabel);

    //if (evt.event() == 443183){
    // std::cout<<"Right event 2!"<<std::endl;
    //}

    // 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 Kalman tracks
      if(tracks.size() != 1){ continue; }

      //if ((evt.event() == 443183)&&(iSlice == 2)){
      //        std::cout<<"Right event and slice!"<<std::endl;
      //}

      //Require BOTH are 3D
      bool all3D = true;

      // get the pid that belongs to this recotrack
      art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);

      int bestidx = -1;
      double bestpid = -1;
      // loop over all of the tracks and find the most muon like one;
      for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
        if(!tracks[itrk]->Is3D()){ 
          all3D = false;
          continue; 
        }
        // get the pid that belongs to this track
        art::Ptr<remid::ReMId> pid = foids.at(itrk);
        if(pid->Value() > bestpid){
          bestpid = pid->Value();
          bestidx = itrk;
        }
      }//End of loop over tracks

      if (!all3D){ continue; }

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

      //ECM I don't think these make sense anymore
      fPassRockMuonCut = false;
      fPassVertX = true;
      fPassVertY = true;
      fPassVertZ = true;
      fPassProtonRemid = true;
      fPassDedxRatio = true;
      fPassEDiffZTest = true;


      //Slightly harsh ReMId cut on muon track
      if(bestidx == -1 || bestpid < 0.8){ continue; }

      if((bestidx != 0) && (bestidx !=1)){
        std::cout<<"BUMMER! A weird track id for muon that shouldn't happen!!"<<bestidx<<std::endl;
        continue;
      }

      int protonidx = -1;
      if (bestidx == 0){
        protonidx = 1;
      }
      if (bestidx == 1){
        protonidx = 0;
      }

     if(protonidx == -1){
        std::cout<<"BUMMER! A weird track id for proton that shouldn't happen!!"<<std::endl;
        continue;
      }     

     //Filling muon track length properties
     fzStart = tracks[bestidx]->Start().Z();
     fzEnd = tracks[bestidx]->Stop().Z();
     fTrackLength = tracks[bestidx]->TotalLength();

     // Track related Info
     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();
    

//"(abs(xStart)>195||abs(yStart)>195||zStart<10)&&abs(xEnd)<150&&abs(yEnd)<150&&zEnd>50&&zEnd<1200"
     if( ( abs(fxStart > 195) ||  abs(fyStart > 195) ||  abs(fzStart < 10) ) &&
         abs(fxEnd) < 150 && abs(fyEnd) < 150 && fzEnd >50 && fzEnd < 1275  
/*&& fzDir > 0.9 */){
       fPassRockMuonCut = true;
     }


     unsigned int ntraj = tracks[bestidx]->NTrajectoryPoints();
     // loop over trajectory points starting from the end
     for(int itrj = ntraj-1; itrj >= 0; --itrj){
       TVector3 trajpoint = tracks[bestidx]->TrajectoryPoint(itrj);
           
     } // 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;

     std::vector< std::vector<int> >   mcells;
     std::vector< std::vector<float> > madcs;
     std::vector< std::vector<float> > mpes;
     std::vector< std::vector<float> > mpecorrs;
         
     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;
         
         std::vector<int>   mCell;
         std::vector<float> mPE;
         std::vector<float> mADC;
         std::vector<float> mPECorr;

         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 minz = xyzplane[2]-dz;
         //      double maxz = xyzplane[2]+dz;
         //      int matchpt = 0;  // Not used for now
         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;
         /* Not used for now
         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;}
         }
         */

         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];
                 mCell.push_back(  chit->Plane() );
                 mPE.push_back(    chit->PE()    );
                 mADC.push_back(   chit->ADC()   );
                 mPECorr.push_back(rhit.PECorr() );

                 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];

                 mCell.push_back(  chit->Plane() );
                 mPE.push_back(    chit->PE()    );
                 mADC.push_back(   chit->ADC()   );
                 mPECorr.push_back(rhit.PECorr() );

                 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); 

           // 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; }
           
           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);
           
           mcells.push_back(mCell);
           madcs.push_back(mADC);
           mpes.push_back(mPE);
           mpecorrs.push_back(mPECorr);

           mcdes.push_back(mcPlaneE);
           mcdebs.push_back(mcPlaneEBirks);
           mcdxs.push_back(mcPathLength);

           } //nhits > 0         
     } // iPlane

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

     art::Ptr<remid::ReMId> pid = foids.at(bestidx);

     //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; }
     if (!(numuE.NDHadCalCat() + numuE.NDHadCalTran() < 0.03)){ continue; }
     //Not much hadronic energy contamination on muon track
     if (numuE.HadTrkE() < 0.025){ fPassHadEOnMuonTrack = true; }

     //Looking at qepid things //ecm
/*
     qeef::QePId qepid;
     if (!fndmnyQepid.isValid()){ continue; }
     std::vector<art::Ptr<qeef::QePId> > qVec = fndmnyQepid.at(iSlice);
     if (qVec.empty()){ continue; }
     qepid = *qVec[0];     

     //Require ratio of average dE/dx  track is good
     if (qepid.Dedx() > 1.75){ fPassDedxRatio = true; }
     //Require Ztest of two methods of calculating QE neutrino energy is good
     if (qepid.EdiffZ() > -1.3){ fPassEDiffZTest = true; }
*/

// ECM need to calculate qepid variables

//Require ratio of average dE/dx of two tracks is good  

     if(dedxs.size()>0){
       double dedx_avg_muon   = std::accumulate(dedxs.begin(),dedxs.end(),0.) / dedxs.size();
//       double dedx_avg_proton = std::accumulate(dedxp.begin(),dedxp.end(),0.) / dedxp.size();
//       std::cerr << std::accumulate(dedxs.begin(),dedxs.end(),0.) << " size " << dedxs.size() << std::endl;
//       std::cerr << std::accumulate(dedxp.begin(),dedxp.end(),0.) << " size " << dedxp.size() << std::endl;
//       std::cerr << "muon" << dedx_avg_muon << " proton " << dedx_avg_proton << " ratio " << dedx_avg_muon /dedx_avg_proton << std::endl;
//       if(dedx_avg_proton /  dedx_avg_muon > 1.75) fPassDedxRatio = true;

//ECM THIS MUST BE A MISTAKE
       if(dedx_avg_muon > 1.75) fPassDedxRatio = true;
     }

     //Require Ztest of two methods of calculating QE neutrino energy is good 
     // std::cerr << numuE.AngleQEE() << " " << numuE.TrkQEE() << " " << numuE.AngleQEError() << std::endl;
     //std::cerr << (numuE.AngleQEE()-numuE.TrkQEE())/numuE.AngleQEError() << std::endl;
     double ediffz = (numuE.AngleQEE()-numuE.TrkQEE())/numuE.AngleQEError();
        if (ediffz > -1.3){ fPassEDiffZTest = 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 two tracks
     int hitDiff = (slicelist[iSlice]->NCell())-(tracks[bestidx]->NCell());
     //if (!(hitDiff < 4)){ continue; }
     if (hitDiff < 4){ fPassOffTrackHits = true; }

     //Check if true signal for sim events
     fTrueNumuCCQE = false;
     fTrueNumuCCWithProton = false;
     fSuperTrueNumuCCWithProton = false;

     fpdg = 0;

     if (bt->HaveTruthInfo()){
       std::vector<cheat::NeutrinoEffPur> funcReturn = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),allhits);
       if (!funcReturn.empty()){
         //      const simb::MCNeutrino& test_neutrino = funcReturn[0].neutrinoInt->GetNeutrino();

         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();
           //      std::cout << "Mu " << fTrueE  << " " << fpdg << std::endl;
         }
       }
     }
     
//     if(ContainedEvent(slicelist[iSlice],tracks[bestidx])){
     if(fPassRockMuonCut )
//     if(1 )
       {   // Containment cut removed for now. Always pass //ecm try to undo 
       if (fPassOffTrackHits && fPassOffTrackEnergy && fPassHadEOnMuonTrack && fPassVertX && fPassVertY && fPassVertZ && fPassProtonRemid && fPassDedxRatio && fPassEDiffZTest){
         fSlicesPassingSelectionCuts->Fill(0.5);
/*
         std::cout << evt.run() << " " 
                   << evt.subRun() << " " 
                   << evt.event() <<  " "
                   << iSlice
                   << std::endl;
*/
         if (fTrueNumuCCWithProton){
           fTrueSlicesPassingSelectionCuts->Fill(0.5);
         }
       }
       // get all the dedx variables filled
       fnDedx = pid->NDedx()-1;
       fnDedxUsed = pid->NMeasurementPlanes();
       // loop over all of the dE/dx measurements
       for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
         fDedx = pid->Dedx(idedx);
         fVert = pid->DedxVertex(idedx);
         fUsed = pid->DedxUsed(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];

         fMultiCell   = mcells[idedx];
         fMultiADC    = madcs[idedx];
         fMultiPE     = mpes[idedx];
         fMultiPECorr = mpecorrs[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();
         // fill the dedx sample tree
         DedxSample->Fill();
       }//Loop over DedxSample       
     }//If passes contaiment         
    } // end loop over slices        

  }


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

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


  //.......................................................................
  bool ReMIdDedxRock::ContainedEvent(art::Ptr<rb::Cluster> slice, art::Ptr<rb::Track> trk)
  {
    
    bool contained = false;
    
    // check that this is a contained event
    
    // get some useful information here
    // 1) how close is the slice to the edge of the detector
    unsigned int slcncellsfromedge = 999999999;
    for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
      const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
      const int plane = chit->Plane();
      unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
      if(cellsfromedge < slcncellsfromedge){
        slcncellsfromedge = cellsfromedge;
       }
    }
    // 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)&&(slcncellsfromedge > 1)&&(slcTime>50000)&&(slcTime<450000)&&(trackNHit > 5)){
      if(firstpl > 1 && lastpl < 212){
        //std::cout<<"In contained event loop! 1"<<std::endl;
        int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
        int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
        // also make sure the track doesn't start too close to the end of the active region
        if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
          // if the track ends in the muon catcher, check that it doesn't clip the corner
          //std::cout<<"In contained event loop! 2"<<std::endl;
          if(trk->Stop().Z() < 1275){ 
            //std::cout<<"In contained event loop! 2a"<<std::endl;
            contained = true; 
          }
          else{
            //std::cout<<"In contained event loop! 2b"<<std::endl;
            /*
            // find the trajectory point closest in z to the transition plane
            double mindist = 99999999;
            unsigned int mintrajid = 0;
            for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
              double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
              if(dist > 0 && dist < mindist){
                mindist = dist;
                mintrajid = itraj;
              }
            }
            // get the direction at this trajectory point
            TVector3 trandir = trk->Dir();
            if(mintrajid < trk->NTrajectoryPoints() -2){
              TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
            }
            */
            //std::cout<<"In contained event loop! 3"<<std::endl;
            //double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
            double ypostranKirk = flivegeom->YPositionAtTransition(trk->Start(),trk->Dir());
            if( ypostranKirk < 55){ 
              //std::cout<<"In contained event loop! 4"<<std::endl;
              contained = true; 
            }
            
          } // end else in muon catcher
        } // end if track ends projections contained
      } // end if track end planes contained 
        // } // end if  expected number of diblocks
    } // end if slice contained
    
    
    return contained;
    
  }
  
 


}


namespace remid{

  DEFINE_ART_MODULE(ReMIdDedxRock)

}