Xeff_module.cc

Go to the documentation of this file.

///
/// \brief 
/// \author Xuebing Bu - xbbu@fnal.gov
///
#include "ShowerLID/EventLID.h"
#include "ShowerLID/ShowerLID.h"
#include "Calibrator/Calibrator.h"

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

// ROOT includes
#include "TVector3.h"
#include "TH1F.h"
#include "TMath.h"
#include "TTree.h"
#include "TProfile2D.h"
#include "time.h"

#include "nusimdata/SimulationBase/MCFlux.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "Simulation/Particle.h"
#include "Simulation/ParticleNavigator.h"
#include "SummaryData/POTSum.h"
#include "SummaryData/SpillData.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/RecoHit.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/Prong.h"
#include "RecoBase/HoughResult.h"

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

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "canvas/Persistency/Common/PtrVector.h"
#include "canvas/Persistency/Common/Assns.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindMany.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include <cmath>
#include <vector>
#include <string>
#include <vector>

//include TMVA
#include "TMVA/Tools.h"
#include "TMVA/Reader.h"
#include "TMVA/MethodCuts.h"

class TVector3;
class TH1F;
class TTree;

namespace rb{class Cluster; class Track;}
namespace sim{class Particle;}
namespace simb{class MCFlux; class MCTruth; class MCNeutrino;}

namespace ncs {
  class Xeff : public art::EDAnalyzer {
  public:
    explicit Xeff(fhicl::ParameterSet const& pset);
    virtual ~Xeff();
    
    void beginJob();
    void analyze(const art::Event& evt);
    void reconfigure(const fhicl::ParameterSet& p);
    void endSubRun(const art::SubRun& sr);
    void endJob(); 
    double fTotalPOT;
    int fTotalSpill;
    TTree *_btree;

  private:

    const rb::RecoHit MakeRecoHit(art::Ptr<rb::CellHit> const& chit,
                                  double const& w);
    
    //folders to get objects from
    std::string fSliceLabel;     ///label for slices
    std::string fKTrackLabel;     ///label for kalman tracks
    std::string fProngLabel;     ///label for prong 
    std::string fInstLabel3D;    ///instance label for prongs 3D
    std::string fInstLabel2D;    ///instance label for prongs 2D
    std::string fMCFluxLabel;    ///label for generator information
    std::string fVertexLabel;    ///label for vertex
    std::string fCellHitLabel;   ///label for cell hits
    std::string flidLabel;

    //save event info to ntuples
    int runnum;
    int subrunnum;
    int evtnum;
    int Nslices;
    double spillPOT;
    int Ndeaddcm;
  
    TTree *_otree;
    //TTree *_evttree;
    int sl_runnum;
    int sl_subrunnum;
    int sl_evtnum;
    int sl_Nslices;
    double sl_spillPOT;
    int sl_Ndeaddcm;

    /*** true info ***/
    float truePt;//true Pt
    float trueE;//true E
    float trueTheta;//true theta
    int ccnc;
    int PDG;
    int origPDG;
    int mode;
    bool isQE;
    int intType;
    int hitNucl;
    float HadX;
    float HadY;
    float HadW;
    float qsqr;
    float OscWt; 
    int lepPDG;
    float lepPx;
    float lepPy;
    float lepPz;
    float lepE;
    int mesonPDG;
    float mesonPx;
    float mesonPy;
    float mesonPz;
    float Vx;//true vertex X
    float Vy;//true vertex Y
    float Vz;//true vertex Z

    /*** Reco Vertex ***/
    int hasVertex;
    int hasOneVertex;
    float recoVx;//reco vertex X
    float recoVy;//reco vertex Y
    float recoVz;//reco vertex Z

    /*** slicer information ***/
    int Ncells;//# of associate cells
    int Nplanes;//# of extended planes
    float Eslice;//reco energy
    int sliceIndex;
    
    /*** longest kalman track ***/
    float KtrackLength;//length of longest track
    int KtrackNcells;//# of cells associated to longest track
    float KtrackEnergy;//energy associated to longest track
    float KtrackNcells_ratio;//# of cells associated to longest track over Ncells
    float KtrackEnergy_ratio;//energy associated to longest track over recoE
    float KtrackNplanes;
    float KtrackTheta;
    float KtrackPhi;
    float KtrackStartX;
    float KtrackStartY;
    float KtrackStartZ;
    float KtrackStopX;
    float KtrackStopY;
    float KtrackStopZ;
    int NKtracks;//# of tracks
    int Ktrack3D;

    /*** mip information ***/
    float Nmip;//# of mip cells
    float Fmip;//# of mip cells over Ncells
    
    /*** prong information ***/
    int Nprongs;
    float Ebalance;
    float Dphi;

    float p1Energy;
    float p1Phi;
    float p1Theta;
    float p1StartX;
    float p1StartY;
    float p1StartZ;
    float p1Length;
    float p1StopX;
    float p1StopY;
    float p1StopZ;
    float p1EnergyXView;//energy for XView cells
    float p1EnergyYView;//energy for YView cells
    int p1NCells;//# of cells
    int p1NCellsXView;//# of cells for XView
    int p1NCellsYView;//# of cells for YView
    int p1Nplanes;
    float p1Fmip;
    //bdt inputs
    float p2Fmip;
    float emaxfrac_6p;
    float efrac_10p;
    float efrac_p2;
    float efrac_p3;
    float efrac_p4;
    float efrac_2sig;

    //leading prong
    int prong1cellNcell;
    float prong1cellX[500];
    float prong1cellY[500];
    float prong1cellZ[500];
    float prong1cellE[500];
    int prong1cellADC[500];
    int prong1cellPlane[500];
    int prong1cellCell[500];
    float prong1cellPECorr[500];
    float prong1cellW[500];

    float lid;
  };
}

