NueAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       NueAna
// Module Type: analyzer
// File:        NueAna_module.cc
//
// dorota.stefan@cern.ch, robert.sulej@cern.ch, tjyang@fnal.gov
//
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "art_root_io/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// LArSoft includes
#include "larcore/Geometry/Geometry.h"
#include "lardataobj/RecoBase/Track.h"
#include "lardataobj/RecoBase/Hit.h"
#include "lardataobj/RecoBase/Cluster.h"
#include "lardataobj/RecoBase/Vertex.h"
#include "lardataobj/RecoBase/SpacePoint.h"
#include "lardataobj/RecoBase/TrackHitMeta.h"
#include "lardataobj/RecoBase/Shower.h"
#include "lardataobj/RecoBase/OpFlash.h"
#include "lardata/ArtDataHelper/TrackUtils.h" // lar::util::TrackPitchInView()
#include "lardataobj/AnalysisBase/Calorimetry.h"
#include "lardata/DetectorInfoServices/DetectorPropertiesService.h"
#include "lardata/Utilities/AssociationUtil.h"
#include "larsim/MCCheater/BackTrackerService.h"
#include "larsim/MCCheater/ParticleInventoryService.h"
#include "nusimdata/SimulationBase/MCTruth.h"
#include "larreco/RecoAlg/PMAlg/Utilities.h"
#include "larreco/RecoAlg/TrackMomentumCalculator.h"
#include "larreco/Calorimetry/CalorimetryAlg.h"
#include "larcoreobj/SummaryData/POTSummary.h"
#include "nusimdata/SimulationBase/MCFlux.h"

// ROOT includes
#include "TTree.h"
#include "TTimeStamp.h"

// C/C++ libraries
#include <memory>
#include <utility>

constexpr int kMaxTrack      = 1000;  //maximum number of tracks
constexpr int kMaxShower     = 1000;  //maximum number of showers
constexpr int kMaxHits       = 40000; //maximum number of hits;
constexpr int kMaxVertices   = 100;    //max number of 3D vertices
constexpr int kMaxPrimaries  = 20000;  //maximum number of primary particles
constexpr int kMaxFlash      = 1000;  //maximum number of flashes
//constexpr int kMaxTrackHits  = 1000;  //maximum number of track trajectory points
namespace dunefd {
        class Hit2D;
  class NueAna;
        class IniSegAlg;
}

