NueSel_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   Brief please!
/// \author  Jianming Bian - bianjm@physics.umn.edu
////////////////////////////////////////////////////////////////////////
#include "RecoJMShower/JMShower.h"
#include "RecoJMShower/EID.h"

extern "C" {
#include <sys/types.h>
#include <sys/stat.h>
}
#include <cmath>
#include <iostream>
#include <fstream>
#include <vector>

// ROOT includes
#include "TFile.h"
#include "TH1D.h"
#include "TLorentzVector.h"
#include "TMath.h"
#include "TMultiLayerPerceptron.h"
#include "TRandom.h"
#include "TTree.h"

// NOvA includes
#include "DAQChannelMap/DAQChannelMapConstants.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/FilterList.h"
//#include "RecoBase/Vertex.h"
#include "RecoBase/RecoHit.h"
#include "Utilities/AssociationUtil.h"
#include "Calibrator/Calibrator.h"

// ART includes
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Principal/Handle.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "fhiclcpp/ParameterSet.h"
//#include "messagefacility/MessageLogger/MessageLogger.h"


namespace jmshower
{

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

    virtual void produce(art::Event& evt);

    virtual void reconfigure(const fhicl::ParameterSet& pset);

    
    void beginJob();
    void endJob();
    void beginRun(art::Run& run);
    double EnergyEstimator(TVector3 beamDir, TLorentzVector elecP4, double etot, int mode);
    
  protected:

    std::string fJMShowerLabel;  ///< Label for jmshower::JMShower input
    std::string fSlicerLabel;
    std::string fLibPath;
    // Osc. weight histograms
    TTree* fEvent;
    //True Nuetrino info
    int    m_evtTrueNuCCNC;
    int    m_evtTrueNuMode;
    int    m_evtTrueNuPdg;
    double m_evtTrueNuEnergy;
    int    m_evtTrueNuIsFid;
    double m_evtTrueNuP4[4];
    double m_evtTrueNuVtx[3];
    int    m_evtTruePdg[20];
    int    m_evtTrueMother[20];
    double m_evtTrueEnergy[20];
    //Reco info
    int  m_evtRun;
    int  m_evtSRun;
    int  m_evtEvent;
    double m_evtNSlice;
    double m_evtNSh;
    double m_evt;
    double m_evtANN;
    double m_evtSh1DedxANN;
    double m_evtSh1DedxANNE;
    double  m_evtEtot;
    double m_evtSigOscW;
    double m_evtBkgOscW;
    int m_evtIsMuon;
    int m_evtIsInvPhoton;
    int m_evtIsMuonSh1;
    double m_evtDetEnergy;
    double m_evtNuEnergy;
    double m_evtSh1Energy;
    double m_evtSh1ExclEnergy;
    double m_evtSh1Length;
    double m_evtSh1Start[3];
    double m_evtSh1Stop[3];
    double m_evtSh1Dir[3];
    double m_evtVtx[3];
    
    double  m_evtSh1EGLLL;
    double  m_evtSh1EGLLT;
    double  m_evtSh1EMuLLL;
    double  m_evtSh1EMuLLT;
    double  m_evtSh1EPi0LLL;
    double  m_evtSh1EPi0LLT;
    double  m_evtSh1EPLLL;
    double  m_evtSh1EPLLT;
    double  m_evtSh1ENLLL;
    double  m_evtSh1ENLLT;
    double  m_evtSh1EPiLLL; 
    double  m_evtSh1EPiLLT;
    double  m_evtSh1VtxGeV; 
    double  m_evtSh1Pi0Mass;
    double  m_evtSh1ShE;
    double  m_evtSh1Gap;
    double  m_evtSh1ELLL;
    double  m_evtSh1ELLT;
    double  m_evtSh1MuLLL;
    double  m_evtSh1MuLLT;
    double  m_evtSh1GLLL;
    double  m_evtSh1InvGLLL;

    double m_evtSh1DistToEdge;
    double m_evtSh1SliceEdge2GeV;
    double m_evtSh1SliceEdge5GeV;
    double m_evtSh1SliceEdge10GeV;
    
  private:
    
    TH1D *ht_nue_rw;
    TH1D *ht_bkg_nue_rw;
    TH1D *ht_bkg_numu_rw;
    
    TMultiLayerPerceptron *mlp_evt;
    
    
    double OscWeight(int pdg, int ccnc, double energy, bool issig) ;
    bool LoadHistograms(int detid);
    void InitializeNTuple();
  };
} // namespace


namespace jmshower 
{
  int Ncount0 = 0, Ncount1 = 0, Ncount2 = 0;
  int Nfill = 0, Nput = 0;
  // mass for electron, photon, muon, pi0, NA, proton, neutron and charged 
  const double xmass[8] = {0.000511, 0, 0.105658, 0.1349766, 0, 0.938272, 0.939565379, 0.139570};
  //......................................................................
  NueSel::NueSel(const fhicl::ParameterSet& pset)
    : fEvent(0)
  {
    reconfigure(pset);

//    produces<std::vector<rb::Vertex> >();
    produces<std::vector<jmshower::EID> >();
    produces< art::Assns<jmshower::EID, rb::Cluster> >();
  }

  //......................................................................
  NueSel::~NueSel()
  {
  }

  //......................................................................
  void NueSel::reconfigure(const fhicl::ParameterSet& pset)
  {
    fJMShowerLabel = pset.get< std::string >("JMShowerLabel");
    fLibPath       = pset.get< std::string >("LibraryPath");
    fSlicerLabel   = pset.get< std::string >("SlicerLabel");
  }