namespace ncs{
  
/*
  //BDT reader
  TMVA::Reader *_bdt=new TMVA::Reader();
  //input variable values
  float inputvars[15]={0.};
  //input variable names
  TString inputs[15]={"_ncells","_energy","_ltrack","_ltrack_ncells_ratio","_fmip",
                      "_nmip","_efrac_20planes","_efrac_slide2","_efrac_slide6",
                      "_efrac_2sigma","_eiso_3sigma","_nprongs","_ebalance",
                      "_nprongs2D","_ebalance2D"};
*/
  //......................................................................
  Xeff::Xeff(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
  {
    mf::LogInfo("Xeff")<<__PRETTY_FUNCTION__<<"\n";
    reconfigure(pset);
  }
  
  //......................................................................
  Xeff::~Xeff()
  {
  }
  
  //......................................................................
  void Xeff::reconfigure(const fhicl::ParameterSet& pset)
  {
    fSliceLabel   = pset.get< std::string >("SliceLabel");
    fKTrackLabel   = pset.get< std::string >("KTrackLabel");
    fProngLabel   = pset.get< std::string >("ProngLabel");
    fInstLabel3D  = pset.get< std::string >("InstLabel3D");
    fInstLabel2D  = pset.get< std::string >("InstLabel2D");
    fMCFluxLabel  = pset.get< std::string >("MCFluxLabel");
    fCellHitLabel = pset.get< std::string >("CellHitLabel");
    fVertexLabel  = pset.get< std::string >("VertexLabel");

    flidLabel     = pset.get< std::string >("lidLabel");
  }

  //......................................................................
  void Xeff::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;    

    fTotalPOT=0.;
    fTotalSpill=0;

/*
    for( int i=0; i<15; ++i ){
      _bdt->AddVariable(inputs[i],&inputvars[i]);
    }
    
    _bdt->BookMVA("BDT method","nonoscBDTweights_20140304.xml");
*/  
  
    //histo
    //Ereco=tfs->make<TH1F>("Ereco","Ereco;Ereco;# of events",100,0.,5.);
    //personal tuples
    _otree = tfs->make<TTree>("Xeff","Xeff particle");
    
    _otree->Branch("sl_runnum",&sl_runnum,"sl_runnum/I");
    _otree->Branch("sl_subrunnum",&sl_subrunnum,"sl_subrunnum/I");
    _otree->Branch("sl_evtnum",&sl_evtnum,"sl_evtnum/I");
    _otree->Branch("sl_Nslices",&sl_Nslices,"sl_Nslices/I");
    _otree->Branch("sl_spillPOT",&sl_spillPOT,"sl_spillPOT/D");
    _otree->Branch("sl_Ndeaddcm",&sl_Ndeaddcm,"sl_Ndeaddcm/I");

    _otree->Branch("truePt",&truePt,"truePt/F");
    _otree->Branch("trueE",&trueE,"trueE/F");
    _otree->Branch("trueTheta",&trueTheta,"trueTheta/F");
    _otree->Branch("ccnc",&ccnc,"ccnc/I");
    _otree->Branch("PDG",&PDG,"PDG/I"); 
    _otree->Branch("origPDG",&origPDG,"origPDG/I"); 
    _otree->Branch("mode",&mode,"mode/I");
    _otree->Branch("isQE",&isQE,"isQE/B");
    _otree->Branch("intType",&intType,"intType/I");
    _otree->Branch("hitNucl",&hitNucl,"hitNucl/I");
    _otree->Branch("HadX",&HadX,"HadX/F");
    _otree->Branch("HadY",&HadY,"HadY/F");
    _otree->Branch("HadW",&HadW,"HadW/F");
    _otree->Branch("qsqr",&qsqr,"qsqr/F");
    _otree->Branch("OscWt",&OscWt,"OscWt/F");
    _otree->Branch("lepE",&lepE,"lepE/F");
    _otree->Branch("lepPx",&lepPx,"lepPx/F");
    _otree->Branch("lepPy",&lepPy,"lepPy/F");
    _otree->Branch("lepPz",&lepPz,"lepPz/F");
    _otree->Branch("lepPDG",&lepPDG,"lepPDG/I");
    _otree->Branch("mesonPDG",&mesonPDG,"mesonPDG/I");
    _otree->Branch("mesonPx",&mesonPx,"mesonPx/F");
    _otree->Branch("mesonPy",&mesonPy,"mesonPy/F");
    _otree->Branch("mesonPz",&mesonPz,"mesonPz/F");
    _otree->Branch("Vx",&Vx,"Vx/F");
    _otree->Branch("Vy",&Vy,"Vy/F");
    _otree->Branch("Vz",&Vz,"Vz/F");
    _otree->Branch("hasVertex",&hasVertex,"hasVertex/I");
    _otree->Branch("hasOneVertex",&hasOneVertex,"hasOneVertex/I");
    _otree->Branch("recoVx",&recoVx,"recoVx/F");
    _otree->Branch("recoVy",&recoVy,"recoVy/F");
    _otree->Branch("recoVz",&recoVz,"recoVz/F");
   
    _otree->Branch("NKtracks",&NKtracks,"NKtracks/I");
    _otree->Branch("KtrackLength",&KtrackLength,"KtrackLength/F");
    _otree->Branch("KtrackNcells",&KtrackNcells,"KtrackNcells/I");
    _otree->Branch("KtrackEnergy",&KtrackEnergy,"KtrackEnergy/F");
    _otree->Branch("KtrackNcells_ratio",&KtrackNcells_ratio,"KtrackNcells_ratio/F");
    _otree->Branch("KtrackEnergy_ratio",&KtrackEnergy_ratio,"KtrackEnergy_ratio/F");
    _otree->Branch("KtrackNplanes",&KtrackNplanes,"KtrackNplanes/F");
    _otree->Branch("KtrackPhi",&KtrackPhi,"KtrackPhi/F");
    _otree->Branch("KtrackTheta",&KtrackTheta,"KtrackTheta/F");
    _otree->Branch("KtrackStartX",&KtrackStartX,"KtrackStartX/F");
    _otree->Branch("KtrackStartY",&KtrackStartY,"KtrackStartY/F");
    _otree->Branch("KtrackStartZ",&KtrackStartZ,"KtrackStartZ/F");
    _otree->Branch("KtrackStopX",&KtrackStopX,"KtrackStopX/F");
    _otree->Branch("KtrackStopY",&KtrackStopY,"KtrackStopY/F");
    _otree->Branch("KtrackStopZ",&KtrackStopZ,"KtrackStopZ/F");
    _otree->Branch("Ktrack3D",&Ktrack3D,"Ktrack3D/I");

    _otree->Branch("Ncells",&Ncells,"Ncells/I");
    _otree->Branch("Nplanes",&Nplanes,"Nplanes/I");
    _otree->Branch("Eslice",&Eslice,"Eslice/F");
    _otree->Branch("sliceIndex",&sliceIndex,"sliceIndex/I");
    
    _otree->Branch("Nmip",&Nmip,"Nmip/F");
    _otree->Branch("Fmip",&Fmip,"Fmip/F");

    _otree->Branch("Nprongs",&Nprongs,"Nprongs/I");
    _otree->Branch("Ebalance",&Ebalance,"Ebalance/F");
    _otree->Branch("Dphi",&Dphi,"Dphi/F");
    
    _otree->Branch("p1Length",&p1Length,"p1Length/F");
    _otree->Branch("p1StopX",&p1StopX,"p1StopX/F");
    _otree->Branch("p1StopY",&p1StopY,"p1StopY/F");
    _otree->Branch("p1StopZ",&p1StopZ,"p1StopZ/F");
    _otree->Branch("p1Energy",&p1Energy,"p1Energy/F");
    _otree->Branch("p1Phi",&p1Phi,"p1Phi/F");
    _otree->Branch("p1Theta",&p1Theta,"p1Theta/F");
    _otree->Branch("p1StartX",&p1StartX,"p1StartX/F");
    _otree->Branch("p1StartY",&p1StartY,"p1StartY/F");
    _otree->Branch("p1StartZ",&p1StartZ,"p1StartZ/F");
    _otree->Branch("p1EnergyXView",&p1EnergyXView,"p1EnergyXView/F");
    _otree->Branch("p1EnergyYView",&p1EnergyYView,"p1EnergyYView/F");
    _otree->Branch("p1NCells",&p1NCells,"p1NCells/I");
    _otree->Branch("p1NCellsXView",&p1NCellsXView,"p1NCellsXView/I");
    _otree->Branch("p1NCellsYView",&p1NCellsYView,"p1NCellsYView/I");
    _otree->Branch("p1Nplanes",&p1Nplanes,"p1Nplanes/I");
    _otree->Branch("p1Fmip",&p1Fmip,"p1Fmip/F");

    _otree->Branch("p2Fmip",&p2Fmip,"p2Fmip/F");
    _otree->Branch("emaxfrac_6p",&emaxfrac_6p,"emaxfrac_6p/F");
    _otree->Branch("efrac_10p",&efrac_10p,"efrac_10p/F");
    _otree->Branch("efrac_p2",&efrac_p2,"efrac_p2/F");
    _otree->Branch("efrac_p3",&efrac_p3,"efrac_p3/F");
    _otree->Branch("efrac_p4",&efrac_p4,"efrac_p4/F");
    _otree->Branch("efrac_2sig",&efrac_2sig,"efrac_2sig/F");

    _otree->Branch("prong1cellNcell",&prong1cellNcell,"prong1cellNcell/I");
    _otree->Branch("prong1cellX",&prong1cellX,"prong1cellX[prong1cellNcell]/F");
    _otree->Branch("prong1cellY",&prong1cellY,"prong1cellY[prong1cellNcell]/F");
    _otree->Branch("prong1cellZ",&prong1cellZ,"prong1cellZ[prong1cellNcell]/F");
    _otree->Branch("prong1cellE",&prong1cellE,"prong1cellE[prong1cellNcell]/F");
    _otree->Branch("prong1cellADC",&prong1cellADC,"prong1cellADC[prong1cellNcell]/I");
    _otree->Branch("prong1cellPlane",&prong1cellPlane,"prong1cellPlane[prong1cellNcell]/I");
    _otree->Branch("prong1cellCell",&prong1cellCell,"prong1cellCell[prong1cellNcell]/I");
    _otree->Branch("prong1cellW",&prong1cellW,"prong1cellW[prong1cellNcell]/F");
    _otree->Branch("prong1cellPECorr",&prong1cellPECorr,"prong1cellPECorr[prong1cellNcell]/F");
    
    _otree->Branch("lid",&lid,"lid/F");

    _btree=new TTree("beamInfo","beam info");
    _btree->Branch("runnum",&runnum,"runnum/I");
    _btree->Branch("subrunnum",&subrunnum,"subrunnum/I");
    _btree->Branch("evtnum",&evtnum,"evtnum/I");
    _btree->Branch("Nslices",&Nslices,"Nslices/I");
    _btree->Branch("spillPOT",&spillPOT,"spillPOT/D");
    _btree->Branch("Ndeaddcm",&Ndeaddcm,"Ndeaddcm/I");

  }

