Xbeam_module.cc

Go to the documentation of this file.

///
/// \brief 
/// \author Xuebing Bu - xbbu@fnal.gov
///
#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 "TProfile.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 Xbeam : public art::EDAnalyzer {
  public:
    explicit Xbeam(fhicl::ParameterSet const& pset);
    virtual ~Xbeam();
    
    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;

  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 fMCFluxLabel;    ///label for generator information
    std::string fCellHitLabel;   ///label for cell hits

    TProfile *slEff_Etrue;

    TTree *_otree;

    /*** true info ***/
    float nu_Vx;//true vertex X
    float nu_Vy;//true vertex Y
    float nu_Vz;//true vertex Z
    float nu_truePt;//true Pt
    float nu_trueE;//true E
    float nu_trueTheta;//true theta
    int nu_ccnc;
    int nu_PDG;
    int nu_origPDG;
    int nu_orig_parentPDG;
    int nu_mode;
    bool nu_isQE;
    int nu_intType;
    int nu_hitNucl;
    float nu_HadX;
    float nu_HadY;
    float nu_HadW;
    float nu_qsqr;
    float nu_OscWt; 
    int meson_PDG;
    float meson_Px;
    float meson_Py;
    float meson_Pz;
    int lep_PDG;
    float lep_Px;
    float lep_Py;
    float lep_Pz;
    float lep_E;

    /*** slicer information ***/
    int Nslices;
    int hasSlice;

  };
}

namespace ncs{
  
  //......................................................................
  Xbeam::Xbeam(fhicl::ParameterSet const& pset)
    : EDAnalyzer(pset)
  {
    mf::LogInfo("Xbeam")<<__PRETTY_FUNCTION__<<"\n";
    reconfigure(pset);
  }
  
  //......................................................................
  Xbeam::~Xbeam()
  {
  }
  
  //......................................................................
  void Xbeam::reconfigure(const fhicl::ParameterSet& pset)
  {
    fSliceLabel   = pset.get< std::string >("SliceLabel");
    fMCFluxLabel  = pset.get< std::string >("MCFluxLabel");
    fCellHitLabel = pset.get< std::string >("CellHitLabel");
  }

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

    fTotalPOT=0.;
    fTotalSpill=0;

    //histo
    slEff_Etrue=tfs->make<TProfile>("slEff_Etrue","",20,0.,10.,0,1.2);
  
    //personal tuples
    _otree = tfs->make<TTree>("Xbeam","Xbeam particle");
    