  bool NueSel::LoadHistograms(int detid)
  {
    std::string fnameOs(fLibPath);
    switch (detid) {
      case novadaq::cnv::kFARDET:
        fnameOs += "extrahist_fd.root";
      break;
      case novadaq::cnv::kNDOS:
        fnameOs += "extrahist_ndos.root";
      break;
      case novadaq::cnv::kNEARDET:
        fnameOs += "extrahist_ndos.root";
      break;
      default:
        throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
                                        << " Bad detector id = " << detid;
    }
/*
    std::string fFnameOs; 
    cet::search_path sp("FW_SEARCH_PATH");
    sp.find_file(fnameOs, fFnameOs);
    struct stat sb;
    if (fFnameOs.empty() || stat(fFnameOs.c_str(), &sb)!=0){
      // Failed to resolve the file name
      mf::LogWarning("NueSel") << "Histogram -- Os " << fFnameOs << " not found.";
      throw cet::exception("FILE_NOT_FOUND") << "Histogram file not found.";
    }
*/
    TFile hfileOs(fnameOs.c_str());
    ht_nue_rw = (TH1D*)hfileOs.Get("ht_nue_rw");
    ht_bkg_nue_rw = (TH1D*)hfileOs.Get("ht_bkg_nue_rw");
    ht_bkg_numu_rw = (TH1D*)hfileOs.Get("ht_bkg_numu_rw");

    std::string fnameWeightsEvt(fLibPath);
    switch (detid) {
      case novadaq::cnv::kFARDET:
        fnameWeightsEvt += "weights_evt_fd.txt";
      break;
      case novadaq::cnv::kNDOS:
        fnameWeightsEvt += "weights_evt_ndos.txt";
      break;
      case novadaq::cnv::kNEARDET:
        fnameWeightsEvt += "weights_evt_ndos.txt";
      break;
      default:
        throw cet::exception("BAD DETID") << __FILE__ << ":" << __LINE__
                                        << " Bad detector id = " << detid;
    }

//    std::cout<<"Trained Evt ANN "<<fnameWeightsEvt.c_str()<<" is loaded."<<std::endl;
    double egLLL;
    double egLLT;
    double epi0LLL;
    double epi0LLT;
    double epLLL;
    double epLLT;
    double enLLL;
    double enLLT;
    double epiLLL;
    double epiLLT;
    double pi0mass;
    double vtxgev;
    int type;
    TTree *trainTree_evt = new TTree("trainTree_evt","m_trainTree_evt");
    trainTree_evt->Branch("egLLL", &egLLL, "egLLL/D");
    trainTree_evt->Branch("egLLT", &egLLT, "egLLT/D");
    trainTree_evt->Branch("epi0LLL", &epi0LLL, "epi0LLL/D");
    trainTree_evt->Branch("epi0LLT", &epi0LLT, "epi0LLT/D");
    trainTree_evt->Branch("epLLL", &epLLL, "epLLL/D");
    trainTree_evt->Branch("epLLT", &epLLT, "epLLT/D");
    trainTree_evt->Branch("enLLL", &enLLL, "enLLL/D");
    trainTree_evt->Branch("enLLT", &enLLT, "enLLT/D");
    trainTree_evt->Branch("epiLLL", &epiLLL, "epiLLL/D");
    trainTree_evt->Branch("epiLLT", &epiLLT, "epiLLT/D");
    trainTree_evt->Branch("pi0mass", &pi0mass, "pi0mass/D");
    trainTree_evt->Branch("vtxgev", &vtxgev, "vtxgev/D");
    trainTree_evt->Branch("type", &type, "type/I");
    mlp_evt = new TMultiLayerPerceptron("egLLL,egLLT,epi0LLL,epi0LLT,epLLL,epLLT,enLLL,enLLT,epiLLL,epiLLT,pi0mass,vtxgev:18:12:6:type",trainTree_evt);
    mlp_evt->LoadWeights(fnameWeightsEvt.c_str());

    return true;
  }


  //......................................................................
  void NueSel::beginJob()
  {
std::cout<<"begin 1 "<<std::endl;
    if(fEvent) return; // Don't do this twice
std::cout<<"begin 2 "<<std::endl;
      art::ServiceHandle<art::TFileService> tfs;

std::cout<<"begin 3 "<<std::endl;
//True neutrino information
    fEvent = tfs->make<TTree>("fEvent", "fEvent");
      fEvent->Branch("evtTrueNuCCNC", &m_evtTrueNuCCNC);
      fEvent->Branch("evtTrueNuMode", &m_evtTrueNuMode);
      fEvent->Branch("evtTrueNuPdg", &m_evtTrueNuPdg);
      fEvent->Branch("evtTrueNuEnergy", &m_evtTrueNuEnergy);
      fEvent->Branch("evtTrueNuIsFid", &m_evtTrueNuIsFid);
      fEvent->Branch("evtTrueNuP4[4]", m_evtTrueNuP4);
      fEvent->Branch("evtTrueNuVtx[3]", m_evtTrueNuVtx);
      fEvent->Branch("evtTrueEnergy[20]", m_evtTrueEnergy);
      fEvent->Branch("evtTruePdg[20]", m_evtTruePdg);
      fEvent->Branch("evtTrueMother[20]", m_evtTrueMother);

//Reconstructed information
      fEvent->Branch("evtRun", &m_evtRun);
      fEvent->Branch("evtSRun", &m_evtSRun);
      fEvent->Branch("evtNSlice", &m_evtNSlice);
      fEvent->Branch("evtNSh", &m_evtNSh);
      fEvent->Branch("evtANN", &m_evtANN);
      fEvent->Branch("evtANNSh1", &m_evtSh1DedxANN);
      fEvent->Branch("evtANNESh1", &m_evtSh1DedxANNE);
      fEvent->Branch("evtEtot", &m_evtEtot);
      fEvent->Branch("evtSigOscW", &m_evtSigOscW);
      fEvent->Branch("evtBkgOscW", &m_evtBkgOscW);
      fEvent->Branch("evtIsMuon", &m_evtIsMuon);
      fEvent->Branch("evtIsMuonSh1", &m_evtIsMuonSh1);
      fEvent->Branch("evtIsInvPhoton", &m_evtIsInvPhoton);
      fEvent->Branch("evtDetEnergy", &m_evtDetEnergy);
      fEvent->Branch("evtNuEnergy", &m_evtNuEnergy);
      fEvent->Branch("evtSh1Energy", &m_evtSh1Energy);
      fEvent->Branch("evtSh1ExclEnergy", &m_evtSh1ExclEnergy);
      fEvent->Branch("evtSh1Length", &m_evtSh1Length);
      fEvent->Branch("evtSh1Start[3]", m_evtSh1Start);
      fEvent->Branch("evtSh1Stop[3]", m_evtSh1Stop);
      fEvent->Branch("evtSh1Dir[3]", m_evtSh1Dir);

      fEvent->Branch("evtVtx[3]", m_evtVtx);
      fEvent->Branch("evtSh1EGLLL",&m_evtSh1EGLLL); 
      fEvent->Branch("evtSh1EGLLT",&m_evtSh1EGLLT);
      fEvent->Branch("evtSh1EMuLLL",&m_evtSh1EMuLLL);
      fEvent->Branch("evtSh1EMuLLT",&m_evtSh1EMuLLT);
      fEvent->Branch("evtSh1EPi0LLL",&m_evtSh1EPi0LLL);
      fEvent->Branch("evtSh1EPi0LLT",&m_evtSh1EPi0LLT);
      fEvent->Branch("evtSh1EPLLL",&m_evtSh1EPLLL);
      fEvent->Branch("evtSh1EPLLT",&m_evtSh1EPLLT);
      fEvent->Branch("evtSh1ENLLL",&m_evtSh1ENLLL);
      fEvent->Branch("evtSh1ENLLT",&m_evtSh1ENLLT);
      fEvent->Branch("evtSh1EPiLLL",&m_evtSh1EPiLLL);
      fEvent->Branch("evtSh1EPiLLT",&m_evtSh1EPiLLT);
      fEvent->Branch("evtSh1VtxGeV",&m_evtSh1VtxGeV);
      fEvent->Branch("evtSh1Pi0Mass",&m_evtSh1Pi0Mass);
      fEvent->Branch("evtSh1ShE",&m_evtSh1ShE);
      fEvent->Branch("evtSh1Gap",&m_evtSh1Gap);
      fEvent->Branch("evtSh1ELLL",&m_evtSh1ELLL);
      fEvent->Branch("evtSh1ELLT",&m_evtSh1ELLT);
      fEvent->Branch("evtSh1MuLLL",&m_evtSh1MuLLL);
      fEvent->Branch("evtSh1MuLLT",&m_evtSh1MuLLT);
      fEvent->Branch("evtSh1GLLL",&m_evtSh1GLLL);
      fEvent->Branch("evtSh1InvGLLL",&m_evtSh1InvGLLL);
std::cout<<"begin 4 "<<std::endl;
      
      fEvent->Branch("evtSh1DistToEdge", &m_evtSh1DistToEdge);
      fEvent->Branch("evtSh1SliceEdge2GeV", &m_evtSh1SliceEdge2GeV);
      fEvent->Branch("evtSh1SliceEdge5GeV", &m_evtSh1SliceEdge5GeV);
      fEvent->Branch("evtSh1SliceEdge10GeV", &m_evtSh1SliceEdge10GeV);
std::cout<<"begin 5 "<<std::endl;
  }
  void NueSel::beginRun(art::Run& /*run*/)
  {
std::cout<<"begin run "<<std::endl;
      art::ServiceHandle<geo::Geometry> geom;
      NueSel::LoadHistograms(geom->DetId());
std::cout<<"after begin run "<<std::endl;
  }