  //......................................................................
  void Xeff::endSubRun(const art::SubRun& sr)
  {
    art::Handle<sumdata::POTSum> pot;
    sr.getByLabel(fMCFluxLabel, pot);
    
    
    if( !pot.failedToGet() ){
      fTotalPOT += pot->totgoodpot;
      fTotalSpill += pot->goodspills;
      
      std::cout << "POT in subrun " << sr.subRun()
                << " = " << pot->totgoodpot <<" with "<<fTotalSpill<<" spills"<< std::endl;
    }
  }

  //......................................................................
  const rb::RecoHit Xeff::MakeRecoHit(art::Ptr<rb::CellHit> const& chit,
                                       double const& w)
  {
    art::ServiceHandle<calib::Calibrator> cal;
    const rb::RecoHit rhit(cal->MakeRecoHit(*chit, w));
    
    return rhit;
  }
  
  //......................................................................
  void Xeff::analyze(const art::Event& evt) 
  {
    sl_runnum=-1;
    sl_subrunnum=-1;
    sl_evtnum=-1;
    sl_Nslices=-1;
    sl_spillPOT=-1.;
    sl_Ndeaddcm=-1;
    Ndeaddcm=-1;

    runnum=-1;
    subrunnum=-1;
    evtnum=-1;
    Nslices=-1;
    spillPOT=0.;
    
    lid=-99.;

    Vx=-999.;
    Vy=-999.;
    Vz=-999.;
    hasVertex=-1;
    hasOneVertex=-1;
    recoVx=-999.;
    recoVy=-999.;
    recoVz=-999.;

    truePt=-1.;
    trueE=-1.;
    ccnc=-1;
    PDG=-1;
    origPDG=-1;
    mode=-1;
    isQE=false;
    intType=-1;
    hitNucl=-1;
    HadX=-999.;
    HadY=-999.;
    HadW=-999.;
    qsqr=-999.;
    OscWt=-999.;

    lepPDG=-1;
    lepPx=-1.;
    lepPy=-1.; 
    lepPz=-1.;
    mesonPDG=-1;
    mesonPx=-1.;
    mesonPy=-1.;
    mesonPz=-1.;
    
    Ncells=0;
    Nplanes=0;
    Eslice=-1.;
    sliceIndex=-1;
       
    KtrackLength=0.;
    KtrackNcells=0;
    KtrackEnergy=0.;
    KtrackEnergy_ratio=0.;
    KtrackNcells_ratio=0.;
    KtrackNplanes=0;
    KtrackTheta=-1.;
    KtrackPhi=-1.;
    KtrackStartX=-999.;
    KtrackStartY=-999.;
    KtrackStartZ=-999.;
    KtrackStopX=-999.;
    KtrackStopY=-999.;
    KtrackStopZ=-999.;
    NKtracks=0.;

    Nmip=-1.;
    Fmip=-1.;
  
    Nprongs=0;
    Ebalance=0.;
    Dphi=-1.;

    p1Energy=0.;
    p1Phi=-1.;
    p1Theta=-1.;
    p1StartX=-999.;
    p1StartY=-999.;
    p1StartZ=-999.;
    p1StopX=-999.;
    p1StopY=-999.;
    p1StopZ=-999.;
    p1Length=-1.;
    p1EnergyXView=0.;
    p1EnergyYView=0.;
    p1NCells=0;
    p1NCellsXView=0;
    p1NCellsYView=0;
    p1Nplanes=0;
    p1Fmip=-1.;
    p2Fmip=0.;
    emaxfrac_6p=0.;
    efrac_10p=0.;
    efrac_p2=0.;
    efrac_p3=0.;
    efrac_p4=0.;
    efrac_2sig=0.;    

    prong1cellNcell=0;
    for( int i=0; i<500; ++i ){
      prong1cellX[i]=-999.;
      prong1cellY[i]=-999.;
      prong1cellZ[i]=-999.;
      prong1cellE[i]=-999.;
      prong1cellADC[i]=0;
      prong1cellPlane[i]=-1;
      prong1cellCell[i]=-1;
      prong1cellW[i]=-1.;
      prong1cellPECorr[i]=-1.;
      
    }
    
    art::ServiceHandle<calib::Calibrator> cal;
    art::ServiceHandle<geo::Geometry> geom;
    
    //event info
    runnum = evt.run();
    subrunnum = evt.subRun();
    evtnum = evt.id().event();

    //get all hits info
    art::Handle<std::vector<rb::CellHit> > chits;
    evt.getByLabel(fCellHitLabel, chits);

    std::vector<int> sq_plane;
    sq_plane.clear();
    std::vector<int> sq_xview;
    sq_xview.clear();
    std::vector<float> sq_tns;
    sq_tns.clear();
    std::vector<int> sq_cell;
    sq_cell.clear();

    int sq_ncells=chits->size();
    for(int i = 0; i < sq_ncells; ++i){
      art::Ptr< rb::CellHit> hit(chits, i);
      double hittime=hit->TNS()/1000.;
       sq_plane.push_back(hit->Plane());
       sq_tns.push_back(hittime);
       sq_cell.push_back(hit->Cell());
       if( (hit->View())==geo::kX ) sq_xview.push_back(1);
       else sq_xview.push_back(0);
    }//loop all hits

    int Nmissingdcm[14]={0};
    for( int ic=0; ic<sq_ncells; ++ic ){
      if( sq_plane[ic]<64 ){
        if( sq_xview[ic]==1 ){
          if( sq_cell[ic]>63 ) Nmissingdcm[0]++;
          else Nmissingdcm[1]++;
        }//XZ view
        else {
          if( sq_cell[ic]>31 ) Nmissingdcm[2]++;
          else  Nmissingdcm[3]++;
        }//YZ view
      }//DB1
      else if( sq_plane[ic]<128 ){
        if( sq_xview[ic]==1 ){
          if( sq_cell[ic]>63 ) Nmissingdcm[4]++;
          else Nmissingdcm[5]++;
        }//XZ view
        else {
          if( sq_cell[ic]>31 ) Nmissingdcm[6]++;
          else Nmissingdcm[7]++;
        }//YZ view
      }//DB2
      else if( sq_plane[ic]<192 ){
        if( sq_xview[ic]==1 ){
          if( sq_cell[ic]>63 ) Nmissingdcm[8]++;
          else Nmissingdcm[9]++;
        }//XZ view
        else{
          if( sq_cell[ic]>31 ) Nmissingdcm[10]++; 
          else Nmissingdcm[11]++;
        }//YZ view
      }//DB3  
      else{
        if( sq_xview[ic]==1 ) Nmissingdcm[12]++;
        else Nmissingdcm[13]++;
      }//muon-catcher
    }//loop all cells

    int ndeaddcms=0;
    for( int id=0; id<14; ++id ){
      if( Nmissingdcm[id]==0 ) ndeaddcms++;
    }

    Ndeaddcm=ndeaddcms; 

    ////////////////////////
    //get reco slicer info// 
    ////////////////////////
    art::Handle< std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSliceLabel, slices);
    
