ROCKMRE_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 "canvas/Persistency/Common/FindOneP.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 ROCKMRE : public art::EDAnalyzer {
  public:
    explicit ROCKMRE(fhicl::ParameterSet const& pset);
    virtual ~ROCKMRE();
    
    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 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 fDataSpillLabel; //label for data spill
    std::string flidLabel;       //label for event LID
    int isMC; //mode for MC

    //save event info to ntuples
    int runnum;
    int subrunnum;
    int evtnum;
    int Nslices;
    double spillPOT;
    int goodSpill;

    int Ndeaddcm;
  
    TTree *_otree;
    
    //true info for MC
    double trueE;
    int ccnc;
    int PDG;
    int origPDG;
    int intType;
    double Vx;
    double Vy;
    double Vz;

    //TTree *_evttree;
    int sl_runnum;
    int sl_subrunnum;
    int sl_evtnum;
    int sl_Nslices;
    int sl_goodSpill;
    double sl_spillPOT;
    int sl_Ndeaddcm;

    /*** Vertex ***/
    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 recoE;//reco energy
    float PlaneEnergy[192];//set the maximal plane # to 192 
    int MinPlane;//min plane index
    int PlaneNcells[192];//# of cells in each plane
    int sliceIndex;
    
    /*** mip information ***/
    float Nmip;//# of mip cells
    float Fmip;//# of mip cells over Ncells
    
    /*** cell information ***/
    int slcellNcell;//# of good cells
    float slcellX[500];
    float slcellY[500];
    float slcellZ[500];
    float slcellE[500];
    int slcellADC[500];
    int slcellPlane[500];
    int slcellCell[500];
    float slcellPECorr[500];
    float slcellW[500];

    /*** prong information ***/
    int Nprongs;
    float Ebalance;
    float Dphi;

    float prongEnergy[20];
    float prongPhi[20];
    float prongTheta[20];
    float prongStartX[20];
    float prongStartY[20];
    float prongStartZ[20];
    float prongLength[20];
    float prongStopX[20];
    float prongStopY[20];
    float prongStopZ[20];
    float prongEnergyXView[20];//energy for XView cells
    float prongEnergyYView[20];//energy for YView cells
    int prongNCells[20];//# of cells
    int prongNCellsXView[20];//# of cells for XView
    int prongNCellsYView[20];//# of cells for YView
    int prongNplanes[20];

    //leading prong
    int prong1cellNcell;
    float prong1Fmip;
    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];

    //2D prongs
    int Nprongs2D;
    int prong2DInXView[20];//prong in X View or not
    float prongEnergy2D[20];//energy
    int prongNCells2D[20];//# of cells
    float prongPhi2D[20];
    float prongTheta2D[20];
    float prongStartX2D[20];
    float prongStartY2D[20];
    float prongStartZ2D[20];
    float prongLength2D[20];
    float prongStopX2D[20];
    float prongStopY2D[20];
    float prongStopZ2D[20];
    float Ebalance2D;

    float lid;
  };
}
 
namespace ncs{
  