  void NueSel::InitializeNTuple()
  {
//fEvent
     m_evtTrueNuCCNC = 0;
     m_evtTrueNuMode = 0;
     m_evtTrueNuPdg = 0;
     m_evtTrueNuEnergy = 0;
     m_evtTrueNuIsFid = 0;
     for(int i=0;i<4;i++){m_evtTrueNuP4[i] = 0;} 
     for(int i=0;i<3;i++){m_evtTrueNuVtx[i] = 0;}
     for(int i=0;i<20;i++){m_evtTruePdg[i] = 0;} 
     for(int i=0;i<20;i++){m_evtTrueMother[i] = 0;} 
     for(int i=0;i<20;i++){m_evtTrueEnergy[i] = 0;} 
     m_evtRun = 0;
     m_evtSRun = 0;
     m_evtEvent = 0;
     m_evtNSlice = 0;
     m_evtNSh = 0;
     m_evt = 0;
     m_evtANN = -5;
     m_evtSh1DedxANN = -5;
     m_evtSh1DedxANNE = -5;
     m_evtEtot = 0;
     m_evtSigOscW = 0;
     m_evtBkgOscW = 0;
     m_evtIsMuon = 0;
     m_evtIsInvPhoton = 0;
     m_evtDetEnergy = 0;
     m_evtTrueNuEnergy = 0;
     m_evtNuEnergy = 0;
     m_evtSh1Energy = 0;
     m_evtSh1ExclEnergy = 0;
     m_evtSh1Length = 0;

     for(int i=0;i<3;i++){ 
       m_evtVtx[i]=0;
       m_evtSh1Start[i]=0;
       m_evtSh1Stop[i]=0;
       m_evtSh1Dir[i]=0;
     }

     m_evtSh1EGLLL = 0;
     m_evtSh1EGLLT = 0;
     m_evtSh1EMuLLL = 0;
     m_evtSh1EMuLLT = 0;
     m_evtSh1EPi0LLL = 0;
     m_evtSh1EPi0LLT = 0;
     m_evtSh1EPLLL = 0;
     m_evtSh1EPLLT = 0;
     m_evtSh1ENLLL = 0;
     m_evtSh1ENLLT = 0;
     m_evtSh1EPiLLL = 0;
     m_evtSh1EPiLLT = 0;
     m_evtSh1VtxGeV = 0;
     m_evtSh1Pi0Mass = 0;
     m_evtSh1ShE = 0;
     m_evtSh1Gap = 0;
     m_evtSh1ELLL = 0;
     m_evtSh1ELLT = 0;
     m_evtSh1MuLLL = 0;
     m_evtSh1MuLLT = 0;
     m_evtSh1GLLL = 0;
     m_evtSh1InvGLLL = 0;
     m_evtSh1DistToEdge = 0;
     m_evtSh1SliceEdge2GeV = 0; 
     m_evtSh1SliceEdge5GeV = 0; 
     m_evtSh1SliceEdge10GeV = 0; 
  }

  //......................................................................
  double NueSel::OscWeight(int pdg, int ccnc, double energy, bool issig) 
  {// Obtain oscillation weight for swapping MC
    double weight = 1.0;
    if(fabs(pdg)==12&&ccnc==0&&issig==true){
       weight = ht_nue_rw->Interpolate(energy); 
    }
    if(fabs(pdg)==14&&ccnc==0&&issig==false) weight = 0.1787396954; 
    return weight;
  }
  //............................
  double NueSel::EnergyEstimator(TVector3 beamDir, TLorentzVector elecP4, double etot, int mode)
  {// Neutrino energy estimator
    double Ee = elecP4.E();
    double Ehad = etot - Ee;
    double EhadCor = 0.282525 + 1.0766*Ehad;
    if(mode!=0&&mode!=1) return etot;
  // EM shower has been calibrated in RecoJMShower, so there is no correction for Ee.
  // The overall correction is for the assymetric convoluton effect. 
    double Enue = (Ee + EhadCor)/0.98343; 
    return Enue;
  }


  //......................................................................
  void NueSel::endJob()
  {
    std::cout<<"Total number of events   "<<Ncount0<<std::endl;
    std::cout<<"Muon like events         "<<Ncount1<<std::endl;
    std::cout<<"Nue CC events            "<<Ncount2<<std::endl;
    std::cout<<"Nfill, Nput            "<<Nfill<<" "<<Nput<<std::endl;
  }

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