    //get tracks, verticies and showers associated with the slices
    //art::FindManyP<rb::Track> fmDTrack(slices, evt, fDTrackLabel);
    art::FindManyP<rb::Track> fmKTrack(slices, evt, fKTrackLabel);
    art::FindManyP<rb::Vertex> fmVertex(slices, evt, fVertexLabel);
    //art::FindManyP<jmshower::JMShower> fmShower(slices, evt, fShowerLabel);
    
    //std::cout<<evt.run()<<" "<<evt.subRun()<<" "<<evt.id().event()<<std::endl;
    //std::cout<<"# of slices: "<<slices->size()<<std::endl;
    
    int nslice=0;
    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
      
      const rb::Cluster& slice = (*slices)[sliceIdx];
      if(slice.IsNoise()) continue;
      ++nslice;
      
    }//loop slices
    
    Nslices=nslice;
    //std::cout<<"# of slice: "<<nslice<<std::endl;
    
    _btree->Fill();
    
    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
      
      const rb::Cluster& slice = (*slices)[sliceIdx];
      if(slice.IsNoise()) continue;
      
      //get MC truth info
      art::ServiceHandle<cheat::BackTracker> bt;
      art::PtrVector<rb::CellHit> allhits;
      for(size_t h = 0; h < chits->size(); ++h)
        allhits.push_back(art::Ptr<rb::CellHit>(chits, h));
      const art::Ptr<rb::Cluster> allslice(slices,sliceIdx);
      std::vector<cheat::NeutrinoEffPur> iacts = bt->SliceToNeutrinoInteractions(allslice, allhits);
      
      if( iacts.empty() ){
        std::cout<<"no matching between slice and nu interactions!"<<std::endl;
        return;
      }
      // Pick the best-matching one
      art::Ptr<simb::MCTruth>& truth = iacts[0].neutrinoInt;
      // Assn interface is needlessly confusing. Have to pack stuff up like this
      art::PtrVector<simb::MCTruth> pv;
      pv.push_back(truth);
      // Get the flux associated with this truth
      //art::PtrVector<simb::MCFlux> fluxs = util::FindManyP<simb::MCFlux>(pv, evt, "geantgen", 0);
      art::FindManyP<simb::MCFlux> fmf(pv, evt, fMCFluxLabel);
      std::vector< art::Ptr<simb::MCFlux> > fluxs = fmf.at(0);
      simb::MCFlux flux;
      if(fluxs.size() == 1){
        flux = *fluxs[0];
      }
      else{
        std::cout<<"There is no assoicated neutrino intereaction !"<<std::endl;
      }
      const simb::MCNeutrino& nu = truth->GetNeutrino();
      
      //select nue CC only events
      if( !(abs(flux.fntype)==12 && abs(nu.Nu().PdgCode())==12 ) ) continue;
      if( !(nu.CCNC()==0) ) continue;
      
      trueE   = nu.Nu().E();
      truePt  = nu.Pt();
      trueTheta = nu.Theta();
      ccnc    = nu.CCNC();
      PDG     = nu.Nu().PdgCode();
      origPDG = flux.fntype;
      mode    = nu.Mode();
      isQE    = nu.InteractionType() == simb::kQE ||
        nu.InteractionType() == simb::kCCQE;
      intType = nu.InteractionType();
      hitNucl = nu.HitNuc();
      HadX    = nu.X();
      HadY    = nu.Y();
      HadW    = nu.W();
      qsqr    = nu.QSqr();
      
      Vx=nu.Nu().Vx();
      Vy=nu.Nu().Vy();
      Vz=nu.Nu().Vz();
      
      if( fabs(nu.Nu().Vx())>200. ) continue;
      if( fabs(nu.Nu().Vy())>200. ) continue;
      if( nu.Nu().Vz()<0. ) continue;

      lepPDG = nu.Lepton().PdgCode();
      lepE = nu.Lepton().E();
      lepPx = nu.Lepton().Px();
      lepPy = nu.Lepton().Py();
      lepPz = nu.Lepton().Pz();
      
      art::Handle<std::vector<simb::MCTruth> > mclist;
      evt.getByLabel(fMCFluxLabel,mclist);
      art::Handle< std::vector<simb::MCFlux> > fllist;
      evt.getByLabel(fMCFluxLabel, fllist);
      if( mclist->size()>0 && fllist->size()>0 ){
        for( unsigned int i_intx=0; i_intx<mclist->size(); ++i_intx){
          art::Ptr<simb::MCTruth> mctruth(mclist,i_intx);
          const simb::MCNeutrino& mcnu = mctruth->GetNeutrino();
          TLorentzVector Tslnu(nu.Nu().Px(),nu.Nu().Py(),
                               nu.Nu().Pz(),nu.Nu().E());
          TLorentzVector Tmcnu(mcnu.Nu().Px(),mcnu.Nu().Py(),
                               mcnu.Nu().Pz(),mcnu.Nu().E());
          if(mcnu.Nu().PdgCode()==nu.Nu().PdgCode() &&
             mcnu.CCNC()==nu.CCNC() && mcnu.Mode()==nu.Mode() &&
             mcnu.Lepton().PdgCode()==nu.Lepton().PdgCode() &&
             Tslnu==Tmcnu ){
            art::Ptr<simb::MCFlux> flux(fllist,i_intx);
            mesonPDG = flux->ftptype;
            mesonPx = flux->ftpx;
            mesonPy = flux->ftpy;
            mesonPz = flux->ftpz;
          }
        }
      }
      
      //missing DCMs
      sl_Ndeaddcm=ndeaddcms;

      ///////////////// 
      //select vertex//
      ///////////////// 
      if (!(fmVertex.isValid())){
        hasVertex=0;
      }//no reco vertex
      else{
        hasVertex=1;
        std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
        if( vert.size() != 1 ){
          hasOneVertex=0;
        }// ! 1-vertex
        else{
          hasOneVertex=1;
          ///////////////////////////////////
          //select the longest kalman track//
          ///////////////////////////////////
          std::vector<art::Ptr<rb::Track>> Ktracks = fmKTrack.at(sliceIdx);
          double longestKTrack = 0;
          int longestKTrack_ncells=0;
          double longestKTrack_energy=0.;
          int longestKTrack_nplanes=0;
          float longestKTrack_phi=0.;
          float longestKTrack_theta=0.; 
          float longestKTrack_startx=0.;
          float longestKTrack_starty=0.;
          float longestKTrack_startz=0.;
          float longestKTrack_stopx=0.;
          float longestKTrack_stopy=0.;
          float longestKTrack_stopz=0.;
          int longestKTrack_is3D=0;
          
          if( fmKTrack.isValid() ){
            for(unsigned int n = 0; n < Ktracks.size(); ++n){
              //if( !(Ktracks[n]->Is3D()) ) continue;
              //tracks3d.push_back(Ktracks[n]);
              
              const double L = Ktracks[n]->TotalLength();
              if(L > longestKTrack){
                longestKTrack = L;
                longestKTrack_ncells=Ktracks[n]->NCell();
                longestKTrack_energy=Ktracks[n]->TotalGeV();
                longestKTrack_nplanes=Ktracks[n]->ExtentPlane();
                TVector3 dir = Ktracks[n]->Dir();
                longestKTrack_phi=dir.Phi();
                longestKTrack_theta=dir.Theta();
                TVector3 start = Ktracks[n]->Start();
                TVector3 stop = Ktracks[n]->Stop();
                longestKTrack_startx = start.X();
                longestKTrack_starty = start.Y();
                longestKTrack_startz = start.Z();
                longestKTrack_stopx = stop.X();
                longestKTrack_stopy = stop.Y();
                longestKTrack_stopz = stop.Z();
                if( Ktracks[n]->Is3D() ) longestKTrack_is3D=1;
              }//for longest track
            }//loop all tracks
          }//has tracks
          
          //mip cells and reco energy for slicer
          double nCells = slice.NCell();
          int lengthPlanes = slice.ExtentPlane();
          if( slice.ExtentPlane()>192 ) lengthPlanes=192;
          
          double recoEnergy=0.;
          double nmip=0.;
          double Ntotal=0.;
          for( unsigned int icell=0; icell<slice.NCell(); ++ icell){
            const art::Ptr<rb::CellHit>& chit = slice.Cell(icell);
            const rb::RecoHit rhit(cal->MakeRecoHit(*chit, slice.W(chit.get())));       
            
            if( !rhit.IsCalibrated() ) continue;
            //std::cout<<"cell "<<icell<<" energy is "<<rhit.GeV()<<std::endl;
            recoEnergy += rhit.GeV();
            
            if( rhit.PECorr()>100. && rhit.PECorr()<245. ) nmip++;
            Ntotal++;
          }
          
          sl_runnum = evt.run();
          sl_subrunnum = evt.subRun();
          sl_evtnum = evt.id().event();
          sl_Nslices=nslice;
          
          NKtracks = Ktracks.size();
          KtrackLength = longestKTrack;
          KtrackEnergy = longestKTrack_energy;
          KtrackNcells = longestKTrack_ncells;
          KtrackEnergy_ratio = longestKTrack_energy/recoEnergy;
          KtrackNcells_ratio = longestKTrack_ncells/nCells;
          KtrackNplanes = longestKTrack_nplanes;
          KtrackPhi = longestKTrack_phi;
          KtrackTheta = longestKTrack_theta;
          KtrackStartX = longestKTrack_startx;
          KtrackStartY = longestKTrack_starty;
          KtrackStartZ = longestKTrack_startz;
          KtrackStopX = longestKTrack_stopx;
          KtrackStopY = longestKTrack_stopy;
          KtrackStopZ = longestKTrack_stopz;
          Ktrack3D = longestKTrack_is3D;
          
          sliceIndex = sliceIdx;
          Eslice = recoEnergy;
          Ncells = nCells;
          Nplanes = lengthPlanes;
          Nmip=nmip;
          Fmip=nmip/Ntotal;
          
          recoVx=vert[0]->GetX();
          recoVy=vert[0]->GetY();
          recoVz=vert[0]->GetZ();
          
          ////////////////////////////////
          //select most-energetic shower//
          ////////////////////////////////
          art::FindManyP<rb::Prong> fmProng3D(vert, evt, art::InputTag(fProngLabel,fInstLabel3D));
          std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
          
          int nprongs=SelectedProng.size();
          if( nprongs>20 ) nprongs=20;
          Nprongs=nprongs;
          std::vector<double> ProngEnergy;
          std::vector<double> ProngEnergyXView;
          std::vector<double> ProngEnergyYView;
          std::vector<double> ProngNCells;
          std::vector<double> ProngNCellsXView;
          std::vector<double> ProngNCellsYView;
          std::vector<double> ProngPhi;
          std::vector<double> ProngTheta;
          std::vector<double> ProngLength;
          std::vector<double> ProngStartX;
          std::vector<double> ProngStartY;
          std::vector<double> ProngStartZ;
          std::vector<double> ProngStopX;
          std::vector<double> ProngStopY;
          std::vector<double> ProngStopZ;
          std::vector<int> ProngNplanes;
          
          //energy balance between first two prongs
          double ratio=1.;
          //phi angle between first two prongs
          double dphi=0.;
          //index for most-energetic prong
          int Index1=-1;
          if( SelectedProng.size()>0 ){
            
            std::vector<int> prongIndex;
            std::vector<float> pcellPECorr;

            for(int ip=0; ip<nprongs; ++ip ){ 
              //reco direction
              TVector3 dir1 = SelectedProng[ip]->Dir();
              TVector3 start1 = SelectedProng[ip]->Start();
              TVector3 Pstop = start1 + (SelectedProng[ip]->TotalLength())*dir1;
              
              //reco momentum
              double Eprong=0.;     
              double EprongX=0.;
              double EprongY=0.;
              
              int npcells=0;;
              int nxpcells=0;
              int nypcells=0;
              
              for( unsigned int icell=0; icell<SelectedProng[ip]->NCell(); ++ icell){
                const art::Ptr<rb::CellHit>& chit = SelectedProng[ip]->Cell(icell);
                const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng[ip]->W(chit.get())));
                
                if( !rhit.IsCalibrated() ) continue;
                Eprong += rhit.GeV();
                npcells++;
                
                if( chit->View() == geo::kX ){
                  EprongX += rhit.GeV();
                  ++nxpcells;
                }
                if( chit->View() == geo::kY ){
                  EprongY += rhit.GeV();
                  ++nypcells;
                }
                prongIndex.push_back(ip);
                pcellPECorr.push_back(rhit.PECorr());   
              }//loop prong associated cells
              ProngTheta.push_back(dir1.Theta());
              ProngLength.push_back(SelectedProng[ip]->TotalLength());
              
              ProngStartX.push_back(start1.X());
              ProngStartY.push_back(start1.Y());
              ProngStartZ.push_back(start1.Z());
              ProngStopX.push_back(Pstop.X());
              ProngStopY.push_back(Pstop.Y());
              ProngStopZ.push_back(Pstop.Z());
              
              ProngEnergy.push_back(Eprong);
              ProngEnergyXView.push_back(EprongX);
              ProngEnergyYView.push_back(EprongY);
              ProngPhi.push_back(dir1.Phi());
              ProngNCells.push_back(npcells);
              ProngNCellsXView.push_back(nxpcells);
              ProngNCellsYView.push_back(nypcells);
 
              ProngNplanes.push_back(SelectedProng[ip]->ExtentPlane());
            }//loop prongs
            
            int Index2=-1;
            double MaxEnergy1=-999.;
            double MaxEnergy2=-999.;
            for( int i=0;i<nprongs; ++i ){
              if( MaxEnergy1<ProngEnergy[i] ){
                MaxEnergy1=ProngEnergy[i];
                Index1=i;
              }
            }//leading
            
            p1Phi=ProngPhi[Index1];
            p1Theta=ProngTheta[Index1];
            p1StartX=ProngStartX[Index1];
            p1StartY=ProngStartY[Index1];
            p1StartZ=ProngStartZ[Index1];
            p1StopX=ProngStopX[Index1];
            p1StopY=ProngStopY[Index1];
            p1StopZ=ProngStopZ[Index1];
            p1Length=ProngLength[Index1];
            
            p1Energy=ProngEnergy[Index1];
            p1EnergyXView=ProngEnergyXView[Index1];
            p1EnergyYView=ProngNCellsYView[Index1];
            p1NCells=ProngNCells[Index1];
            p1NCellsXView=ProngNCellsXView[Index1];
            p1NCellsYView=ProngNCellsYView[Index1];
            p1Nplanes=ProngNplanes[Index1];
            
            int p1_ncells=0;
            int p1_mipcells=0;
            std::vector<float> p1cellX;
            std::vector<float> p1cellY;
            std::vector<float> p1cellZ;
            std::vector<float> p1cellE;
            std::vector<int> p1cellPlane;
            std::vector<int> p1cellCell;
            std::vector<int> p1cellADC;
            std::vector<float> p1cellPE;
            std::vector<float> p1cellPECorr;
            std::vector<int> p1cellXView;
            std::vector<float> p1cellTNS;
            std::vector<float> p1cellW;
            
            for( unsigned int icell=0; icell<SelectedProng[Index1]->NCell(); ++ icell){
              const art::Ptr<rb::CellHit>& chit = SelectedProng[Index1]->Cell(icell);
              const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng[Index1]->W(chit.get())));
              double wx=SelectedProng[Index1]->W(chit.get());
              double wy=wx;
              if( !rhit.IsCalibrated() ) continue;
              ++p1_ncells;
              if( rhit.PECorr()>100. && rhit.PECorr()<240. ) ++p1_mipcells;
              
              if( chit->View()==geo::kX ){
                p1cellW.push_back(wx);
              }//XZ view
              else{
                p1cellW.push_back(wy);
              }//YZ view
              
              p1cellX.push_back(rhit.X());
              p1cellY.push_back(rhit.Y());
              p1cellZ.push_back(rhit.Z());
              p1cellE.push_back(rhit.GeV());
              
              p1cellADC.push_back(chit->ADC());
              p1cellPlane.push_back(chit->Plane());
              p1cellCell.push_back(chit->Cell());
              p1cellPECorr.push_back(rhit.PECorr());
              
            }//loop energetic prong associated cells
            //if( p1_ncells>500 ){
            //   std::cout<<"&&& crazy event with "<<p1_ncells<<" cells for leading prong &&&"<<std::endl;
            //   std::cout<<evt.id().event()<<" "<<sliceIdx<<" "<<nu.Nu().E()<<std::endl;
            // }
            if( p1_ncells>500 ) p1_ncells=500;
            prong1cellNcell=p1_ncells;
            int minplane=999;
            for( int ic=0; ic<p1_ncells; ++ic ){
              prong1cellX[ic]=p1cellX[ic];
              prong1cellY[ic]=p1cellY[ic];
              prong1cellZ[ic]=p1cellZ[ic];
              prong1cellE[ic]=p1cellE[ic];
              
              prong1cellADC[ic]=p1cellADC[ic];
              prong1cellPlane[ic]=p1cellPlane[ic];
              prong1cellCell[ic]=p1cellCell[ic];
              prong1cellPECorr[ic]=p1cellPECorr[ic];
              prong1cellW[ic]=p1cellW[ic];

              if( minplane>p1cellPlane[ic] ){
                 minplane=p1cellPlane[ic];
              }
            }//loop selected cells
            if( p1_ncells>0 ) p1Fmip=(float)p1_mipcells/p1_ncells;
            
            int Nplanes_p1=ProngNplanes[Index1];
 
            if( ProngEnergy[Index1]>0. ){
            std::vector<double> p1PlaneEnergy;
            p1PlaneEnergy.clear();
            double totalE_road20=0.;
            for( int ip=0; ip<Nplanes_p1; ++ip ){
              double eplane=0.;

              double cellX=0.;//energetic cell X in each plane
              double cellY=0.;//energetic cell Y in each plane
              double Ecell=0.;
              double EXcell=0.;
              double EYcell=0.;

              std::vector<double> Xxview;
              std::vector<double> Exview;
              std::vector<double> Yyview;
              std::vector<double> Eyview;
              Xxview.clear();
              Exview.clear();
              Yyview.clear();
              Eyview.clear();

              for( int ic=0; ic<p1_ncells; ++ic ){
                int Dplane=p1cellPlane[ic]-minplane;
                if( Dplane==ip ){
                  eplane += p1cellE[ic];

                  Ecell += p1cellE[ic];
                  EXcell += p1cellE[ic]*p1cellX[ic];
                  EYcell += p1cellE[ic]*p1cellY[ic];

                  if( p1cellPlane[ic]%2 != 0 ){
                    Xxview.push_back(p1cellX[ic]);
                    Exview.push_back(p1cellE[ic]);
                  }//in X view
                  else {
                    Yyview.push_back(p1cellY[ic]);
                    Eyview.push_back(p1cellE[ic]);
                  }//in Y view
                }
              }//loop leading-prong associated cells

              p1PlaneEnergy.push_back(eplane);

              if( Ecell>0.0 ){
                cellX = EXcell/Ecell;
                cellY = EYcell/Ecell;
              }
              int Nxviews=Xxview.size();
              int Nyviews=Yyview.size();
              if( Nxviews>0 ){
                for( int ixv=0; ixv<Nxviews; ++ixv ){
                  if( fabs(Xxview[ixv]-cellX) < 2.*2. ) totalE_road20+=Exview[ixv];
                }
              }//in X view
              if( Nyviews>0 ){
                for( int iyv=0; iyv<Nyviews; ++iyv ){
                  if( fabs(Yyview[iyv]-cellY) < 2.*2. ) totalE_road20+=Eyview[iyv];
                }
              }//in Y view

            }//loop leading-prong associated planes

            double efrac_6plane=-1.;
            double p1_Eplane[20]={0.};//energy deposit upto 20 planes
           
            int nplanes=p1PlaneEnergy.size();
            for( int ip=0; ip<nplanes; ++ip ){
              for( int j=0; j<20; ++j ){
                if( ip<=j ) p1_Eplane[j]+=p1PlaneEnergy[ip];
              }

              if( ip<nplanes-6 ){
                double e_frac6=(p1PlaneEnergy[ip]+p1PlaneEnergy[ip+1]+p1PlaneEnergy[ip+2]+p1PlaneEnergy[ip+3]+p1PlaneEnergy[ip+4]+p1PlaneEnergy[ip+5])/ProngEnergy[Index1];
                if( efrac_6plane<e_frac6 ) efrac_6plane=e_frac6;
              }
              else if( nplanes<=6 ) efrac_6plane=1.;
              
            }//end loop for p1PlaneEnergy

            double efrac_Eplane[20]={0.};
            for( int ip=0; ip<20; ++ip ){
              efrac_Eplane[ip] = p1_Eplane[ip]/ProngEnergy[Index1];
            }

            double efrac_plane2=0.;
            double efrac_plane3=0.;
            double efrac_plane4=0.;

            if( p1PlaneEnergy.size()>1 && ProngEnergy[Index1]>0. )
              efrac_plane2=p1PlaneEnergy[1]/ProngEnergy[Index1];
            if( p1PlaneEnergy.size()>2 && ProngEnergy[Index1]>0. )
              efrac_plane3=p1PlaneEnergy[2]/ProngEnergy[Index1];
            if( p1PlaneEnergy.size()>3 && ProngEnergy[Index1]>0. )
              efrac_plane4=p1PlaneEnergy[3]/ProngEnergy[Index1];

            emaxfrac_6p=efrac_6plane;
            efrac_10p=efrac_Eplane[9];
            efrac_p2=efrac_plane2;
            efrac_p3=efrac_plane3;
            efrac_p4=efrac_plane4;
            efrac_2sig=totalE_road20/ProngEnergy[Index1];

            }//leading-prong energy > 0

            //for sub-leading prong  
            if( nprongs>1 ){
              for( int i=0;i<nprongs; ++i ){
                if( i==Index1 ) continue;
                if( MaxEnergy2<ProngEnergy[i] ){
                  MaxEnergy2=ProngEnergy[i];
                  Index2=i;
                }
              }//sub-leading
              double Eden=ProngEnergy[Index1]+ProngEnergy[Index2];
              double Enum=ProngEnergy[Index1]-ProngEnergy[Index2];
              if( Eden>0.0 ) ratio=Enum/Eden;
              dphi=fabs(ProngPhi[Index1]-ProngPhi[Index2]);
              if( dphi>2.*TMath::Pi() ) dphi=dphi-2.*TMath::Pi();
              if( dphi>TMath::Pi() ) dphi=2.*TMath::Pi()-dphi;
 
              //fraction of mip hits for sub-leading prong
              if( prongIndex.size()>0 ){
                double nmips_p2=0.;
                double ntotal_p2=0.;
                for( unsigned int ic=0; ic<prongIndex.size(); ++ic ){
                  if( prongIndex[ic] != Index2 ) continue;
                    ntotal_p2 += 1.;
                  if( pcellPECorr[ic]>100. && pcellPECorr[ic]<240. )
                    nmips_p2 += 1.;  
                }//
                if( nmips_p2>0. ) p2Fmip = nmips_p2/ntotal_p2;
              }//

            }//has at least two prongs
            
            Ebalance=ratio;
            Dphi=dphi;
          }//has at least one prong
          
         //LID 
         art::FindOneP<slid::EventLID> fslLID(slices, evt, flidLabel);
         cet::maybe_ref<art::Ptr<slid::EventLID> const> rlid(fslLID.at(sliceIdx));
         art::Ptr<slid::EventLID> elid = rlid.ref();
         if( elid ){
           lid = elid->Value();
         }
          
        }//has only one vertex
      }//has a valid reco vertex    
      
      _otree->Fill();
    }//end loop of slice
    
    return;
  }//end analyze
  
  void Xeff::endJob()
  {
    art::ServiceHandle<art::TFileService> tfs;
    TH1* hPOT = tfs->make<TH1F>("hTotalPOT", ";; POT", 1, 0, 1);
    hPOT->Fill(.5, fTotalPOT);
    TH1* hSPILL = tfs->make<TH1F>("hTotalSpill",";; SPILL", 1, 0, 1);
    hSPILL->Fill(.5, fTotalSpill);
  }
  
}//end namespace

namespace ncs{
  
  DEFINE_ART_MODULE(Xeff);
  
}