class dunefd::NueAna : public art::EDAnalyzer {
public:
  explicit NueAna(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  NueAna(NueAna const &) = delete;
  NueAna(NueAna &&) = delete;
  NueAna & operator = (NueAna const &) = delete;
  NueAna & operator = (NueAna &&) = delete;

  // Required functions.
  void analyze(art::Event const & e) override;

  // Selected optional functions.
  void beginJob() override;
  void endJob() override;

  //void beginSubRun(const art::SubRun& sr);
  void endSubRun(const art::SubRun& sr) override;

  void reconfigure(fhicl::ParameterSet const& p) ;

private:

  void ResetVars();
  bool insideFidVol(const double posX, const double posY, const double posZ);

  // Declare member data here.
  TTree *fTree;
  TTree* fPOT;

  // TTree variables

  // Run information
  int run;
  int subrun;
  int event;
  float evttime;
  float taulife;
  short isdata;
  double pot;

  // Track information
  int ntracks_reco;                   //number of reconstructed tracks
  int   trkid[kMaxTrack];             //track id from recob::Track::ID(), this does not have to be the index of track
  float trkstartx[kMaxTrack];         //track start position (cm)
  float trkstarty[kMaxTrack];
  float trkstartz[kMaxTrack];
  float trkendx[kMaxTrack];           //track end position (cm)
  float trkendy[kMaxTrack];
  float trkendz[kMaxTrack];
  float trkstartdcosx[kMaxTrack];     //track start direction cosine
  float trkstartdcosy[kMaxTrack];
  float trkstartdcosz[kMaxTrack];
  float trkenddcosx[kMaxTrack];       //track end direction cosine
  float trkenddcosy[kMaxTrack];
  float trkenddcosz[kMaxTrack];
  float trklen[kMaxTrack];            //track length (cm)
  float trkke[kMaxTrack][3];          //track kinetic energy (in 3 planes)
  float trkpida[kMaxTrack][3];        //track PIDA (in 3 planes)
  int   trkbestplane[kMaxTrack];      //best plane for trkke and trkpida
  float trkmomrange[kMaxTrack];       //Track momentum using range

  //geant information for the track
  int   trkg4id[kMaxTrack];           //geant track id for the track
  int   trkg4pdg[kMaxTrack];          //pdg of geant particle
  float trkg4startx[kMaxTrack];       //start position of geant particle
  float trkg4starty[kMaxTrack];
  float trkg4startz[kMaxTrack];
  float trkg4initdedx[kMaxTrack];     //initial dE/dx of the track using true energy (MeV/cm)

  int nhits;
  int nhits_stored;
  Short_t  hit_plane[kMaxHits];      //plane number
  Short_t  hit_wire[kMaxHits];       //wire number
  Int_t    hit_channel[kMaxHits];    //channel ID
  Short_t  hit_tpc[kMaxHits];        //tpc
  Float_t  hit_peakT[kMaxHits];      //peak time
  Float_t  hit_charge[kMaxHits];     //charge (area)
  Float_t  hit_summedADC[kMaxHits];  //summed ADC
  Float_t  hit_startT[kMaxHits];     //hit start time
  Float_t  hit_endT[kMaxHits];       //hit end time
  Int_t    hit_trkkey[kMaxHits];     //track index if hit is associated with a track
  Float_t  hit_dQds[kMaxHits];       //hit dQ/ds
  Float_t  hit_dEds[kMaxHits];       //hit dE/ds
  Float_t  hit_resrange[kMaxHits];   //hit residual range
  Int_t    hit_shwkey[kMaxHits];     //shower index if hit is associated with a shower

  // vertex information
  int infidvol;
  Short_t  nvtx;                     //number of vertices
  Float_t  vtx[kMaxVertices][3];     //vtx[3] 

  Float_t       vtxrecomc;              // distance between mc and reco vtx
  Float_t       vtxrecomcx;             
  Float_t       vtxrecomcy;             
  Float_t       vtxrecomcz;             

  // shower information
  int nshws;                         //number of showers
  int shwid[kMaxShower];             //recob::Shower::ID()
  Float_t shwdcosx[kMaxShower];      //shower direction cosine
  Float_t shwdcosy[kMaxShower];
  Float_t shwdcosz[kMaxShower];
  Float_t shwstartx[kMaxShower];     //shower start position (cm)
  Float_t shwstarty[kMaxShower];
  Float_t shwstartz[kMaxShower];
  Float_t shwenergy[kMaxShower][3];  //shower energy measured on the 3 planes (GeV)
  Float_t shwdedx[kMaxShower][3];    //shower dE/dx of the initial track measured on the 3 plane (MeV/cm)
  int shwbestplane[kMaxShower];      //recommended plane for energy and dE/dx information
  int   shwg4id[kMaxShower];          //geant track id for the shower

  // flash information
  int    flash_total;                //total number of flashes
  Float_t flash_time[kMaxFlash];     //flash time
  Float_t flash_width[kMaxFlash];    //flash width
  Float_t flash_abstime[kMaxFlash];  //flash absolute time
  Float_t flash_YCenter[kMaxFlash];  //flash y center (cm)
  Float_t flash_YWidth[kMaxFlash];   //flash y width (cm)
  Float_t flash_ZCenter[kMaxFlash];  //flash z center (cm)
  Float_t flash_ZWidth[kMaxFlash];   //flash z width (cm)
  Float_t flash_TotalPE[kMaxFlash];  //flash total pe

  //mctruth information
  Int_t     mcevts_truth;    //number of neutrino Int_teractions in the spill
  Int_t     nuPDG_truth;     //neutrino PDG code
  Int_t     ccnc_truth;      //0=CC 1=NC
  Int_t     mode_truth;      //0=QE/El, 1=RES, 2=DIS, 3=Coherent production
  Float_t  enu_truth;       //true neutrino energy
  Float_t  Q2_truth;        //Momentum transfer squared
  Float_t  W_truth;         //hadronic invariant mass
  Float_t  X_truth;
  Float_t  Y_truth;
  Int_t    hitnuc_truth;    //hit nucleon
  Int_t    target_truth;    //hit nucleus 
  Float_t  nuvtxx_truth;    //neutrino vertex x
  Float_t  nuvtxy_truth;    //neutrino vertex y
  Float_t  nuvtxz_truth;    //neutrino vertex z
  Float_t  nu_dcosx_truth;  //neutrino dcos x
  Float_t  nu_dcosy_truth;  //neutrino dcos y
  Float_t  nu_dcosz_truth;  //neutrino dcos z
  Float_t  lep_mom_truth;   //lepton momentum
  Float_t  lep_dcosx_truth; //lepton dcos x
  Float_t  lep_dcosy_truth; //lepton dcos y
  Float_t  lep_dcosz_truth; //lepton dcos z
  Float_t  t0_truth;        // t0

  // === Storing Geant4 MC Truth Information ===
  int no_primaries;                             //<---Number of primary Geant4 particles in the event
  int geant_list_size;                          //<---Number of Geant4 particles tracked
  int pdg[kMaxPrimaries];                       //<---PDG Code number of this particle
  Float_t Eng[kMaxPrimaries];                   //<---Energy of the particle
  Float_t Px[kMaxPrimaries];                    //<---Px momentum of the particle
  Float_t Py[kMaxPrimaries];                    //<---Py momentum of the particle
  Float_t Pz[kMaxPrimaries];                    //<---Pz momentum of the particle
  Float_t StartPointx[kMaxPrimaries];           //<---X position that this Geant4 particle started at
  Float_t StartPointy[kMaxPrimaries];           //<---Y position that this Geant4 particle started at
  Float_t StartPointz[kMaxPrimaries];           //<---Z position that this Geant4 particle started at
  Float_t EndPointx[kMaxPrimaries];             //<---X position that this Geant4 particle ended at
  Float_t EndPointy[kMaxPrimaries];             //<---Y position that this Geant4 particle ended at
  Float_t EndPointz[kMaxPrimaries];             //<---Z position that this Geant4 particle ended at
  Float_t Startdcosx[kMaxPrimaries];            //<---X direction cosine that Geant4 particle started at
  Float_t Startdcosy[kMaxPrimaries];            //<---Y direction cosine that Geant4 particle started at
  Float_t Startdcosz[kMaxPrimaries];            //<---Z direction cosine that Geant4 particle started at
  int NumberDaughters[kMaxPrimaries];           //<---Number of Daughters this particle has
  int TrackId[kMaxPrimaries];                   //<---Geant4 TrackID number
  int Mother[kMaxPrimaries];                    //<---TrackID of the mother of this particle
  int process_primary[kMaxPrimaries];           //<---Is this particle primary (primary = 1, non-primary = 1)
  std::vector<std::string> G4Process;         //<---The process which created this particle
  std::vector<std::string> G4FinalProcess;    //<---The last process which this particle went under

  //flux information
  Int_t    ptype_flux;        //Parent GEANT code particle ID
  Float_t  pdpx_flux;        //Parent X momentum at decay point (GeV)
  Float_t  pdpy_flux;        //Parent Y momentum at decay point (GeV)
  Float_t  pdpz_flux;        //Parent Z momentum at decay point (GeV)
  Int_t    pntype_flux;      //oscillated neutrino type
  Float_t  vx_flux;          //X position of hadron/muon decay (cm)
  Float_t  vy_flux;          //Y position of hadron/muon decay (cm)
  Float_t  vz_flux;          //Z position of hadron/muon decay (cm)

  //end of ttree variables

  //Module labels to get data products
  std::string fHitsModuleLabel;
  std::string fClusterModuleLabel;
  std::string fTrackModuleLabel;
  std::string fShowerModuleLabel;
  std::string fVertexModuleLabel;
  std::string fGenieGenModuleLabel;
  std::string fPOTModuleLabel;
  std::string fFlashModuleLabel;
  std::string fCalorimetryModuleLabel;

  double fFidVolCut;

  calo::CalorimetryAlg fCalorimetryAlg;

};


dunefd::NueAna::NueAna(fhicl::ParameterSet const & pset)
  : EDAnalyzer(pset)
  , fCalorimetryAlg(pset.get<fhicl::ParameterSet>("CalorimetryAlg"))
{
  reconfigure(pset);
}

//void dunefd::NueAna::beginSubRun(const art::SubRun& sr){
//
//  auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
//  if (potListHandle)
//    pot = potListHandle->totpot;
//  else
//    pot = 0.;
//  if (fPOT) fPOT->Fill();
//}

void dunefd::NueAna::endSubRun(const art::SubRun& sr){

  auto potListHandle = sr.getHandle< sumdata::POTSummary >(fPOTModuleLabel);
  if (potListHandle)
    pot = potListHandle->totpot;
  else
    pot = 0.;
  if (fPOT) fPOT->Fill();
}

void dunefd::NueAna::analyze(art::Event const & evt)
{

  std::cout << " **************************************** analyze ************ " << std::endl;
  // Implementation of required member function here.
  ResetVars();
  art::ServiceHandle<geo::Geometry> geom;
  art::ServiceHandle<cheat::BackTrackerService> bt_serv;
  art::ServiceHandle<cheat::ParticleInventoryService> pi_serv;
  const sim::ParticleList& plist = pi_serv->ParticleList();

  run = evt.run();
  subrun = evt.subRun();
  event = evt.id().event();
  art::Timestamp ts = evt.time();
  TTimeStamp tts(ts.timeHigh(), ts.timeLow());
  evttime = tts.AsDouble();
  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataFor(evt);
  auto const detProp = art::ServiceHandle<detinfo::DetectorPropertiesService const>()->DataFor(evt, clockData);
  taulife = detProp.ElectronLifetime();
  isdata = evt.isRealData();

  // * hits
  std::vector<art::Ptr<recob::Hit> > hitlist;
  auto hitListHandle = evt.getHandle< std::vector<recob::Hit> >(fHitsModuleLabel);
  if (hitListHandle)
    art::fill_ptr_vector(hitlist, hitListHandle);

  // * tracks
  std::vector<art::Ptr<recob::Track> > tracklist;
  auto trackListHandle = evt.getHandle< std::vector<recob::Track> >(fTrackModuleLabel);
  if (trackListHandle)
    art::fill_ptr_vector(tracklist, trackListHandle);

  // * vertices
  std::vector<art::Ptr<recob::Vertex> > vtxlist;
  auto vtxListHandle = evt.getHandle< std::vector<recob::Vertex> >(fVertexModuleLabel);
  if (vtxListHandle)
    art::fill_ptr_vector(vtxlist, vtxListHandle);

  // * showers
  std::vector<art::Ptr<recob::Shower>> shwlist;
  auto shwListHandle = evt.getHandle<std::vector<recob::Shower> >(fShowerModuleLabel);
  if (shwListHandle)
    art::fill_ptr_vector(shwlist, shwListHandle);
  
  // * flashes
  std::vector<art::Ptr<recob::OpFlash> > flashlist;
  auto flashListHandle = evt.getHandle< std::vector<recob::OpFlash> >(fFlashModuleLabel);
  if (flashListHandle)
    art::fill_ptr_vector(flashlist, flashListHandle);

  // * associations
  art::FindManyP<recob::Hit> fmth(trackListHandle, evt, fTrackModuleLabel);
  art::FindManyP<recob::Hit, recob::TrackHitMeta> fmthm(trackListHandle, evt, fTrackModuleLabel);
  art::FindManyP<recob::SpacePoint> fmhs(hitListHandle, evt, fTrackModuleLabel);
  art::FindMany<anab::Calorimetry>  fmcal(trackListHandle, evt, fCalorimetryModuleLabel);

  //hit information
  nhits = hitlist.size();
  nhits_stored = std::min(nhits, kMaxHits);
  for (int i = 0; i < nhits && i < kMaxHits ; ++i){//loop over hits
    hit_channel[i] = hitlist[i]->Channel();
    hit_plane[i]   = hitlist[i]->WireID().Plane;
    hit_wire[i]    = hitlist[i]->WireID().Wire;
    hit_tpc[i]     = hitlist[i]->WireID().TPC;
    hit_peakT[i]   = hitlist[i]->PeakTime();
    hit_charge[i]  = hitlist[i]->Integral();
    hit_summedADC[i] = hitlist[i]->SummedADC();
    hit_startT[i] = hitlist[i]->PeakTimeMinusRMS();
    hit_endT[i] = hitlist[i]->PeakTimePlusRMS();

  }

  //track information
  ntracks_reco=tracklist.size();

  recob::Track::Vector_t larStart;
  recob::Track::Vector_t larEnd;
  trkf::TrackMomentumCalculator trkm;
  for(int i=0; i<std::min(int(tracklist.size()),kMaxTrack);++i){
    recob::Track::Point_t trackStart, trackEnd;
    std::tie(trackStart, trackEnd) = tracklist[i]->Extent(); 
    larStart = tracklist[i]->VertexDirection();
    larEnd = tracklist[i]->EndDirection();

    trkid[i]       = tracklist[i]->ID();
    trkstartx[i]      = trackStart.X();
    trkstarty[i]      = trackStart.Y();
    trkstartz[i]      = trackStart.Z();
    trkendx[i]        = trackEnd.X();
    trkendy[i]        = trackEnd.Y();
    trkendz[i]        = trackEnd.Z();
    trkstartdcosx[i]  = larStart.X();
    trkstartdcosy[i]  = larStart.Y();
    trkstartdcosz[i]  = larStart.Z();
    trkenddcosx[i]    = larEnd.X();
    trkenddcosy[i]    = larEnd.Y();
    trkenddcosz[i]    = larEnd.Z();
    trklen[i]         = tracklist[i]->Length();
    if (!std::isnan(trklen[i])) trkmomrange[i]    = trkm.GetTrackMomentum(trklen[i],13);
    if (fmthm.isValid()){
      auto vhit = fmthm.at(i);
      auto vmeta = fmthm.data(i);
      for (size_t h = 0; h < vhit.size(); ++h){
        if (vhit[h].key()<kMaxHits){
          hit_trkkey[vhit[h].key()] = tracklist[i].key();
          if (vmeta[h]->Dx()){
            hit_dQds[vhit[h].key()] = vhit[h]->Integral()*fCalorimetryAlg.LifetimeCorrection(clockData, detProp, vhit[h]->PeakTime())/vmeta[h]->Dx();
            hit_dEds[vhit[h].key()] = fCalorimetryAlg.dEdx_AREA(clockData, detProp, *vhit[h], vmeta[h]->Dx());
          }
          hit_resrange[vhit[h].key()] = tracklist[i]->Length(vmeta[h]->Index());
        }
      }//loop over all hits
    }//fmthm is valid
    else if (fmth.isValid()){
      std::vector< art::Ptr<recob::Hit> > vhit = fmth.at(i);
      for (size_t h = 0; h < vhit.size(); ++h){
        if (vhit[h].key()<kMaxHits){
          hit_trkkey[vhit[h].key()] = tracklist[i].key();
        }
      }
    }
    if (fmcal.isValid()){
      unsigned maxnumhits = 0;
      std::vector<const anab::Calorimetry*> calos = fmcal.at(i);
      for (auto const& calo : calos){
        if (calo->PlaneID().isValid){
          trkke[i][calo->PlaneID().Plane] = calo->KineticEnergy();
          if (calo->dEdx().size()>maxnumhits){
            maxnumhits = calo->dEdx().size();
            trkbestplane[i] = calo->PlaneID().Plane;
          }
          double pida = 0;
          int used_trkres = 0;
          for (size_t ip = 0; ip<calo->dEdx().size(); ++ip){
            if (calo->ResidualRange()[ip]<30){
              pida += calo->dEdx()[ip]*pow(calo->ResidualRange()[ip],0.42);
              ++used_trkres;
            }
          }
          if (used_trkres) pida/=used_trkres;
          trkpida[i][calo->PlaneID().Plane] = pida;
        }
      }
    }
    if (!isdata&&fmth.isValid()){
      // Find true track for each reconstructed track
      int TrackID = 0;
      std::vector< art::Ptr<recob::Hit> > allHits = fmth.at(i);
      
      std::map<int,double> trkide;
      for(size_t h = 0; h < allHits.size(); ++h){
        art::Ptr<recob::Hit> hit = allHits[h];
        std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
        for(size_t e = 0; e < TrackIDs.size(); ++e){
          trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
        }           
      }
      // Work out which IDE despoited the most charge in the hit if there was more than one.
      double maxe = -1;
      double tote = 0;
      for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
        tote += ii->second;
        if ((ii->second)>maxe){
          maxe = ii->second;
          TrackID = ii->first;
        }
      }
      // Now have trackID, so get PdG code and T0 etc.
      const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
      if (particle){
        trkg4id[i] = TrackID;
        trkg4pdg[i] = particle->PdgCode();
        trkg4startx[i] = particle->Vx();
        trkg4starty[i] = particle->Vy();
        trkg4startz[i] = particle->Vz();
        float sum_energy = 0;
        int numhits = 0;
        //std::map<float,float> hite;
        double x = 0;
        double y = 0;
        double z = 0;
        double mindis = 1e10;
        //find the closest point to the neutrino vertex
        for(size_t h = 0; h < allHits.size(); ++h){
          art::Ptr<recob::Hit> hit = allHits[h];
          if (hit->WireID().Plane==2){
            std::vector<art::Ptr<recob::SpacePoint> > spts = fmhs.at(hit.key());
            if (spts.size()){
              double dis = sqrt(pow(spts[0]->XYZ()[0]-trkg4startx[i],2)+
                                pow(spts[0]->XYZ()[1]-trkg4starty[i],2)+
                                pow(spts[0]->XYZ()[2]-trkg4startz[i],2));
              if (dis<mindis){
                mindis = dis;
                x = spts[0]->XYZ()[0];
                y = spts[0]->XYZ()[1];
                z = spts[0]->XYZ()[2];
              }
            }
          }
        }
        for(size_t h = 0; h < allHits.size(); ++h){
          art::Ptr<recob::Hit> hit = allHits[h];
          if (hit->WireID().Plane==2){
            std::vector<art::Ptr<recob::SpacePoint> > spts = fmhs.at(hit.key());
            if (spts.size()){
              if (sqrt(pow(spts[0]->XYZ()[0]-x,2)+
                       pow(spts[0]->XYZ()[1]-y,2)+
                       pow(spts[0]->XYZ()[2]-z,2))<3){
                std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
                float toten = 0;
                for(size_t e = 0; e < TrackIDs.size(); ++e){
                  //sum_energy += TrackIDs[e].energy;
                  toten+=TrackIDs[e].energy;
                }
                if (toten){
                  sum_energy += toten;
                  ++numhits;
                }
              }
            }
          }
        }

        float pitch = 0;
        float dis1 = sqrt(pow(trkstartx[i]-trkg4startx[i],2)+
                          pow(trkstarty[i]-trkg4starty[i],2)+
                          pow(trkstartz[i]-trkg4startz[i],2));
        float dis2 = sqrt(pow(trkendx[i]-trkg4startx[i],2)+
                          pow(trkendy[i]-trkg4starty[i],2)+
                          pow(trkendz[i]-trkg4startz[i],2));
        if (dis1<dis2){
          try{
            pitch = lar::util::TrackPitchInView(*(tracklist[i]),geo::kZ,0);
          }
          catch(...){
            pitch = 0;
          }
        }
        else{
          try{
            pitch = lar::util::TrackPitchInView(*(tracklist[i]), geo::kZ,tracklist[i]->NumberTrajectoryPoints()-1);
          }
          catch(...){
            pitch = 0;
          }
        }
        if ( pitch && numhits ) {
          trkg4initdedx[i] = sum_energy/(numhits*pitch);
        }
        else{
          trkg4initdedx[i] = 0;
        }
      }//if (particle)
    }//MC
  }