  void NueSel::produce(art::Event& evt)
  {
//----------------------------------------------------
//  Run, subrun and event information
//----------------------------------------------------
    Ncount0++;
    int run   = evt.run();
    int srun  = evt.subRun();
    int event = evt.id().event();
//    std::cout<<"In NueSel event "<<event<<std::endl;


// Neutrino MC truth
    TLorentzVector evtTrueNuP4(0,0,0,0);
    TVector3 evtTrueNuVtx(0,0,0);
    int evtTrueNuCCNC = -1;
    int evtTrueNuMode = -1;
    int evtTrueNuPdg = 0;
    double evtTrueNuEnergy = 0;

/*
    if(!evt.isRealData()){
      for(int i=0;i<1;i++){
       art::Handle< std::vector<simb::MCTruth> > mc;
       try {
         evt.getByLabel("generator", mc);
         if (mc->size()!=1) break;
       }
       catch (...) {
         break;
       }
       if ((*mc)[0].NeutrinoSet()) {
         const simb::MCNeutrino& gen = (*mc)[0].GetNeutrino();
         const simb::MCParticle& nu  = gen.Nu();
         evtTrueNuP4 = nu.Momentum();
         evtTrueNuVtx.SetXYZ(nu.EndX(),nu.EndY(),nu.EndZ());
         evtTrueNuCCNC = gen.CCNC() ;
         evtTrueNuMode = gen.Mode() ;
         evtTrueNuPdg = nu.PdgCode() ;
         evtTrueNuEnergy = nu.Momentum().E() ;
       }
      }
    }
*/


// Showers reconstruction

    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<calib::Calibrator> cal;


// Event information with PID
    std::unique_ptr<std::vector<jmshower::EID> > eidcol(new std::vector<jmshower::EID>);
    std::unique_ptr< art::Assns<jmshower::EID, rb::Cluster> > assns(new art::Assns<jmshower::EID, rb::Cluster>);


//  EID product for analysis.
    art::Handle<std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fSlicerLabel, slicecol);