    _otree->Branch("nu_Vx",&nu_Vx,"nu_Vx/F");
    _otree->Branch("nu_Vy",&nu_Vy,"nu_Vy/F");
    _otree->Branch("nu_Vz",&nu_Vz,"nu_Vz/F");
    _otree->Branch("nu_truePt",&nu_truePt,"nu_truePt/F");
    _otree->Branch("nu_trueE",&nu_trueE,"nu_trueE/F");
    _otree->Branch("nu_trueTheta",&nu_trueTheta,"nu_trueTheta/F");
    _otree->Branch("nu_ccnc",&nu_ccnc,"nu_ccnc/I");
    _otree->Branch("nu_PDG",&nu_PDG,"nu_PDG/I"); 
    _otree->Branch("nu_orig_parentPDG",&nu_orig_parentPDG,"nu_orig_parentPDG/I"); 
    _otree->Branch("nu_origPDG",&nu_origPDG,"nu_origPDG/I"); 
    _otree->Branch("nu_mode",&nu_mode,"nu_mode/I");
    _otree->Branch("nu_isQE",&nu_isQE,"nu_isQE/B");
    _otree->Branch("nu_intType",&nu_intType,"nu_intType/I");
    _otree->Branch("nu_hitNucl",&nu_hitNucl,"nu_hitNucl/I");
    _otree->Branch("nu_HadX",&nu_HadX,"nu_HadX/F");
    _otree->Branch("nu_HadY",&nu_HadY,"nu_HadY/F");
    _otree->Branch("nu_HadW",&nu_HadW,"nu_HadW/F");
    _otree->Branch("nu_qsqr",&nu_qsqr,"nu_qsqr/F");
    _otree->Branch("nu_OscWt",&nu_OscWt,"nu_OscWt/F");
    _otree->Branch("meson_PDG",&meson_PDG,"meson_PDG/I");
    _otree->Branch("meson_Px",&meson_Px,"meson_Px/F");
    _otree->Branch("meson_Py",&meson_Py,"meson_Py/F");
    _otree->Branch("meson_Pz",&meson_Pz,"meson_Pz/F");
    _otree->Branch("lep_PDG",&lep_PDG,"lep_PDG/I");
    _otree->Branch("lep_Px",&lep_Px,"lep_Px/F");
    _otree->Branch("lep_Py",&lep_Py,"lep_Py/F");
    _otree->Branch("lep_Pz",&lep_Pz,"lep_Pz/F");
    _otree->Branch("lep_E",&lep_E,"lep_E/F");

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

  }

  //......................................................................
  void Xbeam::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 Xbeam::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 Xbeam::analyze(const art::Event& evt) 
  {
      nu_Vx = -999.;      
      nu_Vy = -999.;      
      nu_Vz = -999.;      
      nu_truePt = -1.;
      nu_trueE = -1.;
      nu_trueTheta = -1.;

      nu_ccnc=-1;
      nu_PDG=-1;
      nu_origPDG=-1;
      nu_orig_parentPDG=-1;
      nu_mode=-1;
      nu_isQE=false;
      nu_intType=-1;
      nu_hitNucl=-1;
      
      nu_HadX=-1.;
      nu_HadY=-1.;
      nu_HadW=-1.;
      nu_qsqr=-1.;
      nu_OscWt=-1.;

      meson_PDG=-1;
      meson_Px=-1.;
      meson_Py=-1.;
      meson_Pz=-1.;
     
      lep_PDG=-1;
      lep_Px=-1.;
      lep_Py=-1.;
      lep_Pz=-1.;
      lep_E=-1.;

    Nslices=-1;
    hasSlice=-1;
      
    art::ServiceHandle<calib::Calibrator> cal;
    art::ServiceHandle<geo::Geometry> geom;
    //////////////
    //Truth info//
    //////////////
    art::Handle<std::vector<simb::MCTruth> > mclist;
    evt.getByLabel(fMCFluxLabel,mclist);
    art::Handle< std::vector<simb::MCFlux> > fllist;
    evt.getByLabel(fMCFluxLabel, fllist);
    //# of interactions
    //std::cout<<"# of interactions: "<<mclist->size()<<std::endl;
    if( mclist->size()!=1 ){
      std::cout<<"ABNORMAL! "<<mclist->size()<<" interactions."<<std::endl;
      return;
    }

      art::Ptr<simb::MCTruth> mctruth(mclist,0);
      const simb::MCNeutrino& nu = mctruth->GetNeutrino();

      nu_trueE   = nu.Nu().E();
      nu_PDG     = nu.Nu().PdgCode();
      nu_Vx=nu.Nu().Vx();
      nu_Vy=nu.Nu().Vy();
      nu_Vz=nu.Nu().Vz();

      nu_truePt  = nu.Pt();
      nu_trueTheta = nu.Theta();
      nu_ccnc    = nu.CCNC();
      nu_mode    = nu.Mode();
      nu_isQE    = nu.InteractionType() == simb::kQE ||
        nu.InteractionType() == simb::kCCQE;
      nu_intType = nu.InteractionType();
      nu_hitNucl = nu.HitNuc();
      nu_HadX    = nu.X();
      nu_HadY    = nu.Y();
      nu_HadW    = nu.W();
      nu_qsqr    = nu.QSqr();

      lep_PDG = nu.Lepton().PdgCode();
      lep_Px = nu.Lepton().Px();
      lep_Py = nu.Lepton().Py();
      lep_Pz = nu.Lepton().Pz();
      lep_E = nu.Lepton().E();


      art::Ptr<simb::MCFlux> flux(fllist,0);

      nu_origPDG = flux->fntype;
      nu_orig_parentPDG = flux->fptype;
      meson_PDG = flux->ftptype;
      meson_Px = flux->ftpx;
      meson_Py = flux->ftpy;
      meson_Pz = flux->ftpz;

      //std::cout<<i_intx<<" "<<nu.CCNC()<<" "<<flux->fntype<<" "<<nu.Nu().PdgCode()<<std::endl;
      //std::cout<<flux->fntype<<" "<<flux->fptype<<" "<<flux->ftptype<<std::endl;
      //std::cout<<"meson is "<<flux->ftptype<<std::endl;
    //CC nue only
    if( nu.CCNC() != 0 ) return;
    if( nu.Nu().PdgCode() != 12 ) return;
    if( flux->fntype != 12 ) return;

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

    bool matchSlice=false;
 
    for(unsigned int sliceIdx = 0; sliceIdx < slices->size(); ++sliceIdx){
 
      const rb::Cluster& slice = (*slices)[sliceIdx];
      if(slice.IsNoise()) continue;
      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;
      }
      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();
      if( flux.fntype==12 && nu.Nu().PdgCode()==12 && nu.CCNC()==0 ) matchSlice=true;
    }
    if( matchSlice ) hasSlice=1;
    else hasSlice=0;

    slEff_Etrue->Fill(nu.Nu().E(),matchSlice?1:0);

    _otree->Fill();
    
    return;
  }//end analyze
  
  void Xbeam::endJob()
  {
  }
  
}//end namespace

namespace ncs{
  
  DEFINE_ART_MODULE(Xbeam);
  
}