  //vertex information
  nvtx = vtxlist.size();
  for (int i = 0; i < nvtx && i < kMaxVertices ; ++i){//loop over hits
    Double_t xyz[3] = {};
    vtxlist[i]->XYZ(xyz);
    for (size_t j = 0; j<3; ++j) vtx[i][j] = xyz[j];
  }

  //shower information
  if (shwListHandle.isValid()){
  art::FindManyP<recob::Hit> fmsh(shwListHandle, evt, fShowerModuleLabel);

  nshws = shwlist.size();
  
  for (int i = 0; i<std::min(int(shwlist.size()),kMaxShower); ++i){
    shwid[i] = shwlist[i]->ID();
    shwdcosx[i] = shwlist[i]->Direction().X(); 
    shwdcosy[i] = shwlist[i]->Direction().Y(); 
    shwdcosz[i] = shwlist[i]->Direction().Z(); 
    shwstartx[i] = shwlist[i]->ShowerStart().X();
    shwstarty[i] = shwlist[i]->ShowerStart().Y();
    shwstartz[i] = shwlist[i]->ShowerStart().Z();
    for (size_t j = 0; j<(shwlist[i]->Energy()).size(); ++j){
      shwenergy[i][j] = shwlist[i]->Energy()[j];
    }      
    for (size_t j = 0; j<(shwlist[i]->dEdx()).size(); ++j){
      shwdedx[i][j] = shwlist[i]->dEdx()[j];
    }
    shwbestplane[i] = shwlist[i]->best_plane();
    if (fmsh.isValid()){
      auto vhit = fmsh.at(i);
      for (size_t h = 0; h < vhit.size(); ++h){
        if (vhit[h].key()<kMaxHits){
          hit_shwkey[vhit[h].key()] = shwlist[i].key();
        }
      }
    }
    if (!isdata&&fmsh.isValid()){
      // Find true track for each reconstructed track
      int TrackID = 0;
      std::vector< art::Ptr<recob::Hit> > allHits = fmsh.at(i);
      std::map<int,double> trkide;
      for(size_t h = 0; h < allHits.size(); ++h){
        art::Ptr<recob::Hit> hit = allHits[h];
        std::vector<sim::TrackIDE> TrackIDs = bt_serv->HitToTrackIDEs(clockData, hit);
        for(size_t e = 0; e < TrackIDs.size(); ++e){
          trkide[TrackIDs[e].trackID] += TrackIDs[e].energy;
        }           
      }
      // Work out which IDE despoited the most charge in the hit if there was more than one.
      double maxe = -1;
      double tote = 0;
      for (std::map<int,double>::iterator ii = trkide.begin(); ii!=trkide.end(); ++ii){
        tote += ii->second;
        if ((ii->second)>maxe){
          maxe = ii->second;
          TrackID = ii->first;
        }
      }
      // Now have trackID, so get PdG code and T0 etc.
      const simb::MCParticle *particle = pi_serv->TrackIdToParticle_P(TrackID);
      if (particle){
        shwg4id[i] = TrackID;
      }
    }
  }
  }
  // flash information
  flash_total = flashlist.size();
  for ( int f = 0; f < std::min(flash_total,kMaxHits); ++f ) {
    flash_time[f]      = flashlist[f]->Time();
    flash_width[f]     = flashlist[f]->TimeWidth();
    flash_abstime[f]   = flashlist[f]->AbsTime();
    flash_YCenter[f]   = flashlist[f]->YCenter();
    flash_YWidth[f]    = flashlist[f]->YWidth();
    flash_ZCenter[f]   = flashlist[f]->ZCenter();
    flash_ZWidth[f]    = flashlist[f]->ZWidth();
    flash_TotalPE[f]   = flashlist[f]->TotalPE();
  }



