FillSandbox_module.cc

Go to the documentation of this file.

/////////////////////////////////////////////////////////
/// \brief  Sandbox to make objects for numu analysis
///         in CAF
/// \author Kirk Bays
////////////////////////////////////////////////////////
#include <string>

#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "fhiclcpp/ParameterSet.h"

#include "NovaDAQConventions/DAQConventions.h"
#include "NumuSandbox/NumuSandFxs.h"
#include "NumuSandbox/NumuSandObj.h"
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "GeometryObjects/CellGeo.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/Track.h"
#include "RecoBase/Vertex.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/RecoHit.h"
#include "ReMId/ReMId.h"
#include "Utilities/AssociationUtil.h"
#include "Utilities/func/MathUtil.h"
#include "NumuEnergy/NumuE.h" // New

#include "MCCheater/BackTracker.h"
#include "TFile.h"
#include "TMVA/Reader.h"

namespace numusand {

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

    void produce(art::Event & evt);
    void beginRun(art::Run& run);

  private:

    numusand::NumuSandFxs*     fNumuSandFxs;
    std::string                fSlicerLabel;
    std::string                fTrackLabel;
    std::string                fEnergyLabel;
    std::string                fPIDLabel;
    std::string                fVtxLabel;

  };
}


namespace numusand
{

  //----------------------------------------------------------------------
  FillSandbox::FillSandbox(fhicl::ParameterSet const & pset):
    fSlicerLabel (pset.get< std::string >("SlicerLabel")),
    fTrackLabel  (pset.get< std::string >("TrackLabel")),
    fEnergyLabel (pset.get< std::string >("EnergyLabel")),
    fPIDLabel    (pset.get< std::string >("PIDLabel")),
    fVtxLabel    (pset.get< std::string >("VtxLabel"))
  {
    produces< std::vector<numusand::NumuSandObj> >();
    produces< art::Assns<numusand::NumuSandObj, rb::Cluster> >();
  }
  

  //----------------------------------------------------------------------
  FillSandbox::~FillSandbox()
  {
  }

  //-------------------------------------------------------------------
  void FillSandbox::beginRun(art::Run& run)
  {
    fNumuSandFxs = new NumuSandFxs();
  } //end beginRun