    for(unsigned int sliceIdx = 0; sliceIdx < slicecol->size(); ++sliceIdx){//Reconstruct EID slice by slice
      if(rb::IsFiltered(evt,slicecol,sliceIdx)) continue;
      const rb::Cluster& slice = (*slicecol)[sliceIdx];
      if(slice.IsNoise()) continue;

      double evtGapTNS = 1e8;
      for(unsigned int iic = 0; iic < slicecol->size(); ++iic){
        if(iic==sliceIdx)continue;
        const rb::Cluster& slice1 = (*slicecol)[iic];
        if(slice1.IsNoise()) continue;
//        if(slice1.NCell()<20)continue;
        if(evtGapTNS > fabs(slice1.MeanTNS()-slice.MeanTNS())){
          evtGapTNS = fabs(slice1.MeanTNS()-slice.MeanTNS());
        } 
      }


      NueSel::InitializeNTuple();
      m_evtRun = run;
      m_evtSRun = srun;
      m_evtEvent = event;

      m_evtTrueNuCCNC = evtTrueNuCCNC;
      m_evtTrueNuMode = evtTrueNuMode;
      m_evtTrueNuPdg = evtTrueNuPdg;


      m_evtTrueNuEnergy = evtTrueNuEnergy;
      m_evtTrueNuP4[0] = evtTrueNuP4.Px();
      m_evtTrueNuP4[1] = evtTrueNuP4.Py();
      m_evtTrueNuP4[2] = evtTrueNuP4.Pz();
      m_evtTrueNuP4[3] = evtTrueNuP4.E();
      m_evtTrueNuVtx[0] = evtTrueNuVtx[0];
      m_evtTrueNuVtx[1] = evtTrueNuVtx[1];    
      m_evtTrueNuVtx[2] = evtTrueNuVtx[2];

      const double wx = (slice.NXCell() > 0) ? slice.W(slice.XCell(0).get()) : 0;
      const double wy = (slice.NYCell() > 0) ? slice.W(slice.YCell(0).get()) : 0;
      //calculate slce reco energy
      double slicegev = 0;
      for(unsigned int i=0;i<slice.NCell();i++){
        const art::Ptr<rb::CellHit>& chit = slice.Cell(i);
        double w = wx;
        if(geom->Plane(chit->Plane())->View()==geo::kY)w=wy;
//        const rb::RecoHit rhit = cal->MakeRecoHit(chit.get(), w);
        const rb::RecoHit rhit(cal->MakeRecoHit(*chit, w ));


        double gev=0;
        if(rhit.IsCalibrated())gev = rhit.GeV();
        slicegev+=gev;        
      }

      art::FindManyP<jmshower::JMShower> jms(slicecol, evt, fJMShowerLabel);
      std::vector<art::Ptr<jmshower::JMShower> >  showercol = jms.at(sliceIdx);
      int evtNSh = showercol.size();

      jmshower::EID eid;
      if(evtNSh<1){ // No shower be recontructed
        eid.SetANN(-5.);
        eid.SetANNE(-5.);
        eid.SetANNEL(-5.);
        eid.SetANNEPi0(-5.);
        eid.SetANNEPi0EL(-5.);
        eid.SetANNNoCos(-5.);
        eid.SetANNENoCos(-5.);
        eid.SetVtx(evtTrueNuVtx);
        eid.SetNuEnergy(slicegev*1.9);
        eid.SetDetEnergy(0);
        eid.SetEvtNCell(slice.NCell());
        eid.SetEvtSumCosTheta(0);
        eid.SetEvtSumP(0);
        eid.SetEvtSumPt(0);
        eid.SetEvtSumP0(0);
        eid.SetEvtSumPt0(0);

        eid.SetEvtEtot(0);
        eid.SetEvtMinX(99999);
        eid.SetEvtMinY(99999);
        eid.SetEvtMinZ(99999);
        eid.SetEvtMaxX(-99999);
        eid.SetEvtMaxY(-99999);
        eid.SetEvtMaxZ(-99999);
        eid.SetEvtGapTNS(evtGapTNS);
        eid.SetEvtNCellToEdge(10000);

        eid.SetSh1Energy(0);
        eid.SetSh1ExclEnergy(0);
        eid.SetSh1TotalLength(0);
        TVector3 v3null(0,0,0);
        eid.SetSh1Start(v3null);
        eid.SetSh1Stop(v3null);
        eid.SetSh1Dir(v3null);
        eid.SetSh1VtxDoca(0);
        eid.SetSh1Gap(0);
        eid.SetSh1NPlane(0);
        eid.SetSh1XNPlane(0);
        eid.SetSh1YNPlane(0);
        eid.SetSh1NCell(0);
        eid.SetSh1XNCell(0);
        eid.SetSh1YNCell(0);
        eid.SetSh1PID(-1);
        eid.SetSh1NCellToEdge(10000);

        eid.SetSh1Sh2Dang(-5);
        eid.SetSh2Energy(0);
        eid.SetSh2ExclEnergy(0);
        eid.SetSh2TotalLength(0);
        eid.SetSh2Start(v3null);
        eid.SetSh2Stop(v3null);
        eid.SetSh2Dir(v3null);
        eid.SetSh2VtxDoca(0);
        eid.SetSh2Gap(0);
        eid.SetSh2NPlane(0);
        eid.SetSh2XNPlane(0);
        eid.SetSh2YNPlane(0);
        eid.SetSh2NCell(0);
        eid.SetSh2XNCell(0);
        eid.SetSh2YNCell(0);
        eid.SetSh2PID(-1);
        eid.SetSh2NCellToEdge(10000);


        eid.SetIsMuon(0);
        eid.SetIsInvPhoton(0);
        eid.SetSigOscW(0);
        eid.SetBkgOscW(0);

        eid.SetDedx0(-5);
        eid.SetDedx1(-5);
        eid.SetDedx2(-5);
        eid.SetDedx3(-5);
        eid.SetDedx4(-5);
        eid.SetDedx5(-5);

        eid.SetEGLLL(-10);
        eid.SetEGLLT(-10);
        eid.SetEMuLLL(-10);
        eid.SetEMuLLT(-10);
        eid.SetEPi0LLL(-10);
        eid.SetEPi0LLT(-10);
        eid.SetEPLLL(-10);
        eid.SetEPLLT(-10);
        eid.SetENLLL(-10);
        eid.SetENLLT(-10);
        eid.SetEPiLLL(-10);
        eid.SetEPiLLT(-10);
        eid.SetVtxGeV(0);
        eid.SetPi0Mass(0);
        eid.SetShE(0);
        eid.SetGap(0);
        eid.SetELLL(-10);
        eid.SetELLT(-10);
        eid.SetEelLLL(-10);
        eid.SetEelLLT(-10);
        eid.SetMuLLL(-10);
        eid.SetMuLLT(-10);
        eid.SetGLLL(-10);
        eid.SetInvGLLL(-10);

        eid.SetSh1DistToEdge(0);
        eid.SetEEdge2Cell(0);
        eid.SetEEdge5Cell(0);
        eid.SetEEdge10Cell(0);


        eidcol->push_back(eid);
        Nfill++;
        fEvent->Fill(); // Fill NTuple

      }else{

        double evtMinX = 99999;
        double evtMinY = 99999;
        double evtMinZ = 99999;
        double evtMaxX = -99999;
        double evtMaxY = -99999;
        double evtMaxZ = -99999;
        for(unsigned int i=0;i<showercol.size();i++){

          if((showercol[i]->Start())[0]<evtMinX)evtMinX = (showercol[i]->Start())[0];
          if((showercol[i]->Stop())[0]<evtMinX)evtMinX = (showercol[i]->Stop())[0];
          if((showercol[i]->Start())[1]<evtMinY)evtMinY = (showercol[i]->Start())[1];
          if((showercol[i]->Stop())[1]<evtMinY)evtMinY = (showercol[i]->Stop())[1];
          if((showercol[i]->Start())[2]<evtMinZ)evtMinZ = (showercol[i]->Start())[2];
          if((showercol[i]->Stop())[2]<evtMinZ)evtMinZ = (showercol[i]->Stop())[2];

          if((showercol[i]->Start())[0]>evtMaxX)evtMaxX = (showercol[i]->Start())[0];
          if((showercol[i]->Stop())[0]>evtMaxX)evtMaxX = (showercol[i]->Stop())[0];
          if((showercol[i]->Start())[1]>evtMaxY)evtMaxY = (showercol[i]->Start())[1];
          if((showercol[i]->Stop())[1]>evtMaxY)evtMaxY = (showercol[i]->Stop())[1];
          if((showercol[i]->Start())[2]>evtMaxZ)evtMaxZ = (showercol[i]->Start())[2];
          if((showercol[i]->Stop())[2]>evtMaxZ)evtMaxZ = (showercol[i]->Stop())[2];
       }


        std::vector<TLorentzVector> showerpGam;// Set shower momentum with photon assumption
        showerpGam.clear();

        for(unsigned int i = 0; i < showercol.size(); i++){
//          std::cout<<"ish, N hit, energy "<<i<<" "<<showercol[i]->NCell()<<" "<<showercol[i]->Energy()<<std::endl;
          double showereraw = showercol[i]->Energy();;
          TLorentzVector showerptrk(showercol[i]->Dir()[0]*showereraw,showercol[i]->Dir()[1]*showereraw,showercol[i]->Dir()[2]*showereraw,showereraw);
          showerpGam.push_back(showerptrk);

        }
//       Looping over showers
        int ish1 = 0;
        int ish1e = 0;
        double sh1ann = -1;
        double sh1eann = -1;
        double sh1eanne = -1;
        double sh1eannnocos = -1;
        double sh1eannenocos = -1;
        double sh1gev = 0;
        double sh1egev = 0;
        double sh1elength = 0;
        double evtetot = 0;
        std::vector<TLorentzVector> showerP4Col;
        std::vector<TLorentzVector> showerP40Col;
        std::vector<int>  showerpid;
        showerP4Col.clear();
        showerP40Col.clear();

        double sh1egLLL = -10;  
        double sh1egLLT = -10;
        double sh1emuLLL = -10;
        double sh1emuLLT = -10;
        double sh1epi0LLL = -10;
        double sh1epi0LLT = -10;
        double sh1epLLL = -10;
        double sh1epLLT = -10;
        double sh1enLLL = -10;
        double sh1enLLT = -10;
        double sh1epiLLL = -10;
        double sh1epiLLT = -10;
        double sh1vtxgev = 0;
        double sh1pi0mass = 0;
        double sh1shE = 0;
        double sh1gap = 0;
        double sh1eLLL = -10;
        double sh1eLLT = -10;
        double sh1muLLL = -10;
        double sh1muLLT = -10;
        double sh1eelLLL = -10;
        double sh1eelLLT = -10;

 
        int evtNCellToEdge = 10000;

        for(unsigned int i = 0; i < showercol.size(); i++){
          if(evtNCellToEdge>showercol[i]->NCellToEdge())evtNCellToEdge=showercol[i]->NCellToEdge();
          double showerenergy = showercol[i]->Energy();
          double showerlength = showercol[i]->TotalLength();
//        ANN deterimantion
          double showerann = showercol[i]->ANN(0);
          double showeranne = showercol[i]->ANN(2);
          double showerannnocos = showercol[i]->ANN(3);
          double showerannenocos = showercol[i]->ANN(4);

          double ellmax = -9999;
          double elllmax = -9999;
          double elltmax = -9999;
          for(int itype = int(jmshower::kNull)+1; itype != int(jmshower::kLast); ++itype){
            if(itype == int(jmshower::kElectron)||itype == int(jmshower::kElectronQE)||itype == int(jmshower::kElectronRES)||itype == int(jmshower::kElectronDIS)||itype == int(jmshower::kElectronCOH)){
              if(ellmax<showercol[i]->DedxLongLL(itype)+showercol[i]->DedxTransLL(itype)){
                ellmax = showercol[i]->DedxLongLL(itype)+showercol[i]->DedxTransLL(itype);
                elllmax = showercol[i]->DedxLongLL(itype);
                elltmax = showercol[i]->DedxTransLL(itype);
              }
            }
          }

          double egLLL = elllmax - showercol[i]->DedxLongLL(int(jmshower::kGamma));
          double egLLT = elltmax - showercol[i]->DedxTransLL(int(jmshower::kGamma));
          double emuLLL = elllmax - showercol[i]->DedxLongLL(int(jmshower::kMuon));
          double emuLLT = elltmax - showercol[i]->DedxTransLL(int(jmshower::kMuon));
          double epi0LLL = elllmax - showercol[i]->DedxLongLL(int(jmshower::kPi0));
          double epi0LLT = elltmax - showercol[i]->DedxTransLL(int(jmshower::kPi0));
          double epLLL = elllmax - showercol[i]->DedxLongLL(int(jmshower::kProton));
          double epLLT = elltmax - showercol[i]->DedxTransLL(int(jmshower::kProton));
          double enLLL = elllmax - showercol[i]->DedxLongLL(int(jmshower::kNeutron));
          double enLLT = elltmax - showercol[i]->DedxTransLL(int(jmshower::kNeutron));
          double epiLLL = elllmax - showercol[i]->DedxLongLL(int(jmshower::kPion));
          double epiLLT = elltmax - showercol[i]->DedxTransLL(int(jmshower::kPion));




          double vtxgev = showercol[i]->VtxGeV();
          double pi0mass = std::max(showercol[i]->Pi0Mgg(),0.);
          double shE = 0;
          if(showercol[i]->SliceGeV()>1e-10)shE=showercol[i]->Energy()/showercol[i]->SliceGeV();
          double gap = showercol[i]->Gap();
          double eLLL = showercol[i]->DedxLongLL(int(jmshower::kElectron));
          double eLLT = showercol[i]->DedxTransLL(int(jmshower::kElectron));
          double muLLL = showercol[i]->DedxLongLL(int(jmshower::kMuon));
          double muLLT = showercol[i]->DedxTransLL(int(jmshower::kMuon));
          double eelLLL = showercol[i]->DedxLongLL(int(jmshower::kElectronSG));
          double eelLLT = showercol[i]->DedxTransLL(int(jmshower::kElectronSG));




          if(showerann > sh1ann){// Looking for the the shower with the largest electron ANN then set it's ANN as event ANN 
            ish1 = i;
            sh1ann = showerann;
            sh1gev = showerenergy;

          }

          if(showerenergy > sh1egev){// Looking for the leading shower then set it's ANN as event ANN  
            ish1e = i;
            sh1eann = showerann;
            sh1eanne = showeranne;
            sh1eannnocos = showerannnocos;
            sh1eannenocos = showerannenocos;
            sh1egev = showerenergy;
            sh1elength = showerlength;

            sh1egLLL = egLLL; 
            sh1egLLT = egLLT;
            sh1emuLLL = emuLLL;
            sh1emuLLT = emuLLT;
            sh1epi0LLL = epi0LLL;
            sh1epi0LLT = epi0LLT;
            sh1epLLL = epLLL;
            sh1epLLT = epLLT;
            sh1enLLL = enLLL;
            sh1enLLT = enLLT;
            sh1epiLLL = epiLLL;
            sh1epiLLT = epiLLT;
            sh1vtxgev = vtxgev;
            sh1pi0mass = pi0mass;
            sh1shE = shE;
            sh1gap = gap;
            sh1eLLL = eLLL;
            sh1eLLT = eLLT;
            sh1muLLL = muLLL;
            sh1muLLT = muLLT;

            sh1eelLLL = eelLLL;
            sh1eelLLT = eelLLT;


          }
          evtetot+=showerenergy;

//assign P4 based on shower PID, calculated by longitudinal+transverse likelihoods
          TVector3 p3(showercol[i]->Dir()[0]*showerenergy,showercol[i]->Dir()[1]*showerenergy,showercol[i]->Dir()[2]*showerenergy);
          double maxll = -9999;
          int lltype = -1;
          for( int itype = int(jmshower::kNull)+1; itype != int(jmshower::kLast); ++itype){
            if(itype!=4){
              if(showercol[i]->DedxLongLL(itype)+showercol[i]->DedxTransLL(itype)>maxll){
                lltype = itype;
                maxll = showercol[i]->DedxLongLL(itype)+showercol[i]->DedxTransLL(itype);
              }
            }
          }
          double m = xmass[0];
          if(lltype>=0&&lltype<8) m=xmass[lltype];

          double showerEtot0 = sqrt(p3.X()*p3.X()+p3.Y()*p3.Y()+p3.Z()*p3.Z());
          double P3 = sqrt(showerEtot0*showerEtot0+2*showerEtot0*m);
          double E = showerEtot0+m;
          TLorentzVector showerP40(showercol[i]->Dir()[0]*showerEtot0,showercol[i]->Dir()[1]*showerEtot0,showercol[i]->Dir()[2]*showerEtot0, showerEtot0);
          TLorentzVector showerP4(showercol[i]->Dir()[0]*P3, showercol[i]->Dir()[1]*P3, showercol[i]->Dir()[2]*P3, E);

//         TLorentzVector showerP4(showercol[i]->Dir()[0]*showerenergy, showercol[i]->Dir()[1]*showerenergy,showercol[i]->Dir()[2]*showerenergy,showerenergy);
          showerP40Col.push_back(showerP40); // 4-Momentum of showers without pid
          showerP4Col.push_back(showerP4); // 4-Momentum of showers
          showerpid.push_back(lltype);// PID of showers

        }

// looking for the secondary shower

        int ish2 = -1;
        double sh2gev = 0;

        for(unsigned int i = 0; i < showercol.size(); i++){
          if(i==(unsigned int)ish1e)continue;
          if(showercol[i]->Energy() > sh2gev){
            ish2 = i;
            sh2gev = showercol[i]->Energy();
          }
        }

// calculate total momentum, Pt and direction for the event 
        TLorentzVector evtSumP4(0,0,0,0);
        TLorentzVector evtSumP40(0,0,0,0);

        for(unsigned int ii=0;ii<showerP4Col.size();ii++){
          evtSumP4+=showerP4Col[ii];
        }

        for(unsigned int ii=0;ii<showerP40Col.size();ii++){
          evtSumP40+=showerP40Col[ii];
        }

        TVector3 nuDir = geom->NuMIBeamDirection(); //(-0.845e-3, 52.552e-3, 0.998618);//Beam dir
        double evtSumTheta = evtSumP4.Angle(nuDir);
        double evtSumPt = evtSumP4.Vect().Mag()*sin(evtSumTheta);
        double evtSumP = evtSumP4.Vect().Mag();
        double evtSumTheta0 = evtSumP40.Angle(nuDir);
        double evtSumPt0 = evtSumP40.Vect().Mag()*sin(evtSumTheta0);
        double evtSumP0 = evtSumP40.Vect().Mag();


        if(showercol.size()==1){
          evtetot+=showercol[0]->VtxGeV();
        }

        bool ismuon = showercol[ish1]->IsMuon(); // Muon like shower?
        bool ismuone = showercol[ish1e]->IsMuon(); // Muon like shower?

        if(sh1ann>0.84&&ismuon==true)Ncount1++;
//        if( ismuon == true ) return; // remove muon

        TVector3 beamDir = evtTrueNuP4.Vect();  // True Neutrino 4-momentum

// Osc  illation weigt for sigal
        double sigw = OscWeight(evtTrueNuPdg,evtTrueNuCCNC, evtTrueNuEnergy, true); // Signal 
        double bkgw = OscWeight(evtTrueNuPdg,evtTrueNuCCNC, evtTrueNuEnergy, false); // Background 

//   Neutrino energy reconstruction 
        double evtNuEnergy = showercol[ish1e]->NueEnergy(); 
//        double evtNuEnergy = NueSel::EnergyEstimator(beamDir, showerP4Col[ish1], evtetot, 1); 

 
        m_evtANN = sh1ann;    // Event ANN
        m_evtEtot = evtetot;  // Total detected energy
        m_evtSigOscW = sigw;  // Oscillation weight for swapping MC, obtain the expected signal energy distribution
        m_evtBkgOscW = bkgw;  // Oscillation weight for beam MC, obtain the expected background energy distribution
        m_evtNSh = evtNSh; // Number of shower in a event
        m_evtIsMuon = int(ismuon); // Is muon like?
        m_evtIsInvPhoton = int(showercol[ish1e]->IsInvPhoton()); // Is back-going photon?
   
        m_evtSh1DedxANN = sh1eann;    // Event ANN
        m_evtSh1DedxANNE = sh1eanne;    // Event ANN
        m_evtIsMuonSh1 = int(ismuone); // Is muon like?

        m_evtDetEnergy = evtetot;  // Total detected energy
        m_evtTrueNuEnergy = evtTrueNuEnergy; //True neutrnio energy
        m_evtNuEnergy = evtNuEnergy; //Recontructed neutrino energy
        m_evtSh1Energy = sh1egev; // Energy of the most energetic shower 
        m_evtSh1ExclEnergy = showercol[ish1e]->ExclEnergy(); // Energy of the most energetic shower 
        m_evtSh1Length = sh1elength; // Length of the most energetic shower
        m_evtSh1Start[0] = showercol[ish1e]->Start()[0]; 
        m_evtSh1Start[1] = showercol[ish1e]->Start()[1]; 
        m_evtSh1Start[2] = showercol[ish1e]->Start()[2]; 
        m_evtSh1Stop[0] = showercol[ish1e]->Stop()[0];
        m_evtSh1Stop[1] = showercol[ish1e]->Stop()[1];
        m_evtSh1Stop[2] = showercol[ish1e]->Stop()[2];
        m_evtSh1Dir[0] = showercol[ish1e]->Dir()[0];
        m_evtSh1Dir[1] = showercol[ish1e]->Dir()[1];
        m_evtSh1Dir[2] = showercol[ish1e]->Dir()[2];



        m_evtVtx[0] = evtTrueNuVtx[0]; // Event vertex
        m_evtVtx[1] = evtTrueNuVtx[1];
        m_evtVtx[2] = evtTrueNuVtx[2];



// Input variables for EID 
        m_evtSh1EGLLL=sh1egLLL;
        m_evtSh1EGLLT=sh1egLLT;
        m_evtSh1EMuLLL=sh1emuLLL;
        m_evtSh1EMuLLT=sh1emuLLT;
        m_evtSh1EPi0LLL=sh1epi0LLL;
        m_evtSh1EPi0LLT=sh1epi0LLT;
        m_evtSh1EPLLL=sh1epLLL;
        m_evtSh1EPLLT=sh1epLLT;
        m_evtSh1ENLLL=sh1enLLL;
        m_evtSh1ENLLT=sh1enLLT;
        m_evtSh1EPiLLL=sh1epiLLL;
        m_evtSh1EPiLLT=sh1epiLLT;
        m_evtSh1VtxGeV=sh1vtxgev;
        m_evtSh1Pi0Mass=sh1pi0mass;
        m_evtSh1ShE=sh1shE;
        m_evtSh1Gap=sh1gap;
        m_evtSh1ELLL=sh1eLLL;
        m_evtSh1ELLT=sh1eLLT;
        m_evtSh1MuLLL=sh1muLLL;
        m_evtSh1MuLLT=sh1muLLT;
        m_evtSh1GLLL=showercol[ish1e]->DedxLongLL(int(jmshower::kGamma));
        m_evtSh1InvGLLL=showercol[ish1e]->DedxInvLongLL(int(jmshower::kGamma));
        m_evtSh1DistToEdge = showercol[ish1e]->DistToEdge();
        m_evtSh1SliceEdge2GeV = showercol[ish1e]->SliceEdge2GeV();
        m_evtSh1SliceEdge5GeV = showercol[ish1e]->SliceEdge5GeV();
        m_evtSh1SliceEdge10GeV = showercol[ish1e]->SliceEdge10GeV();

        if(sh1ann>0.84&&sh1gev>0&&ismuon==false)Ncount2++;
//        std::cout<<"N NueCC "<<Ncount2<<std::endl;
        Nfill++;
        fEvent->Fill(); // Fill NTuple


        eid.SetANN(sh1eann);
        eid.SetANNE(sh1eanne);
        eid.SetANNEL(showercol[ish1e]->ANN(5));
        eid.SetANNEPi0(showercol[ish1e]->ANN(6));
        eid.SetANNEPi0EL(showercol[ish1e]->ANN(7));
//        std::cout<<"in nue sel ANNEL, ANNEPi0, ANNEPi0EL "<<eid.ANNEL()<<" "<<eid.ANNEPi0()<<" "<<eid.ANNEPi0EL()<<std::endl;
        eid.SetANNNoCos(sh1eannnocos);
        eid.SetANNENoCos(sh1eannenocos);



        eid.SetVtx(evtTrueNuVtx);
        eid.SetNuEnergy(evtNuEnergy);
        eid.SetDetEnergy(evtetot);

        eid.SetEvtNCell(slice.NCell());
//        eid.SetEvtSumCosTheta(cos(evtSumTheta));
        eid.SetEvtSumCosTheta(cos(showercol[ish1e]->Theta()));
        eid.SetEvtSumP(evtSumP);
        eid.SetEvtSumPt(evtSumPt);
        eid.SetEvtSumP0(evtSumP0);
        eid.SetEvtSumPt0(evtSumPt0);

        eid.SetEvtEtot(showercol[ish1e]->SliceEtot());
        eid.SetEvtMinX(evtMinX);
        eid.SetEvtMinY(evtMinY);
        eid.SetEvtMinZ(evtMinZ);
        eid.SetEvtMaxX(evtMaxX);
        eid.SetEvtMaxY(evtMaxY);
        eid.SetEvtMaxZ(evtMaxZ);
        eid.SetEvtGapTNS(evtGapTNS);
        eid.SetEvtNCellToEdge(evtNCellToEdge);

        eid.SetSh1Energy(sh1egev);
        eid.SetSh1ExclEnergy(showercol[ish1e]->ExclEnergy());
        eid.SetSh1TotalLength(sh1elength);
        eid.SetSh1Start(showercol[ish1e]->Start());
        eid.SetSh1Stop(showercol[ish1e]->Stop());
        eid.SetSh1Dir(showercol[ish1e]->Dir());

        eid.SetSh1VtxDoca(showercol[ish1e]->VtxDoca());
        eid.SetSh1Gap(showercol[ish1e]->Gap());
        eid.SetSh1NPlane(showercol[ish1e]->NPlane());
        eid.SetSh1XNPlane(showercol[ish1e]->XNPlane());
        eid.SetSh1YNPlane(showercol[ish1e]->YNPlane());
        eid.SetSh1NCell(showercol[ish1e]->NCell());
        eid.SetSh1XNCell(showercol[ish1e]->NXCell());
        eid.SetSh1YNCell(showercol[ish1e]->NYCell());
        eid.SetSh1PID(showerpid[ish1e]);
        eid.SetSh1NCellToEdge(showercol[ish1e]->NCellToEdge());

        if(ish2!=-1){
          double evtSh1Sh2Dang = 180*showerP4Col[ish1e].Angle(showerP4Col[ish2].Vect())/3.14159265;
          eid.SetSh1Sh2Dang(evtSh1Sh2Dang);
          eid.SetSh2Energy(showercol[ish2]->Energy());
          eid.SetSh2ExclEnergy(showercol[ish2]->ExclEnergy());
          eid.SetSh2TotalLength(showercol[ish2]->TotalLength());
          eid.SetSh2Start(showercol[ish2]->Start());
          eid.SetSh2Stop(showercol[ish2]->Stop());
          eid.SetSh2Dir(showercol[ish2]->Dir());

          eid.SetSh2VtxDoca(showercol[ish2]->VtxDoca());
          eid.SetSh2Gap(showercol[ish2]->Gap());
          eid.SetSh2NPlane(showercol[ish2]->NPlane());
          eid.SetSh2XNPlane(showercol[ish2]->XNPlane());
          eid.SetSh2YNPlane(showercol[ish2]->YNPlane());
          eid.SetSh2NCell(showercol[ish2]->NCell());
          eid.SetSh2XNCell(showercol[ish2]->NXCell());
          eid.SetSh2YNCell(showercol[ish2]->NYCell());
          eid.SetSh2PID(showerpid[ish2]);
          eid.SetSh2NCellToEdge(showercol[ish2]->NCellToEdge());

        }else{
          eid.SetSh1Sh2Dang(-5);
          eid.SetSh2Energy(0);
          eid.SetSh2ExclEnergy(0);
          eid.SetSh2TotalLength(0);
          TVector3 v3null(0,0,0);
          eid.SetSh2Start(v3null);
          eid.SetSh2Stop(v3null);
          eid.SetSh2Dir(v3null);
          eid.SetSh2VtxDoca(0);
          eid.SetSh2Gap(0);
          eid.SetSh2NPlane(0);
          eid.SetSh2XNPlane(0);
          eid.SetSh2YNPlane(0);
          eid.SetSh2NCell(0);
          eid.SetSh2XNCell(0);
          eid.SetSh2YNCell(0);
          eid.SetSh2PID(-1);
          eid.SetSh2NCellToEdge(10000);
        }


        eid.SetIsMuon(ismuone);
        eid.SetIsInvPhoton(showercol[ish1e]->IsInvPhoton());
        eid.SetSigOscW(sigw);
        eid.SetBkgOscW(bkgw);

        eid.SetDedx0(showercol[ish1e]->PlaneDedx(0));
        eid.SetDedx1(showercol[ish1e]->PlaneDedx(1));
        eid.SetDedx2(showercol[ish1e]->PlaneDedx(2));
        eid.SetDedx3(showercol[ish1e]->PlaneDedx(3));
        eid.SetDedx4(showercol[ish1e]->PlaneDedx(4));
        eid.SetDedx5(showercol[ish1e]->PlaneDedx(5));

        eid.SetEGLLL(sh1egLLL);
        eid.SetEGLLT(sh1egLLT);
        eid.SetEMuLLL(sh1emuLLL);
        eid.SetEMuLLT(sh1emuLLT);
        eid.SetEPi0LLL(sh1epi0LLL);
        eid.SetEPi0LLT(sh1epi0LLT);
        eid.SetEPLLL(sh1epLLL);
        eid.SetEPLLT(sh1epLLT);
        eid.SetENLLL(sh1enLLL);
        eid.SetENLLT(sh1enLLT);
        eid.SetEPiLLL(sh1epiLLL);
        eid.SetEPiLLT(sh1epiLLT);
        eid.SetVtxGeV(sh1vtxgev);
        eid.SetPi0Mass(sh1pi0mass);
        eid.SetShE(sh1shE);
        eid.SetGap(sh1gap);
        eid.SetELLL(sh1eLLL);
        eid.SetELLT(sh1eLLT);
        eid.SetEelLLL(sh1eelLLL);
        eid.SetEelLLT(sh1eelLLT);
        eid.SetMuLLL(sh1muLLL);
        eid.SetMuLLT(sh1muLLT);
        eid.SetGLLL(showercol[ish1e]->DedxLongLL(int(jmshower::kGamma)));
        eid.SetInvGLLL(showercol[ish1e]->DedxInvLongLL(int(jmshower::kGamma)));


        eid.SetSh1DistToEdge(showercol[ish1e]->DistToEdge()); // Shower distance to edge
        eid.SetEEdge2Cell(showercol[ish1e]->SliceEdge2GeV());  //Energy deposit within 2 cells to detector edges
        eid.SetEEdge5Cell(showercol[ish1e]->SliceEdge5GeV());  //Energy deposit within 5 cells to detector edges
        eid.SetEEdge10Cell(showercol[ish1e]->SliceEdge10GeV()); //Energy deposit within 10 cells to detector edges
        eidcol->push_back(eid);
      }

      util::CreateAssn(*this,evt,*eidcol,art::Ptr<rb::Cluster>(slicecol,sliceIdx),*assns);
    }


    //unsigned int iSlice = showercol[ish1e]->ISlice();
    //std::cout << "iSlice: " << iSlice << std::endl;
    //util::CreateAssn(*this,evt,*eidcol,slicelist[iSlice],*assns);
    //assert(iSlice >=0 && iSlice < slicecol->size());


    evt.put(std::move(eidcol));
    evt.put(std::move(assns));
    Nput++;
//    if(Ncount0%100==0)std::cout<<Ncount0<<" events processed, "<<Ncount2<<" Pass NueCC Sel  "<<Ncount1<<" muon like"<<std::endl;
    return;
  }

} //namespace jmshower


////////////////////////////////////////////////////////////////////////
namespace jmshower
{
  DEFINE_ART_MODULE(NueSel)
}