  if (!isdata){

    // * MC truth information
    std::vector<art::Ptr<simb::MCTruth> > mclist;
    auto mctruthListHandle = evt.getHandle< std::vector<simb::MCTruth> >(fGenieGenModuleLabel);
    if (mctruthListHandle)
      art::fill_ptr_vector(mclist, mctruthListHandle);

    std::vector<art::Ptr<simb::MCFlux> > fluxlist;
    auto mcfluxListHandle = evt.getHandle< std::vector<simb::MCFlux> >(fGenieGenModuleLabel);
    if (mcfluxListHandle)
      art::fill_ptr_vector(fluxlist, mcfluxListHandle);
    

    mcevts_truth=mclist.size();
    if (mcevts_truth){
      art::Ptr<simb::MCTruth> mctruth = mclist[0];
      if (mctruth->Origin() == simb::kBeamNeutrino){
        nuPDG_truth  = mctruth->GetNeutrino().Nu().PdgCode();
        ccnc_truth   = mctruth->GetNeutrino().CCNC();
        mode_truth   = mctruth->GetNeutrino().Mode();
        Q2_truth     = mctruth->GetNeutrino().QSqr();
        W_truth      = mctruth->GetNeutrino().W();
        X_truth      = mctruth->GetNeutrino().X();
        Y_truth      = mctruth->GetNeutrino().Y();
        hitnuc_truth = mctruth->GetNeutrino().HitNuc();
        target_truth = mctruth->GetNeutrino().Target();
        enu_truth    = mctruth->GetNeutrino().Nu().E();
        nuvtxx_truth = mctruth->GetNeutrino().Nu().Vx();
        nuvtxy_truth = mctruth->GetNeutrino().Nu().Vy();
        nuvtxz_truth = mctruth->GetNeutrino().Nu().Vz();
        if (mctruth->GetNeutrino().Nu().P()){
          nu_dcosx_truth = mctruth->GetNeutrino().Nu().Px()/mctruth->GetNeutrino().Nu().P();
          nu_dcosy_truth = mctruth->GetNeutrino().Nu().Py()/mctruth->GetNeutrino().Nu().P();
          nu_dcosz_truth = mctruth->GetNeutrino().Nu().Pz()/mctruth->GetNeutrino().Nu().P();
        }
        lep_mom_truth = mctruth->GetNeutrino().Lepton().P();
        if (mctruth->GetNeutrino().Lepton().P()){
          lep_dcosx_truth = mctruth->GetNeutrino().Lepton().Px()/mctruth->GetNeutrino().Lepton().P();
          lep_dcosy_truth = mctruth->GetNeutrino().Lepton().Py()/mctruth->GetNeutrino().Lepton().P();
          lep_dcosz_truth = mctruth->GetNeutrino().Lepton().Pz()/mctruth->GetNeutrino().Lepton().P();
        }
        
        if (mctruth->NParticles()){
          simb::MCParticle particle = mctruth->GetParticle(0);
          t0_truth = particle.T();
        }
        
        
        float mindist2 = 9999; // cm;
        TVector3 nuvtx(nuvtxx_truth, nuvtxy_truth, nuvtxz_truth);
        infidvol = insideFidVol(nuvtxx_truth, nuvtxy_truth, nuvtxz_truth); 
        //find the closest reco vertex to the neutrino mc truth
        if (infidvol)
          {
            // vertex is when at least two tracks meet
            for(size_t i = 0; i < vtxlist.size(); ++i){ // loop over vertices
              Double_t xyz[3] = {};
              vtxlist[i]->XYZ(xyz);
              TVector3 vtxreco(xyz);
              float dist2 = pma::Dist2(vtxreco, nuvtx);
              if (dist2 < mindist2)
                {
                  mindist2 = dist2;
                  vtxrecomc = std::sqrt(dist2);
                  vtxrecomcx = vtxreco.X() - nuvtxx_truth;
                  vtxrecomcy = vtxreco.Y() - nuvtxy_truth;
                  vtxrecomcz = vtxreco.Z() - nuvtxz_truth;
                }
            }
            
            // two endpoints of tracks are somehow also vertices...
            for (size_t i = 0; i < tracklist.size(); ++i){ // loop over tracks
              float dist2 = pma::Dist2(tracklist[i]->Vertex(), nuvtx);
              if (dist2 < mindist2)
                {
                  mindist2 = dist2;
                  vtxrecomc = std::sqrt(dist2);
                  vtxrecomcx = tracklist[i]->Vertex().X() - nuvtxx_truth;
                  vtxrecomcy = tracklist[i]->Vertex().Y() - nuvtxy_truth;
                  vtxrecomcz = tracklist[i]->Vertex().Z() - nuvtxz_truth;
                  
                }
              dist2 = pma::Dist2(tracklist[i]->End(), nuvtx);
              if (dist2 < mindist2)
                {
                  mindist2 = dist2;
                  vtxrecomc = std::sqrt(dist2);
                  vtxrecomcx = tracklist[i]->End().X() - nuvtxx_truth;
                  vtxrecomcy = tracklist[i]->End().Y() - nuvtxy_truth;
                  vtxrecomcz = tracklist[i]->End().Z() - nuvtxz_truth;
                  
                }
            }
          }
      }//is neutrino
    }
    
    if (fluxlist.size()){
      ptype_flux  = fluxlist[0]->fptype;
      pdpx_flux   = fluxlist[0]->fpdpx;
      pdpy_flux   = fluxlist[0]->fpdpy;
      pdpz_flux   = fluxlist[0]->fpdpz;
      pntype_flux = fluxlist[0]->fntype;
      vx_flux     = fluxlist[0]->fvx;
      vy_flux     = fluxlist[0]->fvy;
      vz_flux     = fluxlist[0]->fvz;
    }

    //save g4 particle information
    std::vector<const simb::MCParticle* > geant_part;
    
    // ### Looping over all the Geant4 particles from the BackTrackerService ###
    for(size_t p = 0; p < plist.size(); ++p) 
      {
        // ### Filling the vector with MC Particles ###
        geant_part.push_back(plist.Particle(p)); 
      }
    
    //std::cout<<"No of geant part= "<<geant_part.size()<<std::endl;
    
    // ### Setting a string for primary ###
    std::string pri("primary");
    
    int primary=0;
    int geant_particle=0;
    
    // ############################################################
    // ### Determine the number of primary particles from geant ###
    // ############################################################
    for( unsigned int i = 0; i < geant_part.size(); ++i ){
      geant_particle++;
      // ### Counting the number of primary particles ###
      if(geant_part[i]->Process()==pri)
        { primary++;}
    }//<---End i loop
    
    
    // ### Saving the number of primary particles ###
    no_primaries=primary;
    // ### Saving the number of Geant4 particles ###
    geant_list_size=geant_particle;
    
    // ### Looping over all the Geant4 particles ###
    for( unsigned int i = 0; i < geant_part.size(); ++i ){
   
      // ### If this particle is primary, set = 1 ###
      if(geant_part[i]->Process()==pri)
        {process_primary[i]=1;}
      // ### If this particle is not-primary, set = 0 ###
      else
        {process_primary[i]=0;}
      
      // ### Saving the particles mother TrackID ###
      Mother[i]=geant_part[i]->Mother();
      // ### Saving the particles TrackID ###
      TrackId[i]=geant_part[i]->TrackId();
      // ### Saving the PDG Code ###
      pdg[i]=geant_part[i]->PdgCode();
      // ### Saving the particles Energy ###
      Eng[i]=geant_part[i]->E();
      
      // ### Saving the Px, Py, Pz info ###
      Px[i]=geant_part[i]->Px();
      Py[i]=geant_part[i]->Py();
      Pz[i]=geant_part[i]->Pz();
      
      // ### Saving the Start and End Point for this particle ###
      StartPointx[i]=geant_part[i]->Vx();
      StartPointy[i]=geant_part[i]->Vy();
      StartPointz[i]=geant_part[i]->Vz();
      EndPointx[i]=geant_part[i]->EndPosition()[0];
      EndPointy[i]=geant_part[i]->EndPosition()[1];
      EndPointz[i]=geant_part[i]->EndPosition()[2];

      // ### Saving the processes for this particle ###
      //std::cout<<"finding proc"<<std::endl;
      G4Process.push_back( geant_part[i]->Process() );
      G4FinalProcess.push_back( geant_part[i]->EndProcess() );
      //std::cout<<"found proc"<<std::endl;
//      std::cout << "ID " << TrackId[i] << ", pdg " << pdg[i] << ", Start X,Y,Z " << StartPointx[i] << ", " << StartPointy[i] << ", " << StartPointz[i]
//              << ", End XYZ " << EndPointx[i] << ", " << EndPointy[i] << ", " << EndPointz[i] << ", Start Proc " << G4Process[i] << ", End Proc " << G4FinalProcess[i]
//              << std::endl;

      // ### Saving the Start direction cosines for this particle ###
      Startdcosx[i] = geant_part[i]->Momentum(0).Px() / geant_part[i]->Momentum(0).P();
      Startdcosy[i] = geant_part[i]->Momentum(0).Py() / geant_part[i]->Momentum(0).P();
      Startdcosz[i] = geant_part[i]->Momentum(0).Pz() / geant_part[i]->Momentum(0).P();
      // ### Saving the number of Daughters for this particle ###
      NumberDaughters[i]=geant_part[i]->NumberDaughters();
      
    } //geant particles


  }//is neutrino
  fTree->Fill();
}