  //-------------------------------------------------------------------
  void FillSandbox::produce(art::Event& evt)
  {

    art::Handle<std::vector<rb::Cluster> > slicevec;
    evt.getByLabel(fSlicerLabel,slicevec);

    if(slicevec->empty()) {
      mf::LogWarning ("No Slices")<<"No Slices in the input file";
      return;
    }

    art::ServiceHandle<cheat::BackTracker> bt;
    art::ServiceHandle<geo::Geometry>  geom;
    art::ServiceHandle<geo::LiveGeometry>  livegeom;

    art::FindManyP<rb::Track> trackAssnList(slicevec, evt, fTrackLabel);

    // set up associations between slices and energy
    art::FindManyP<numue::NumuE> eAss(slicevec,evt,fEnergyLabel);
    
    std::unique_ptr< std::vector<numusand::NumuSandObj> > outputObjects(new std::vector<numusand::NumuSandObj>);
    std::unique_ptr< art::Assns<numusand::NumuSandObj, rb::Cluster> > assoc(new art::Assns<numusand::NumuSandObj, rb::Cluster>);

    for(size_t i = 0; i < slicevec->size(); ++i){

      art::Ptr<rb::Cluster> slice(slicevec,i);
      
      //get the total no. of cell hits
      art::PtrVector<rb::CellHit> allslchits = slice->AllCells();
      
      if(slice->IsNoise()) continue;

      numusand::NumuSandObj sandbox; // initialization in constructor

      float missE;

      fNumuSandFxs->getMissingE(slice, missE);

      sandbox.SetMissingE(missE);

      sandbox.SetNDropped(livegeom->NDropouts());
      sandbox.SetNBad(livegeom->NBadDCMBC());

      float xmin = 9999, ymin = 9999, zmin = 9999, xmax = -9999, ymax = -9999, zmax = -9999;
      float xmin2 = 9999,  ymin2 = 9999,  zmin2 = 9999,  xmint = 9999,  ymint = 9999,  zmint = 9999;
      float xmax2 = -9999, ymax2 = -9999, zmax2 = -9999, xmaxt = -9999, ymaxt = -9999, zmaxt = -9999;
      int ntx = 0, nty = 0;


      // Temporary code.  Not a great place to put these and not much time to
      // re-order anything.  Calculate vtx E within 20/40/60 cm to origin.
      art::FindManyP<rb::Vertex> fmElastic(slicevec, evt,fVtxLabel);
      if(fmElastic.isValid()){
        std::vector< art::Ptr<rb::Vertex> > elastics = fmElastic.at(i);
        if (elastics.size() == 1){
          // Now, we can get the vtx, and loop over cells to get vars
          TVector3 vtx = elastics[0]->GetXYZ();
          float vtxE20 = 0;
          float vtxE40 = 0;
          float vtxE60 = 0;

          for (unsigned int hIdx = 0; hIdx < slice->NCell(); hIdx++){
            rb::RecoHit rhit = slice->RecoHit(hIdx);
            if (!rhit.IsCalibrated()) continue;
            bool isX = slice->Cell(hIdx)->View()==geo::kX;
            // Calcualte dist from hit to vtx.
            // Only calculate using coordinate cell measures! (DocDB 14330)
            double d =
              isX ? util::sqr(rhit.X()-vtx[0]) :util::sqr(rhit.Y()-vtx[1]);
            d += util::sqr(rhit.Z()-vtx[2]);
            d = sqrt(d);
            if (d < 20){
              vtxE20 += rhit.GeV();
              vtxE40 += rhit.GeV();
              vtxE60 += rhit.GeV();
            }
            else if (d < 40){
              vtxE40 += rhit.GeV();
              vtxE60 += rhit.GeV();
            }
            else if (d < 60)
              vtxE60 += rhit.GeV();
          } // End loop over NCells
          sandbox.SetVtxE20(vtxE20);
          sandbox.SetVtxE40(vtxE40);
          sandbox.SetVtxE60(vtxE60);
        } // End NVtx == 1
      } // End isValid



      for(unsigned int j = 0; j < slice->NCell(); ++j) {
        const rb::CellHit* h = slice->Cell(j).get();
        double xyz[3];
        geom->Plane(h->Plane())->Cell(h->Cell())->GetCenter(xyz);
        const double x = xyz[0];
        const double y = xyz[1];
        const double z = xyz[2];

        if(h->View()==geo::kX) {
          if(x < xmin) { xmin2 = xmin; xmin = x;}
          else if(x < xmin2) xmin2 = x;
          if(z < zmin) { zmin2 = zmin; zmin = z;}
          else if(z < zmin2) zmin2 = z;
          if(x > xmax) { xmax2 = xmax; xmax = x;}
          else if(x > xmax2) xmax2 = x;
          if(z > zmax) { zmax2 = zmax; zmax = z;}
          else if (z > zmax2) zmax2 = zmax;

          if(h->PE() > 100) { // above threshold only TODO: add tracks and recohits, use pecorr
            ntx++;
            if(x < xmint) xmint = x;
            if(z < zmint) zmint = z;

            if(x > xmaxt) xmaxt = x;
            if(z > zmaxt) zmaxt = z;
          }
        }

        else if(h->View()==geo::kY) {
          if(y < ymin) { ymin2 = ymin; ymin = y;}
          else if (y < ymin2) ymin2 = y;
          if(z < zmin) { zmin2 = zmin; zmin = z;}
          else if (z < zmin2) zmin2 = z;

          if(y > ymax) { ymax2 = ymax; ymax = y;}
          else if(y > ymax2) ymax2 = y;
          if(z > zmax) { zmax2 = zmax; zmax = z;}
          else if(z > zmax2) zmax2 = z;

          if(h->PE() > 100) { // above threshold only TODO: add tracks and recohits, use pecorr
            nty++;
            if(y < ymint) ymint = y;
            if(z < zmint) zmint = z;

            if(y > ymaxt) ymaxt = y;
            if(z > zmaxt) zmaxt = z;
          }
        }

      }
      sandbox.SetXMin2(xmin2);
      sandbox.SetYMin2(ymin2);
      sandbox.SetZMin2(zmin2);
      sandbox.SetXMax2(xmax2);
      sandbox.SetYMax2(ymax2);
      sandbox.SetZMax2(zmax2);

      if(ntx > 0 && nty > 0) {
        sandbox.SetXMint(xmint);
        sandbox.SetYMint(ymint);
        sandbox.SetZMint(zmint);
        sandbox.SetXMaxt(xmaxt);
        sandbox.SetYMaxt(ymaxt);
        sandbox.SetZMaxt(zmaxt);
      }

      else { // -5 default can be confused as a real value; +=9999 too far from real values
        sandbox.SetXMint(-999);
        sandbox.SetYMint(-999);
        sandbox.SetZMint(-999);
        sandbox.SetXMaxt(-999);
        sandbox.SetYMaxt(-999);
        sandbox.SetZMaxt(-999);
      }


      if(!trackAssnList.isValid()) {
        std::cout << "NumuSandbox: No Kalman Tracks!" << std::endl;
        outputObjects->push_back(sandbox);
        util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
        continue;
      }

      const std::vector< art::Ptr<rb::Track> > sliceTracks = trackAssnList.at(i);

      art::FindOneP<remid::ReMId> remidVec(sliceTracks, evt, fPIDLabel);
      if(!remidVec.isValid()) {
        //        std::cout << "NumuSandbox: No ReMId information!" << std::endl;
        outputObjects->push_back(sandbox);
        util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
        continue;
      }


      int bestIdx = -1;  // for index in list of the highest PID track
      int pimuIdx = -5;  // for index in list of a pimu decay event
      int Idx2 = -1;
      int idx3d = -1;
      double bestpimuE = -999.;  // keep track of how much energy due to a non-primary muon was in a track
                         // this is for if a pi->mu muon contributes light to more than once track
                         // we only want to associate with the track that saw 'most' of it.
                         // note we assume all non-primary muons are from pi->mu decay.
      double bestval = -999., val2 = -999;
      int nmutrks = 0;

      // ----------------------------------
      // WORKING AREA
      
      float avededxtrk1 = 0.0;
      float avededxtrk2 = 0.0;
      float avededxtrk1Last4Cells = 0.0;
      float avededxtrk2Last4Cells = 0.0;
      float avededxtrk1Last6Cells = 0.0;
      float avededxtrk2Last6Cells = 0.0;
      float avededxtrk1Last8Cells = 0.0;
      float avededxtrk2Last8Cells = 0.0;
      float scattangletrk1 = 0.0;
      float scattangletrk2 = 0.0;
      float actvtxx = 0.0;
      float slccalE = 0.0;
      
      if((sliceTracks.size() == 2)){
        // Loop over tracks in slice
        for(size_t j = 0; j < sliceTracks.size(); ++j){
          const art::Ptr<rb::Track> track = sliceTracks[j];
          art::Ptr<remid::ReMId> trackRemid = remidVec.at(j);
          
          if(!track->Is3D()) continue;
          
          if(trackRemid->Value() > bestval) { // best track so far
            val2 = bestval;
            Idx2 = bestIdx;
            bestval = trackRemid->Value();
            bestIdx = j;
          }
          else if(trackRemid->Value() > val2) { // not best, but better second best track
            val2 = trackRemid->Value();
            Idx2 = j;
          }
          
        } // end loop over tracks
        
        // Making sure is OK
        if((!remidVec.isValid()) ||
           (bestval < 0.0) || (val2 < 0.0)){
          continue;
        }
        
        // Get info on muon like track
        //if ((sliceTracks[bestval]->NCell()>0) && (bestval>0.75)){
        // get info about this track
        const art::Ptr<rb::Track> track1 = sliceTracks[bestIdx]; 
        //
        // Fill here with dEdx info and other variables
        //
        fNumuSandFxs->getAveTrackdEdx(track1,avededxtrk1);
        sandbox.SetAvedEdxTrack1(avededxtrk1);
        
        fNumuSandFxs->getAveTrackdEdxLast4Cells(track1,avededxtrk1Last4Cells);
        sandbox.SetAvedEdxTrack1Last4Cells(avededxtrk1Last4Cells);
        
        fNumuSandFxs->getAveTrackdEdxLast6Cells(track1,avededxtrk1Last6Cells);
        sandbox.SetAvedEdxTrack1Last6Cells(avededxtrk1Last6Cells);
        
        fNumuSandFxs->getAveTrackdEdxLast8Cells(track1,avededxtrk1Last8Cells);
        sandbox.SetAvedEdxTrack1Last8Cells(avededxtrk1Last8Cells);
        
        scattangletrk1 = fNumuSandFxs->getScatt(track1);
        if(track1->TotalLength() > 0) sandbox.SetScattAngleTrack1(scattangletrk1/track1->TotalLength());
        
        //}
        //else if((sliceTracks[val2]->NCell()>0)){
        // get info about this track
        const art::Ptr<rb::Track> track2 = sliceTracks[Idx2]; 
        //
        // Fill here with dEdx info and other variables
        //
        fNumuSandFxs->getAveTrackdEdx(track2,avededxtrk2);
        sandbox.SetAvedEdxTrack2(avededxtrk2);
        
        fNumuSandFxs->getAveTrackdEdxLast4Cells(track2,avededxtrk2Last4Cells);
        sandbox.SetAvedEdxTrack2Last4Cells(avededxtrk2Last4Cells);
        
        fNumuSandFxs->getAveTrackdEdxLast6Cells(track2,avededxtrk2Last6Cells);
        sandbox.SetAvedEdxTrack2Last6Cells(avededxtrk2Last6Cells);
        
        fNumuSandFxs->getAveTrackdEdxLast8Cells(track2,avededxtrk2Last8Cells);
        sandbox.SetAvedEdxTrack2Last8Cells(avededxtrk2Last8Cells);
        
        scattangletrk2 = fNumuSandFxs->getScatt(track2);
        if(track2->TotalLength() > 0) sandbox.SetScattAngleTrack2(scattangletrk2/track2->TotalLength());
        //}
        
        ///////////////////////////////////
        ///--------------------------------
        ///  ADD HERE VARS FOR QE PID?
        ///-------------------------------
        /////////////////////////////////// 
                
        // WORKING AREA 
        //art::PtrVector<rb::CellHit> allslchits = slice->AllCells();
        art::PtrVector<rb::CellHit> alltrk1hits = track1->AllCells();
        art::PtrVector<rb::CellHit> alltrk2hits = track2->AllCells();
        
        fNumuSandFxs->getActivityVtx(alltrk1hits,alltrk2hits,allslchits,actvtxx);
        slccalE =  slice->CalorimetricEnergy();
        
        if(slice->CalorimetricEnergy() <= 0){
          //continue;
          std::cout<<"Calorimetric energy is zero!!"<<std::endl;
        }
        else{
          sandbox.SetActVtx(actvtxx/slccalE);
          //std::cout<<"TMVAvars[0]: "<< actvtxx/slccalE <<std::endl;
          //TMVAvars[0] = actvtx/slccalE;
        }

        //---- 
        // Another working area
        // from QEEventFinder/QeFinder_module.cc
        // get numu energy estimates
        float trkCCE = -5.0;
        float RecoMuonE = -5.0;
        float Reco2ndTrkE = -5.0;
        
        if(eAss.isValid()){
          std::vector<art::Ptr<numue::NumuE> > numuEnergy = eAss.at(i);
          if(numuEnergy.size() > 0){
            // get the energy variables
            trkCCE = numuEnergy[0]->TrkCCE();
          }//end if valid size
          
          art::FindOneP<rb::Energy>   trackEnergy(sliceTracks, evt, fEnergyLabel);
          if(trackEnergy.isValid()){
            RecoMuonE = trackEnergy.at(bestIdx)->E();
          }
          //Reco2ndTrkE = track2->TotalGeV(); // sloppy approach
          Reco2ndTrkE = track2->CalorimetricEnergy(); // sloppy approach
          
          float tmpHadE = trkCCE - RecoMuonE;
          
          float offTrkEnFrac = (tmpHadE - Reco2ndTrkE)/trkCCE;
          
          sandbox.SetOffTrkFra(offTrkEnFrac);

          //std::cout<<"TMVAvars[0]: "<< offTrkEnFrac <<std::endl;
          
        }// eAss is valid
        else{
          std::cout<<"Event Association is not Valid!"<<std::endl;
        }
        
      }//end slice tracks = 2
      
      // ----------------------------------
      
      // Loop over tracks in slice
      for(size_t j = 0; j < sliceTracks.size(); ++j){
        const art::Ptr<rb::Track> track = sliceTracks[j];
        art::Ptr<remid::ReMId> trackRemid = remidVec.at(j);

        if(!track->Is3D()) {
          continue;
        }

        idx3d++;
        
        if(trackRemid->Value() > bestval) { // best track so far
          val2 = bestval;
          Idx2 = bestIdx;
          bestval = trackRemid->Value();
          bestIdx = j;
        }
        else if(trackRemid->Value() > val2) { // not best, but better second best track
          val2 = trackRemid->Value();
          Idx2 = j;
        }

        if(trackRemid->Value() > 0.6) nmutrks++;

        if(bt->HaveTruthInfo()) {
          art::PtrVector<rb::CellHit> trkhits = track->AllCells();
          const std::vector<const sim::Particle*> parts = bt->HitsToParticle(trkhits);
          const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(trkhits);

          // check particle that contributes most to track for pi->mu decay
          if( !parts.empty() && !trids.empty() ) {
            if(abs(parts[0]->PdgCode())==13 && parts[0]->Process()=="Decay") {
              double tempE = trids[0].energyFrac * track->TotalGeV();
              if(tempE > bestpimuE) {
                bestpimuE = tempE;
                pimuIdx = idx3d+1;
              }
            }
          }

          // check particle that contributes second most; track can contain both pi and mu
          if(parts.size() > 1 && trids.size() > 1) {
            if(abs(parts[1]->PdgCode())==13 && parts[1]->Process()=="Decay") {
              double tempE = trids[1].energyFrac * track->TotalGeV();
              if(tempE > bestpimuE) {
                bestpimuE = tempE;
                pimuIdx = idx3d+1;
              }
            }
          }
        }
      } //end loop over tracks

      // ===============================================================
      //New: Get info on muon like track
      //if( (bestval > 0.0) && (val2 > 0.0) ){
      if((bestval > 0.0)){
        
        const art::Ptr<rb::Track> mutrack = sliceTracks[bestIdx];
        art::PtrVector<rb::CellHit> mutrkhits = mutrack->AllCells();
        
        //Getting all the hits in the slice that are not on our best track
        ////rb::Cluster vertexHits = MakeVertexCluster(allslchits,mutrkhits);
        rb::Cluster vertexCluster;
        int nhadhits = -5;
        for (unsigned int sliceHit = 0; sliceHit < allslchits.size(); ++sliceHit){
          
          bool match = 0;
          
          for (unsigned int trackHit = 0; trackHit < mutrkhits.size(); ++trackHit){
            match = (allslchits[sliceHit] == mutrkhits[trackHit]);
            if (match) break;
          }//End of loop over hits in track   
          
          if (!match) vertexCluster.Add(allslchits[sliceHit]);
          
        }//End of loop over hits in the slice
        
        nhadhits = vertexCluster.NCell();
        //std::cout<<" ==== No. hits in had system: "<<vertexCluster.NCell()<<std::endl;
        
        sandbox.SetNHadHits(nhadhits);
        
        // New
        unsigned int minCellsFromEdge = 99999999;
        //int  cellsFromEdge = -1;
        
        // TO DO: case when nhadhits is zero?
        // Loop over vertex cluster
        for(unsigned int hitIdx = 0; hitIdx < vertexCluster.NCell(); ++hitIdx){
          const art::Ptr<rb::CellHit>& chit = vertexCluster.Cell(hitIdx);
          const int planeNum = chit->Plane();
          const unsigned int cellsFromEdge = std::min((unsigned int)chit->Cell(), geom->Plane(planeNum)->Ncells() - 1 - chit->Cell());
          
          if(cellsFromEdge < minCellsFromEdge){
            minCellsFromEdge = cellsFromEdge;
          } //end if
          
        } //end for        
        //std::cout<<" ==== Had ncells from edge: "<< minCellsFromEdge  <<std::endl;
        
        sandbox.SetHadCellsFromEdge(minCellsFromEdge);
        
      }//end if best val
      // ===============================================================
      
      sandbox.SetNMuTrks(nmutrks);
      
      if(bt->HaveTruthInfo()) {  // no truth info case should correspond to unfilled

        if(bestIdx!=-1) {

          const art::Ptr<rb::Track> track1 = sliceTracks[bestIdx];
          const std::vector<const sim::Particle*> parts1 = bt->HitsToParticle(track1->AllCells());
          const std::vector< cheat::TrackIDE > trids1 = bt->HitsToTrackIDE(track1->AllCells());

          if(!parts1.empty()){
            sandbox.SetPDGBest(parts1[0]->PdgCode());
          }
        }

        if(Idx2!=-1) {

          const art::Ptr<rb::Track> track2 = sliceTracks[Idx2];
          const std::vector<const sim::Particle*> parts2 = bt->HitsToParticle(track2->AllCells());
          const std::vector< cheat::TrackIDE > trids2 = bt->HitsToTrackIDE(track2->AllCells());

          if(!parts2.empty()){
            sandbox.SetPDGSecondBest(parts2[0]->PdgCode());
          }
        }

        int nprotons = 0;
        if(pimuIdx == -5) pimuIdx = 0;  // differentiate unfilled with not found
        art::PtrVector<rb::CellHit> slchits = slice->AllCells();
        const std::vector<const sim::Particle*> parts = bt->HitsToParticle(slchits);
        const std::vector< cheat::TrackIDE > trids = bt->HitsToTrackIDE(slchits);
        if(!parts.empty() && !trids.empty()) {
          for(unsigned int i = 0; i < parts.size(); i++) {
            if(trids.size() > i) {
              if(parts[i]->PdgCode()==2212 && trids[i].energyFrac * slice->TotalGeV() > 0.05) nprotons++;
            }
          }
        }
        sandbox.SetNProtons(nprotons);
      }

      sandbox.SetPiMuDecay(pimuIdx);

      outputObjects->push_back(sandbox);
      util::CreateAssn(*this,evt,*(outputObjects.get()),slice,*(assoc.get()));
    } //end loop over slices

    evt.put(std::move(outputObjects));
    evt.put(std::move(assoc));

  } //end produce

} //end namespace


DEFINE_ART_MODULE(numusand::FillSandbox)