  //......................................................................
  ROCKMRE::ROCKMRE(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
  {
    mf::LogInfo("ROCKMRE")<<__PRETTY_FUNCTION__<<"\n";
    reconfigure(pset);
  }
  
  //......................................................................
  ROCKMRE::~ROCKMRE()
  {
  }
  
  //......................................................................
  void ROCKMRE::reconfigure(const fhicl::ParameterSet& pset)
  {
    fSliceLabel   = pset.get< std::string >("SliceLabel");
    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");

    isMC          = pset.get<int>("MCevts");
    fDataSpillLabel = pset.get<std::string> ("DataSpillLabel");
  }

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

    fTotalPOT=0.;
    fTotalSpill=0;

    _otree = tfs->make<TTree>("ROCKMRE","ROCKMRE particle");
    
    _otree->Branch("trueE",&trueE,"trueE/F");
    _otree->Branch("ccnc",&ccnc,"ccnc/I");
    _otree->Branch("PDG",&PDG,"PDG/I");
    _otree->Branch("origPDG",&origPDG,"origPDG/I");
    _otree->Branch("intType",&intType,"intType/I");
    _otree->Branch("Vx",&Vx,"Vx/F");
    _otree->Branch("Vy",&Vy,"Vy/F");
    _otree->Branch("Vz",&Vz,"Vz/F");
    
    _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_goodSpill",&sl_goodSpill,"sl_goodSpill/I");
    _otree->Branch("sl_spillPOT",&sl_spillPOT,"sl_spillPOT/D");

    _otree->Branch("sl_Ndeaddcm",&sl_Ndeaddcm,"sl_Ndeaddcm/I");

    _otree->Branch("recoVx",&recoVx,"recoVx/F");
    _otree->Branch("recoVy",&recoVy,"recoVy/F");
    _otree->Branch("recoVz",&recoVz,"recoVz/F");

    _otree->Branch("slcellNcell",&slcellNcell,"slcellNcell/I");
    _otree->Branch("slcellX",&slcellX,"slcellX[slcellNcell]/F");
    _otree->Branch("slcellY",&slcellY,"slcellY[slcellNcell]/F");
    _otree->Branch("slcellZ",&slcellZ,"slcellZ[slcellNcell]/F");
    _otree->Branch("slcellE",&slcellE,"slcellE[slcellNcell]/F");
    _otree->Branch("slcellADC",&slcellADC,"slcellADC[slcellNcell]/I");
    _otree->Branch("slcellPlane",&slcellPlane,"slcellPlane[slcellNcell]/I");
    _otree->Branch("slcellCell",&slcellCell,"slcellCell[slcellNcell]/I");
    _otree->Branch("slcellW",&slcellW,"slcellW[slcellNcell]/F");
    _otree->Branch("slcellPECorr",&slcellPECorr,"slcellPECorr[slcellNcell]/F");

    _otree->Branch("Ncells",&Ncells,"Ncells/I");
    _otree->Branch("Nplanes",&Nplanes,"Nplanes/I");
    _otree->Branch("recoE",&recoE,"recoE/F");
    _otree->Branch("MinPlane",&MinPlane,"MinPlane/I");
    _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("prongLength",&prongLength,"prongLength[Nprongs]/F");
    _otree->Branch("prongStopX",&prongStopX,"prongStopX[Nprongs]/F");
    _otree->Branch("prongStopY",&prongStopY,"prongStopY[Nprongs]/F");
    _otree->Branch("prongStopZ",&prongStopZ,"prongStopZ[Nprongs]/F");

    _otree->Branch("prongEnergy",&prongEnergy,"prongEnergy[Nprongs]/F");
    _otree->Branch("prongPhi",&prongPhi,"prongPhi[Nprongs]/F");
    _otree->Branch("prongTheta",&prongTheta,"prongTheta[Nprongs]/F");
    _otree->Branch("prongStartX",&prongStartX,"prongStartX[Nprongs]/F");
    _otree->Branch("prongStartY",&prongStartY,"prongStartY[Nprongs]/F");
    _otree->Branch("prongStartZ",&prongStartZ,"prongStartZ[Nprongs]/F");
    _otree->Branch("prongEnergyXView",&prongEnergyXView,"prongEnergyXView[Nprongs]/F");
    _otree->Branch("prongEnergyYView",&prongEnergyYView,"prongEnergyYView[Nprongs]/F");
    _otree->Branch("prongNCells",&prongNCells,"prongNCells[Nprongs]/I");
    _otree->Branch("prongNCellsXView",&prongNCellsXView,"prongNCellsXView[Nprongs]/I");
    _otree->Branch("prongNCellsYView",&prongNCellsYView,"prongNCellsYView[Nprongs]/I");
    _otree->Branch("prongNplanes",&prongNplanes,"prongNplanes[Nprongs]/I");

    _otree->Branch("prong1cellNcell",&prong1cellNcell,"prong1cellNcell/I");
    _otree->Branch("prong1Fmip",&prong1Fmip,"prong1Fmip/F");
    _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("Nprongs2D",&Nprongs2D,"Nprongs2D/I");
    _otree->Branch("prong2DInXView",&prong2DInXView,"prong2DInXView[Nprongs2D]/I");
    _otree->Branch("prongEnergy2D",&prongEnergy2D,"prongEnergy2D[Nprongs2D]/F");
    _otree->Branch("prongNCells2D",&prongNCells2D,"prongNCells2D[Nprongs2D]/I");
    _otree->Branch("prongLength2D",&prongLength2D,"prongLength2D[Nprongs2D]/F");
    _otree->Branch("prongStopX2D",&prongStopX2D,"prongStopX2D[Nprongs2D]/F");
    _otree->Branch("prongStopY2D",&prongStopY2D,"prongStopY2D[Nprongs2D]/F");
    _otree->Branch("prongStopZ2D",&prongStopZ2D,"prongStopZ2D[Nprongs2D]/F");
    _otree->Branch("prongPhi2D",&prongPhi2D,"prongPhi2D[Nprongs2D]/F");
    _otree->Branch("prongTheta2D",&prongTheta2D,"prongTheta2D[Nprongs2D]/F");
    _otree->Branch("prongStartX2D",&prongStartX2D,"prongStartX2D[Nprongs2D]/F");
    _otree->Branch("prongStartY2D",&prongStartY2D,"prongStartY2D[Nprongs2D]/F");
    _otree->Branch("prongStartZ2D",&prongStartZ2D,"prongStartZ2D[Nprongs2D]/F"); 
    _otree->Branch("Ebalance2D",&Ebalance2D,"Ebalance2D/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("goodSpill",&goodSpill,"goodSpill/I");

    _btree->Branch("Ndeaddcm",&Ndeaddcm,"Ndeaddcm/I");

  }

  //......................................................................
  void ROCKMRE::endSubRun(const art::SubRun& sr)
  {
    art::Handle<sumdata::POTSum> pot;
    if( isMC ) sr.getByLabel(fMCFluxLabel, pot);
    else sr.getByLabel(fDataSpillLabel, 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 ROCKMRE::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 ROCKMRE::analyze(const art::Event& evt) 
  {

    lid=-1.;

    sl_runnum=-1;
    sl_subrunnum=-1;
    sl_evtnum=-1;
    sl_Nslices=-1;
    sl_goodSpill=-1;
    sl_spillPOT=-1.;
    sl_Ndeaddcm=-1;
    
    runnum=-1;
    subrunnum=-1;
    evtnum=-1;
    Nslices=-1;
    spillPOT=0.;
    goodSpill=-1;
    Ndeaddcm=-1;
    
    recoVx=-999.;
    recoVy=-999.;
    recoVz=-999.;

    Ncells=0;
    Nplanes=0;
    recoE=-1.;
    MinPlane=-1;
    sliceIndex=-1;
    for( int i=0; i<192; ++i ){
      PlaneEnergy[i] = 0.;
      PlaneNcells[i] = 0;
    }
    
    Nmip=-1.;
    Fmip=-1.;
  
    Nprongs=0;
    Ebalance=0.;
    Dphi=-1.;

    for( int i=0; i<20; ++i ){
      prongEnergy[i]=0.;
      prongPhi[i]=-1.;
      prongTheta[i]=-1.;
      prongStartX[i]=-999.;
      prongStartY[i]=-999.;
      prongStartZ[i]=-999.;
      prongStopX[i]=-999.;
      prongStopY[i]=-999.;
      prongStopZ[i]=-999.;
      prongLength[i]=-1.;
      prongEnergyXView[i]=0.;
      prongEnergyYView[i]=0.;
      prongNCells[i]=0;
      prongNCellsXView[i]=0;
      prongNCellsYView[i]=0;
      prongNplanes[i]=-1;

      prong2DInXView[i]=0;
      prongEnergy2D[i]=0.;
      prongNCells2D[i]=0;
      prongLength2D[i]=-1.;
      prongPhi2D[i]=-1.;
      prongTheta2D[i]=-1.;
      prongStartX2D[i]=-999.;
      prongStartY2D[i]=-999.;
      prongStartZ2D[i]=-999.;
      prongStopX2D[i]=-999.;
      prongStopY2D[i]=-999.;
      prongStopZ2D[i]=-999.;
    }
    Nprongs2D=0;
    Ebalance2D=0.;
    
    slcellNcell=0;
    prong1cellNcell=0;
    prong1Fmip=-1.;
    for( int i=0; i<500; ++i ){
      slcellX[i]=-999.;
      slcellY[i]=-999.;
      slcellZ[i]=-999.;
      slcellE[i]=-999.;
      slcellADC[i]=0;
      slcellPlane[i]=-1;
      slcellCell[i]=-1;
      slcellW[i]=-1.;
      slcellPECorr[i]=-1.;

      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;
   
    //get data spill info
    art::Handle<sumdata::SpillData> spillPot;
    if( isMC ) evt.getByLabel(fMCFluxLabel, spillPot);
    else evt.getByLabel(fDataSpillLabel, spillPot);
   
    if( spillPot.failedToGet() ) return;
    spillPOT = spillPot->spillpot;
    if( isMC ) spillPOT=spillPOT/1.0e+12;
    if( spillPot->goodbeam ) goodSpill=1; 
    //std::cout<<"beam info: "<<goodSpill<<" "<<spillPOT<<std::endl;
    
    //event info
    runnum = evt.run();
    subrunnum = evt.subRun();
    evtnum = evt.id().event();

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

    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=chits0->size(); 
    for(int i = 0; i < sq_ncells; ++i){
      art::Ptr< rb::CellHit> hit(chits0, 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
    
    //find dead dcms
    int Ndbdcm_intime[14]={0};
    int Ndbdcm_outtime[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 ){
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[0]++;
               else Ndbdcm_outtime[0]++;
             }//dcm 1
             else{
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[1]++;
               else Ndbdcm_outtime[1]++;
             }//dcm 2
           }//XZ view
           else {
             if( sq_cell[ic]>31 ){
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[2]++;
               else Ndbdcm_outtime[2]++;
             }//dcm 3
             else {
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[3]++;
               else Ndbdcm_outtime[3]++;
             }//dcm 4
           }//YZ view
         }//DB1
         else if( sq_plane[ic]<128 ){
           if( sq_xview[ic]==1 ){
             if( sq_cell[ic]>63 ){
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[4]++;
               else Ndbdcm_outtime[4]++;
             }//dcm 1
             else{
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[5]++;
               else Ndbdcm_outtime[5]++;
             }//dcm 2
           }//XZ view
           else {
              if( sq_cell[ic]>31 ){
                if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[6]++;
                else Ndbdcm_outtime[6]++;
              }//dcm 3
              else {
                if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[7]++;
                else Ndbdcm_outtime[7]++;
              }//dcm 4
           }//YZ view 
         }//DB2
         else if( sq_plane[ic]<192 ){
           if( sq_xview[ic]==1 ){
             if( sq_cell[ic]>63 ){
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[8]++;
               else Ndbdcm_outtime[8]++;
             }//dcm 1
             else{
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[9]++;
               else Ndbdcm_outtime[9]++;
             }//dcm 2
           }//XZ view
           else {
             if( sq_cell[ic]>31 ){
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[10]++;
               else Ndbdcm_outtime[10]++;
             }//dcm 3
             else {
               if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[11]++;
               else Ndbdcm_outtime[11]++;
             }//dcm 4
           }//YZ view
         }//DB3
         else {
           if( sq_xview[ic]==1 ){
             if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[12]++;
             else Ndbdcm_outtime[12]++;
           }//XZ view
           else {
             if( sq_tns[ic]>218. && sq_tns[ic]<228. ) Ndbdcm_intime[13]++;
             else Ndbdcm_outtime[13]++;
           }//YZ view
         }//MC
    }//loop all cells
  
    int ndeaddcms=0;
    std::vector<std::string> dcmlocation;
    dcmlocation.clear();
    for( int id=0; id<14; ++id ){
      if( Ndbdcm_outtime[id]==0 && Ndbdcm_intime[id]==0 ){ 
        ndeaddcms++;
      }
    }//loop all 14 dcms associated cells
    
    Ndeaddcm=ndeaddcms;

    //get reco slicer info  
    art::Handle< std::vector<rb::Cluster> > slices;
    evt.getByLabel(fSliceLabel, slices);
    art::FindManyP<rb::Vertex> fmVertex(slices, evt, fVertexLabel);
    
    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 missing dcms: "<<ndeaddcms<<std::endl;
    //std::cout<<"# of slices: "<<nslice<<std::endl;

    _btree->Fill();
    
    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){

      const rb::Cluster& slice = (*slices)[sliceIdx];
      if(slice.IsNoise()) continue;
      
      if (!(fmVertex.isValid())) continue;
      std::vector<art::Ptr<rb::Vertex>> vert = fmVertex.at(sliceIdx);
      
      if (vert.size() != 1) continue;
      
     //std::cout<<"having vertex "<<std::endl;

      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;
        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;
      
      sl_Ndeaddcm=Ndeaddcm;
      if( !isMC ){
        sl_goodSpill=goodSpill; 
        sl_spillPOT = spillPOT;
      }
      
     //std::cout<<"# of cells "<<nCells<<std::endl;
     if( nCells>300 ) 
       std::cout<<"xbbu "<<nCells<<" "<<evt.run()<<" "<<evt.subRun()<<" "<<evt.id().event()<<" "<<sliceIdx<<std::endl;
     if( nCells<10 || nCells>500 ) continue;

/*
     if( isMC ){
         art::Handle<std::vector<rb::CellHit> > chits0;
         evt.getByLabel(fCellHitLabel, chits0);
         art::ServiceHandle<cheat::BackTracker> bt;
         art::PtrVector<rb::CellHit> allhits;
        for(size_t h = 0; h < chits0->size(); ++h)
          allhits.push_back(art::Ptr<rb::CellHit>(chits0, 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;
        }
        art::Ptr<simb::MCTruth>& truth = iacts[0].neutrinoInt;
        art::PtrVector<simb::MCTruth> pv;
        pv.push_back(truth);
        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();

        trueE   = nu.Nu().E();
        ccnc    = nu.CCNC();
        PDG     = nu.Nu().PdgCode();
        origPDG = flux.fntype;
        intType = nu.InteractionType();
        Vx=nu.Nu().Vx();
        Vy=nu.Nu().Vy();
        Vz=nu.Nu().Vz();
      }//for MC only
*/

      sliceIndex = sliceIdx;
      recoE = recoEnergy;
      Ncells = nCells;
      Nplanes = lengthPlanes;
      MinPlane = slice.MinPlane();
      Nmip=nmip;
      Fmip=nmip/Ntotal;
      
      recoVx=vert[0]->GetX();
      recoVy=vert[0]->GetY();
      recoVz=vert[0]->GetZ();
      
      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();
      }

      //reco prongs
      art::FindManyP<rb::Prong> fmProng3D(vert, evt, art::InputTag(fProngLabel,fInstLabel3D));
      std::vector<art::Ptr<rb::Prong>> SelectedProng = fmProng3D.at(0);
      
      //std::cout<<"# of prongs "<<SelectedProng.size()<<std::endl;
      if( SelectedProng.size()<1 ) continue;
      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> ProngPhi;


      for(int ip=0; ip<nprongs; ++ip ){ 
          //reco direction
          TVector3 dir1 = SelectedProng[ip]->Dir();
          //recoEta[ip]=dir1.Eta();
          prongPhi[ip]=dir1.Phi();
          prongTheta[ip]=dir1.Theta();
          
          TVector3 start1 = SelectedProng[ip]->Start();
          prongStartX[ip] = start1.X();
          prongStartY[ip] = start1.Y();
          prongStartZ[ip] = start1.Z();
          
          prongLength[ip] = SelectedProng[ip]->TotalLength();
          
          TVector3 Pstop = start1 + (SelectedProng[ip]->TotalLength())*dir1;
          prongStopX[ip] = Pstop.X();
          prongStopY[ip] = Pstop.Y();
          prongStopZ[ip] = Pstop.Z();
          
          prongNplanes[ip] = SelectedProng[ip]->ExtentPlane();

          //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;
            }
            
          }//loop prong associated cells
          ProngEnergy.push_back(Eprong);
          ProngEnergyXView.push_back(EprongX);
          ProngEnergyYView.push_back(EprongY);
          ProngPhi.push_back(dir1.Phi());
          prongEnergy[ip]=Eprong;
          prongEnergyXView[ip]=EprongX;
          prongEnergyYView[ip]=EprongY;
          prongNCells[ip]=npcells;
          prongNCellsXView[ip]=nxpcells;
          prongNCellsYView[ip]=nypcells;
        }//loop prongs
      
      int Index1=-1;
      double MaxEnergy1=-999.;
      if( nprongs>0 ){
        for( int i=0;i<nprongs; ++i ){
          if( MaxEnergy1<ProngEnergy[i] ){
            MaxEnergy1=ProngEnergy[i];
            Index1=i;
          }
        }//leading
      }

      if( Index1>=0 ){

        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
        prong1cellNcell=p1_ncells;
        if(p1_ncells>0 ) prong1Fmip = (float)p1_mipcells/p1_ncells;
        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];
        }//loop selected cells
      }//for most energecti prong
      
      //slice associated cell information
      int ncells_slice=0.;
      std::vector<float> slicecellX;
      std::vector<float> slicecellY;
      std::vector<float> slicecellZ;
      std::vector<float> slicecellE;
      std::vector<int> slicecellPlane;
      std::vector<int> slicecellCell;
      std::vector<int> slicecellADC;
      std::vector<float> slicecellPE;
      std::vector<float> slicecellPECorr;
      std::vector<int> slicecellXView;
      std::vector<float> slicecellTNS;
      std::vector<float> slicecellW; 
      std::vector<int> slicecellNoise;

      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;

        ncells_slice++;

             slicecellW.push_back(slice.W(chit.get()));

          slicecellX.push_back(rhit.X());
          slicecellY.push_back(rhit.Y());
          slicecellZ.push_back(rhit.Z());
          slicecellE.push_back(rhit.GeV());

          slicecellADC.push_back(chit->ADC());
          slicecellPlane.push_back(chit->Plane());
          slicecellCell.push_back(chit->Cell());
          slicecellPECorr.push_back(rhit.PECorr());

      }

      slcellNcell=ncells_slice;

      for( int ic=0; ic<ncells_slice; ++ic ){

          slcellX[ic]=slicecellX[ic];
          slcellY[ic]=slicecellY[ic];
          slcellZ[ic]=slicecellZ[ic];
          slcellE[ic]=slicecellE[ic];

          slcellADC[ic]=slicecellADC[ic];
          slcellPlane[ic]=slicecellPlane[ic];
          slcellCell[ic]=slicecellCell[ic];
          slcellPECorr[ic]=slicecellPECorr[ic];
          slcellW[ic]=slicecellW[ic];

        }//loop selected cells
      
      //study 2D prongs
      art::FindManyP<rb::Prong> fmProng2D(vert, evt, art::InputTag(fProngLabel,fInstLabel2D));
      std::vector<art::Ptr<rb::Prong>> SelectedProng2D = fmProng2D.at(0);
      
      std::vector<int> Prong2DInXView;
      std::vector<double> ProngEnergy2D;
      int nprongs2D=SelectedProng2D.size();
      if( nprongs2D>20 ) nprongs2D=20;
      Nprongs2D = nprongs2D;
      if( nprongs2D>0 ){
      for(int ip=0; ip<nprongs2D; ++ip ){
        prongNCells2D[ip]=SelectedProng2D[ip]->NCell();
        int prongview=SelectedProng2D[ip]->View();
        
        prong2DInXView[ip]=prongview;//0 for XZ view 
        prongLength2D[ip]=SelectedProng2D[ip]->TotalLength();
        
        TVector3 dir1 = SelectedProng2D[ip]->Dir();
        prongPhi2D[ip]=dir1.Phi();
        prongTheta2D[ip]=dir1.Theta();
        
        TVector3 start1 = SelectedProng2D[ip]->Start();
        prongStartX2D[ip] = start1.X();
        prongStartY2D[ip] = start1.Y();
        prongStartZ2D[ip] = start1.Z();
        
        TVector3 Pstop = start1 + (SelectedProng2D[ip]->TotalLength())*dir1;
        prongStopX2D[ip] = Pstop.X();
        prongStopY2D[ip] = Pstop.Y();
        prongStopZ2D[ip] = Pstop.Z();
        
        Prong2DInXView.push_back(prongview);
        //reco momentum
        double Eprong=0.;
        for( unsigned int icell=0; icell<SelectedProng2D[ip]->NCell(); ++ icell){
          const art::Ptr<rb::CellHit>& chit = SelectedProng2D[ip]->Cell(icell);
          const TVector3 mean = SelectedProng2D[ip]->MeanXYZ();
          const rb::RecoHit rhit(cal->MakeRecoHit(*chit, SelectedProng2D[ip]->W(chit.get())));

          if( !rhit.IsCalibrated() ) continue;
          Eprong += rhit.GeV();
        }
        
        prongEnergy2D[ip]=Eprong;
        ProngEnergy2D.push_back(Eprong);
      }
      
      int index3D=-1;
      double max3D=-99.;
      for( int i=0; i<nprongs; ++i ){
        if( max3D<ProngEnergy[i] ){
          max3D=ProngEnergy[i];
          index3D=i;
        }
      }//loop 3D prongs
      int index2D=-1;
      double max2D=-99.;
      for( int i=0;i<nprongs2D; ++i ){
        if( max2D<ProngEnergy2D[i] ){
          max2D=ProngEnergy2D[i];
          index2D=i;
        }
      }//loop 2D prongs
      double ebalance2D=1.;
      if( nprongs>0 && nprongs2D>0 ){
        double Eden=0.;
        double Enum=0.;
        if( Prong2DInXView[index2D]==0 ){//XZ view
          Eden=ProngEnergy2D[index2D]+ProngEnergyXView[index3D];
          Enum=ProngEnergy2D[index2D]-ProngEnergyXView[index3D];
        }
        else {//YZ view
          Eden=ProngEnergy2D[index2D]+ProngEnergyYView[index3D];
          Enum=ProngEnergy2D[index2D]-ProngEnergyYView[index3D];
        }
        if( Eden>0. ) ebalance2D=fabs(Enum)/Eden;
      }
      Ebalance2D=ebalance2D;
    }//having 2D prongs      

      _otree->Fill();
      
      //_evttree->Fill();
    }//end loop of slice
    
    //_evttree->Fill();
    return;
  }//end analyze
  
 void ROCKMRE::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(ROCKMRE);

}