void dunefd::NueAna::beginJob()
{
  // Implementation of optional member function here.
  art::ServiceHandle<art::TFileService> tfs;
  fTree = tfs->make<TTree>("nueana","analysis tree");
  fTree->Branch("run",&run,"run/I");
  fTree->Branch("subrun",&subrun,"subrun/I");
  fTree->Branch("event",&event,"event/I");
  fTree->Branch("evttime",&evttime,"evttime/F");
  fTree->Branch("taulife",&taulife,"taulife/F");
  fTree->Branch("isdata",&isdata,"isdata/S");
  fTree->Branch("ntracks_reco",&ntracks_reco,"ntracks_reco/I");
  fTree->Branch("trkid",trkid,"trkid[ntracks_reco]/I");
  fTree->Branch("trkstartx",trkstartx,"trkstartx[ntracks_reco]/F");
  fTree->Branch("trkstarty",trkstarty,"trkstarty[ntracks_reco]/F");
  fTree->Branch("trkstartz",trkstartz,"trkstartz[ntracks_reco]/F");
  fTree->Branch("trkendx",trkendx,"trkendx[ntracks_reco]/F");
  fTree->Branch("trkendy",trkendy,"trkendy[ntracks_reco]/F");
  fTree->Branch("trkendz",trkendz,"trkendz[ntracks_reco]/F");
  fTree->Branch("trkstartdcosx",trkstartdcosx,"trkstartdcosx[ntracks_reco]/F");
  fTree->Branch("trkstartdcosy",trkstartdcosy,"trkstartdcosy[ntracks_reco]/F");
  fTree->Branch("trkstartdcosz",trkstartdcosz,"trkstartdcosz[ntracks_reco]/F");
  fTree->Branch("trkenddcosx",trkenddcosx,"trkenddcosx[ntracks_reco]/F");
  fTree->Branch("trkenddcosy",trkenddcosy,"trkenddcosy[ntracks_reco]/F");
  fTree->Branch("trkenddcosz",trkenddcosz,"trkenddcosz[ntracks_reco]/F");
  fTree->Branch("trklen",trklen,"trklen[ntracks_reco]/F");
  fTree->Branch("trkbestplane",trkbestplane,"trkbestplane[ntracks_reco]/I");
  fTree->Branch("trkmomrange",trkmomrange,"trkmomrange[ntracks_reco]/F");
  fTree->Branch("trkke",trkke,"trkke[ntracks_reco][3]/F");
  fTree->Branch("trkpida",trkpida,"trkpida[ntracks_reco][3]/F");
  fTree->Branch("trkg4id",trkg4id,"trkg4id[ntracks_reco]/I");
  fTree->Branch("trkg4pdg",trkg4pdg,"trkg4pdg[ntracks_reco]/I");
  fTree->Branch("trkg4startx",trkg4startx,"trkg4startx[ntracks_reco]/F");
  fTree->Branch("trkg4starty",trkg4starty,"trkg4starty[ntracks_reco]/F");
  fTree->Branch("trkg4startz",trkg4startz,"trkg4startz[ntracks_reco]/F");
  fTree->Branch("trkg4initdedx",trkg4initdedx,"trkg4initdedx[ntracks_reco]/F");
  fTree->Branch("nshws",&nshws,"nshws/I");
  fTree->Branch("shwid",shwid,"shwid[nshws]/I");
  fTree->Branch("shwdcosx",shwdcosx,"shwdcosx[nshws]/F");
  fTree->Branch("shwdcosy",shwdcosy,"shwdcosy[nshws]/F");
  fTree->Branch("shwdcosz",shwdcosz,"shwdcosz[nshws]/F");
  fTree->Branch("shwstartx",shwstartx,"shwstartx[nshws]/F");
  fTree->Branch("shwstarty",shwstarty,"shwstarty[nshws]/F");
  fTree->Branch("shwstartz",shwstartz,"shwstartz[nshws]/F");
  fTree->Branch("shwenergy",shwenergy,"shwenergy[nshws][3]/F");
  fTree->Branch("shwdedx",shwdedx,"shwdedx[nshws][3]/F");
  fTree->Branch("shwbestplane",shwbestplane,"shwbestplane[nshws]/I");
  fTree->Branch("shwg4id",shwg4id,"shwg4id[nshws]/I");
  fTree->Branch("flash_total"  ,&flash_total ,"flash_total/I");
  fTree->Branch("flash_time"   ,flash_time   ,"flash_time[flash_total]/F");
  fTree->Branch("flash_width"  ,flash_width  ,"flash_width[flash_total]/F");
  fTree->Branch("flash_abstime",flash_abstime,"flash_abstime[flash_total]/F");
  fTree->Branch("flash_YCenter",flash_YCenter,"flash_YCenter[flash_total]/F");
  fTree->Branch("flash_YWidth" ,flash_YWidth ,"flash_YWidth[flash_total]/F");
  fTree->Branch("flash_ZCenter",flash_ZCenter,"flash_ZCenter[flash_total]/F");
  fTree->Branch("flash_ZWidth" ,flash_ZWidth ,"flash_ZWidth[flash_total]/F");
  fTree->Branch("flash_TotalPE",flash_TotalPE,"flash_TotalPE[flash_total]/F");
  fTree->Branch("nhits",&nhits,"nhits/I");
  fTree->Branch("nhits_stored",&nhits_stored,"nhits_stored/I");
  fTree->Branch("hit_plane",hit_plane,"hit_plane[nhits_stored]/S");
  fTree->Branch("hit_tpc",hit_tpc,"hit_tpc[nhits_stored]/S");
  fTree->Branch("hit_wire",hit_wire,"hit_wire[nhits_stored]/S");
  fTree->Branch("hit_channel",hit_channel,"hit_channel[nhits_stored]/I");
  fTree->Branch("hit_peakT",hit_peakT,"hit_peakT[nhits_stored]/F");
  fTree->Branch("hit_charge",hit_charge,"hit_charge[nhits_stored]/F");
  fTree->Branch("hit_summedADC",hit_summedADC,"hit_summedADC[nhits_stored]/F");
  fTree->Branch("hit_startT",hit_startT,"hit_startT[nhits_stored]/F");
  fTree->Branch("hit_endT",hit_endT,"hit_endT[nhits_stored]/F");
  fTree->Branch("hit_trkkey",hit_trkkey,"hit_trkkey[nhits_stored]/I");
  fTree->Branch("hit_dQds",hit_dQds,"hit_dQds[nhits_stored]/F");
  fTree->Branch("hit_dEds",hit_dEds,"hit_dEds[nhits_stored]/F");
  fTree->Branch("hit_resrange",hit_resrange,"hit_resrange[nhits_stored]/F");
  fTree->Branch("hit_shwkey",hit_shwkey,"hit_shwkey[nhits_stored]/I");
  fTree->Branch("infidvol",&infidvol,"infidvol/I");
  fTree->Branch("nvtx",&nvtx,"nvtx/S");
  fTree->Branch("vtx",vtx,"vtx[nvtx][3]/F");
  fTree->Branch("vtxrecomc",&vtxrecomc,"vtxrecomc/F");
  fTree->Branch("vtxrecomcx",&vtxrecomcx,"vtxrecomcx/F");
  fTree->Branch("vtxrecomcy",&vtxrecomcy,"vtxrecomcy/F");
  fTree->Branch("vtxrecomcz",&vtxrecomcz,"vtxrecomcz/F");
  fTree->Branch("mcevts_truth",&mcevts_truth,"mcevts_truth/I");
  fTree->Branch("nuPDG_truth",&nuPDG_truth,"nuPDG_truth/I");
  fTree->Branch("ccnc_truth",&ccnc_truth,"ccnc_truth/I");
  fTree->Branch("mode_truth",&mode_truth,"mode_truth/I");
  fTree->Branch("enu_truth",&enu_truth,"enu_truth/F");
  fTree->Branch("Q2_truth",&Q2_truth,"Q2_truth/F");
  fTree->Branch("W_truth",&W_truth,"W_truth/F");
  fTree->Branch("X_truth",&X_truth,"X_truth/F");
  fTree->Branch("Y_truth",&Y_truth,"Y_truth/F");
  fTree->Branch("hitnuc_truth",&hitnuc_truth,"hitnuc_truth/I");
  fTree->Branch("target_truth",&target_truth,"target_truth/I");
  fTree->Branch("nuvtxx_truth",&nuvtxx_truth,"nuvtxx_truth/F");
  fTree->Branch("nuvtxy_truth",&nuvtxy_truth,"nuvtxy_truth/F");
  fTree->Branch("nuvtxz_truth",&nuvtxz_truth,"nuvtxz_truth/F");
  fTree->Branch("nu_dcosx_truth",&nu_dcosx_truth,"nu_dcosx_truth/F");
  fTree->Branch("nu_dcosy_truth",&nu_dcosy_truth,"nu_dcosy_truth/F");
  fTree->Branch("nu_dcosz_truth",&nu_dcosz_truth,"nu_dcosz_truth/F");
  fTree->Branch("lep_mom_truth",&lep_mom_truth,"lep_mom_truth/F");
  fTree->Branch("lep_dcosx_truth",&lep_dcosx_truth,"lep_dcosx_truth/F");
  fTree->Branch("lep_dcosy_truth",&lep_dcosy_truth,"lep_dcosy_truth/F");
  fTree->Branch("lep_dcosz_truth",&lep_dcosz_truth,"lep_dcosz_truth/F");
  fTree->Branch("t0_truth",&t0_truth,"t0_truth/F");
  fTree->Branch("no_primaries",&no_primaries,"no_primaries/I");
  fTree->Branch("geant_list_size",&geant_list_size,"geant_list_size/I");
  fTree->Branch("pdg",pdg,"pdg[geant_list_size]/I");
  fTree->Branch("Eng",Eng,"Eng[geant_list_size]/F");
  fTree->Branch("Px",Px,"Px[geant_list_size]/F");
  fTree->Branch("Py",Py,"Py[geant_list_size]/F");
  fTree->Branch("Pz",Pz,"Pz[geant_list_size]/F");
  fTree->Branch("StartPointx",StartPointx,"StartPointx[geant_list_size]/F");
  fTree->Branch("StartPointy",StartPointy,"StartPointy[geant_list_size]/F");
  fTree->Branch("StartPointz",StartPointz,"StartPointz[geant_list_size]/F");
  fTree->Branch("EndPointx",EndPointx,"EndPointx[geant_list_size]/F");
  fTree->Branch("EndPointy",EndPointy,"EndPointy[geant_list_size]/F");
  fTree->Branch("EndPointz",EndPointz,"EndPointz[geant_list_size]/F");
  fTree->Branch("Startdcosx",Startdcosx,"Startdcosx[geant_list_size]/F");
  fTree->Branch("Startdcosy",Startdcosy,"Startdcosy[geant_list_size]/F");
  fTree->Branch("Startdcosz",Startdcosz,"Startdcosz[geant_list_size]/F");
  fTree->Branch("NumberDaughters",NumberDaughters,"NumberDaughters[geant_list_size]/I");
  fTree->Branch("Mother",Mother,"Mother[geant_list_size]/I");
  fTree->Branch("TrackId",TrackId,"TrackId[geant_list_size]/I");
  fTree->Branch("process_primary",process_primary,"process_primary[geant_list_size]/I");
  fTree->Branch("G4Process",&G4Process);//,"G4Process[geant_list_size]");
  fTree->Branch("G4FinalProcess",&G4FinalProcess);//,"G4FinalProcess[geant_list_size]");
  fTree->Branch("ptype_flux",&ptype_flux,"ptype_flux/I");
  fTree->Branch("pdpx_flux",&pdpx_flux,"pdpx_flux/F");
  fTree->Branch("pdpy_flux",&pdpy_flux,"pdpy_flux/F");
  fTree->Branch("pdpz_flux",&pdpz_flux,"pdpz_flux/F");
  fTree->Branch("pntype_flux",&pntype_flux,"pntype_flux/I");
  fTree->Branch("vx_flux",&vx_flux,"vx_flux/F");
  fTree->Branch("vy_flux",&vy_flux,"vy_flux/F");
  fTree->Branch("vz_flux",&vz_flux,"vz_flux/F");

  fPOT = tfs->make<TTree>("pottree","pot tree");
  fPOT->Branch("pot",&pot,"pot/D");
  fPOT->Branch("run",&run,"run/I");
  fPOT->Branch("subrun",&subrun,"subrun/I");
}

void dunefd::NueAna::ResetVars(){

  G4Process.clear();
  G4FinalProcess.clear();

  run = -9999;
  subrun = -9999;
  event = -9999;
  evttime = -9999;
  taulife = 0;
  isdata = -9999;

  ntracks_reco = 0;
  for (int i = 0; i < kMaxTrack; ++i){
    trkid[i] = -9999;
    trkstartx[i] = -9999;
    trkstarty[i] = -9999;
    trkstartz[i] = -9999;
    trkendx[i] = -9999;
    trkendy[i] = -9999;
    trkendz[i] = -9999;
    trkstartdcosx[i] = -9999;
    trkstartdcosy[i] = -9999;
    trkstartdcosz[i] = -9999;
    trkenddcosx[i] = -9999;
    trkenddcosy[i] = -9999;
    trkenddcosz[i] = -9999;
    trklen[i] = -9999;
    trkbestplane[i] = -9999;
    trkmomrange[i] = -9999;
    for (int j = 0; j<3; ++j){
      trkke[i][j] = -9999;
      trkpida[i][j] = -9999;
    }
    trkg4id[i] = -9999;
    trkg4pdg[i] = -9999;
    trkg4startx[i] = -9999;
    trkg4starty[i] = -9999;
    trkg4startz[i] = -9999;
    trkg4initdedx[i] = -9999;
  }

  nshws = 0;
  for (int i = 0; i<kMaxShower; ++i){
    shwid[i] = -9999;
    shwdcosx[i] = -9999;
    shwdcosy[i] = -9999;
    shwdcosz[i] = -9999;
    shwstartx[i] = -9999;
    shwstarty[i] = -9999;
    shwstartz[i] = -9999;
    for (int j = 0; j<3; ++j){
      shwenergy[i][j] = -9999;
      shwdedx[i][j] = -9999;
    }
    shwbestplane[i] = -9999;
    trkg4id[i] = -9999;
  }

  flash_total = 0;
  for (int f = 0; f < kMaxFlash; ++f) {
    flash_time[f]    = -9999;
    flash_width[f]   = -9999;
    flash_abstime[f] = -9999;
    flash_YCenter[f] = -9999;
    flash_YWidth[f]  = -9999;
    flash_ZCenter[f] = -9999;
    flash_ZWidth[f]  = -9999;
    flash_TotalPE[f] = -9999;
  }

  nhits = 0;
  nhits_stored = 0;
  for (int i = 0; i<kMaxHits; ++i){
    hit_plane[i] = -9999;
    hit_wire[i] = -9999;
    hit_tpc[i] = -9999;
    hit_channel[i] = -9999;
    hit_peakT[i] = -9999;
    hit_charge[i] = -9999;
    hit_summedADC[i] = -9999;
    hit_startT[i] = -9999;
    hit_endT[i] = -9999;
    hit_trkkey[i] = -9999;
    hit_dQds[i] = -9999;
    hit_dEds[i] = -9999;
    hit_resrange[i] = -9999;
    hit_shwkey[i] = -9999;
  }

  infidvol = 0;
  nvtx = 0;
  for (int i = 0; i<kMaxVertices; ++i){
    vtx[i][0] = -9999;
    vtx[i][1] = -9999;
    vtx[i][2] = -9999;
  }
  vtxrecomc = 9999;
  vtxrecomcx = 9999;
  vtxrecomcy = 9999;
  vtxrecomcz = 9999;

  mcevts_truth = -9999; 
  nuPDG_truth = -9999;  
  ccnc_truth = -9999;    
  mode_truth = -9999;   
  enu_truth = -9999;     
  Q2_truth = -9999;       
  W_truth = -9999;        
  X_truth = -9999;
  Y_truth = -9999;
  hitnuc_truth = -9999;
  target_truth = -9999;
  nuvtxx_truth = -9999;   
  nuvtxy_truth = -9999;   
  nuvtxz_truth = -9999;   
  nu_dcosx_truth = -9999; 
  nu_dcosy_truth = -9999; 
  nu_dcosz_truth = -9999; 
  lep_mom_truth = -9999;  
  lep_dcosx_truth = -9999;
  lep_dcosy_truth = -9999;
  lep_dcosz_truth = -9999;
  t0_truth = -9999;

  no_primaries = -99999;
  geant_list_size=-9999;
  for (int i = 0; i<kMaxPrimaries; ++i){
    pdg[i] = -99999;
    Eng[i] = -99999;
    Px[i] = -99999;
    Py[i] = -99999;
    Pz[i] = -99999;
    StartPointx[i] = -99999;
    StartPointy[i] = -99999;
    StartPointz[i] = -99999;
    EndPointx[i] = -99999;
    EndPointy[i] = -99999;
    EndPointz[i] = -99999;
    Startdcosx[i]= -99999;
    Startdcosy[i]= -99999;
    Startdcosz[i]= -99999;
    NumberDaughters[i] = -99999;
    Mother[i] = -99999;
    TrackId[i] = -99999;
    process_primary[i] = -99999;
  }

  ptype_flux = -99999;
  pdpx_flux = -99999;
  pdpy_flux = -99999;
  pdpz_flux = -99999;
  pntype_flux = -99999;
  vx_flux = -99999;
  vy_flux = -99999;
  vz_flux = -99999;
}

void dunefd::NueAna::endJob()
{
  // Implementation of optional member function here.
}

void dunefd::NueAna::reconfigure(fhicl::ParameterSet const & p)
{
  fHitsModuleLabel     =   p.get< std::string >("HitsModuleLabel");
  fTrackModuleLabel    =   p.get< std::string >("TrackModuleLabel");
  fShowerModuleLabel    =   p.get< std::string >("ShowerModuleLabel");
  fClusterModuleLabel  =   p.get< std::string >("ClusterModuleLabel");
  fVertexModuleLabel   =   p.get< std::string >("VertexModuleLabel");
  fGenieGenModuleLabel =   p.get< std::string >("GenieGenModuleLabel");
  fFidVolCut           =   p.get< double >("FidVolCut");
  fPOTModuleLabel      =   p.get< std::string >("POTModuleLabel"); 
  fFlashModuleLabel    =   p.get< std::string >("FlashModuleLabel");
  fCalorimetryModuleLabel = p.get< std::string >("CalorimetryModuleLabel");
  return;
}

/***********************************************************************/

bool dunefd::NueAna::insideFidVol(const double posX, const double posY, const double posZ) 
{
        
        art::ServiceHandle<geo::Geometry> geom;
        double vtx[3] = {posX, posY, posZ};
        bool inside = false;

        geo::TPCID idtpc = geom->FindTPCAtPosition(vtx);

        if (geom->HasTPC(idtpc))
        {               
                const geo::TPCGeo& tpcgeo = geom->GetElement(idtpc);
                double minx = tpcgeo.MinX(); double maxx = tpcgeo.MaxX();
                double miny = tpcgeo.MinY(); double maxy = tpcgeo.MaxY();
                double minz = tpcgeo.MinZ(); double maxz = tpcgeo.MaxZ();

                for (size_t c = 0; c < geom->Ncryostats(); c++)
                {
                        const geo::CryostatGeo& cryostat = geom->Cryostat(c);
                        for (size_t t = 0; t < cryostat.NTPC(); t++)
                        {       
                                const geo::TPCGeo& tpcg = cryostat.TPC(t);
                                if (tpcg.MinX() < minx) minx = tpcg.MinX();
                                if (tpcg.MaxX() > maxx) maxx = tpcg.MaxX(); 
                                if (tpcg.MinY() < miny) miny = tpcg.MinY();
                                if (tpcg.MaxY() > maxy) maxy = tpcg.MaxY();
                                if (tpcg.MinZ() < minz) minz = tpcg.MinZ();
                                if (tpcg.MaxZ() > maxz) maxz = tpcg.MaxZ();
                        }
                }       

                
                //x
                double dista = fabs(minx - posX);
                double distb = fabs(posX - maxx); 
                if ((posX > minx) && (posX < maxx) &&
                        (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
                //y
                dista = fabs(maxy - posY);
                distb = fabs(posY - miny);
                if (inside && (posY > miny) && (posY < maxy) &&
                        (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
                else inside = false;

                //z
                dista = fabs(maxz - posZ);
                distb = fabs(posZ - minz);
                if (inside && (posZ > minz) && (posZ < maxz) &&
                        (dista > fFidVolCut) && (distb > fFidVolCut)) inside = true;
                else inside = false;
        }
                
        return inside;
}

/***********************************************************************/


DEFINE_ART_MODULE(dunefd::NueAna)