ReMIdTrain_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////           
// \file    ReMIdTrain_module.cc
// \brief   A module for creating trees which can be used for kNN training and generating sample histograms for ReMId
// \version $Id: ReMIDTrain_module.cc,v 2013/02/21 16:30:29 raddatz Exp $  
// \author  Nicholas Raddatz - raddatz@physics.umn.edu                                                    
//////////////////////////////////////////////////////////////////// 
#include <string>

// ROOT includes
#include "TVector3.h"
#include "TTree.h"
#include "TMath.h"

// Framework includes
#include "fhiclcpp/ParameterSet.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "Utilities/AssociationUtil.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Core/EDAnalyzer.h"

// NOvA includes
#include "Geometry/Geometry.h"
#include "Geometry/LiveGeometry.h"
#include "RunHistory/service/RunHistoryService.h"
#include "ReMId/ReMId.h"
#include "GeometryObjects/Geo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "GeometryObjects/CellGeo.h"
#include "MCCheater/BackTracker.h"
#include "RecoBase/Track.h"
#include "Simulation/Simulation.h"


namespace remid {

  // A module to analyze ReMId objects for pid training
  class ReMIdTrain : public art::EDAnalyzer {
  public:
    
    explicit ReMIdTrain(fhicl::ParameterSet const& pset);
    virtual ~ReMIdTrain();

    void beginJob();
    void analyze(art::Event const& evt);
    void endJob();

  private:
    std::string fSliceLabel;         ///<Where to find reconstructed slices
    std::string fTrackLabel;         ///<Where to find recondtructed tracks
    std::string fPidLabel;           ///<Where to find the pids
    bool fCreateTrainTrees;          ///<Create kNN training tree?
    bool fCreateSampleTrees;         ///<Create trees for sample histograms?

    TTree* TreeS;  ///< Tree for holding signal input variables for kNN training
    TTree* TreeB;  ///< Tree for holding background input variables for kNN training

    // variables used in filling the kNN training trees
    float fTrackLength;
    float fDedxSep;
    float fScatSep;
    float fnuE;
    float fvisE;
    float frecoE;
    float fMeasFrac;

    TTree* DedxSample;  ///< Tree for holding dE/dx variables needed for making signal sample histograms 
    TTree* ScatSample;  ///< Tree for holding scattering variables needed for making signal sample histograms

    // variables used in filling the sample histogram trees
    float fDedx;
    int fnDedx;
    int fnDedxUsed;
    float fDistDedx;
    float fCosScat;
    float fDistLastScat;
    float fScatVar;
    int fnScat;
    float fDistScat;
    int fnPlane;
    int fpdg;
    int fintType;
    bool fVert;
    bool fUsed;

    art::ServiceHandle<geo::Geometry> fgeom;
    art::ServiceHandle<nova::dbi::RunHistoryService> frh;
    art::ServiceHandle<geo::LiveGeometry>  flivegeom;

    unsigned int fdibfirst = 0;
    unsigned int fdiblast = 0;

    bool ContainedEvent(art::Ptr<rb::Cluster> slice, art::Ptr<rb::Track> trk);

    static bool SortByHits(const art::Ptr<rb::Track>&  a, const art::Ptr<rb::Track>&  b);


  };

}


namespace remid{

  //.......................................................................

  struct Counts
  {
    Counts()
    {
      nHits[0] = nHits[1] = 0;
      nOnlyHits[0] = nOnlyHits[1] = 0;
    }
    int nHits[2]; // Indexed by view
    int nOnlyHits[2]; // Hit by only one particle
  };

  //.......................................................................

  ReMIdTrain::ReMIdTrain(fhicl::ParameterSet const& pset):
    EDAnalyzer(pset)
  {
    fSliceLabel = (pset.get< std::string >("SliceLabel"));
    fTrackLabel = (pset.get< std::string >("TrackLabel"));
    fPidLabel   = (pset.get< std::string >("PidLabel"  ));
    fCreateTrainTrees  = (pset.get< bool >("CreateTrainTrees"));
    fCreateSampleTrees = (pset.get< bool >("CreateSampleTrees"));
  }

  //.......................................................................

  ReMIdTrain::~ReMIdTrain()
  {
  }

  //.......................................................................

  void ReMIdTrain::beginJob()
  {
    art::ServiceHandle<art::TFileService> tfs;

    if(fCreateTrainTrees){
      TreeS = tfs->make<TTree>("fTreeS","kNN Input Variables for Signal");
      TreeS->Branch("trackLength", &fTrackLength,"trackLength/F");
      TreeS->Branch("dedxSep", &fDedxSep, "dedxSep/F");
      TreeS->Branch("scatSep", &fScatSep, "scatSep/F");
      TreeS->Branch("nuE", &fnuE, "nuE/F");
      TreeS->Branch("intType", &fintType,"intType/I");    
      TreeS->Branch("visE", &fvisE, "visE/F");   
      TreeS->Branch("recoE", &frecoE, "recoE/F");
      TreeS->Branch("measFrac", &fMeasFrac, "measFrac/F");

      TreeB = tfs->make<TTree>("fTreeB","kNN Input Variables for Background");
      TreeB->Branch("trackLength", &fTrackLength,"trackLength/F");
      TreeB->Branch("dedxSep", &fDedxSep, "dedxSep/F");
      TreeB->Branch("scatSep", &fScatSep, "scatSep/F");
      TreeB->Branch("nuE", &fnuE, "nuE/F");
      TreeB->Branch("intType", &fintType, "intType/I");
      TreeB->Branch("pdg", &fpdg, "pdg/I");
      TreeB->Branch("visE", &fvisE, "visE/F");   
      TreeB->Branch("recoE", &frecoE, "recoE/F");
      TreeB->Branch("measFrac", &fMeasFrac, "measFrac/F");


    }

    if(fCreateSampleTrees){
      DedxSample = tfs->make<TTree>("fDedxSample","Variables for Creating Dedx Sample Histograms");
      DedxSample->Branch("dedx", &fDedx, "dedx/F");
      DedxSample->Branch("nDedx", &fnDedx, "nDedx/I");
      DedxSample->Branch("nDedxUsed", &fnDedxUsed, "nDedxUsed/I");
      DedxSample->Branch("distFromEnd", &fDistDedx, "distFromEnd/F");
      DedxSample->Branch("trackLength", &fTrackLength, "trackLength/F");
      DedxSample->Branch("nPlane", &fnPlane, "nPlane/I");
      DedxSample->Branch("pdg", &fpdg, "pdg/I");
      DedxSample->Branch("intType", &fintType, "intType/I");
      DedxSample->Branch("vert", &fVert, "vert/O");
      DedxSample->Branch("used", &fUsed, "used/O");
      DedxSample->Branch("nuE", &fnuE, "nuE/F");

      ScatSample = tfs->make<TTree>("fScatSample","Variables for Creating Scattering Sample Histograms");
      ScatSample->Branch("cosScat", &fCosScat, "cosSact/F");
      ScatSample->Branch("distFromLastScat", &fDistLastScat, "distFromLastScat/F");
      ScatSample->Branch("scatVar", &fScatVar, "scatVar/F");
      ScatSample->Branch("nScat", &fnScat, "nScat/I");
      ScatSample->Branch("distFromEnd", &fDistScat, "distFromEnd/F");
      ScatSample->Branch("trackLength", &fTrackLength, "trackLength/F");
      ScatSample->Branch("nPlane", &fnPlane, "nPlane/I");
      ScatSample->Branch("pdg", &fpdg, "pdg/I");
      ScatSample->Branch("intType", &fintType, "intType/I");

    }

  }

  //.......................................................................

  void ReMIdTrain::analyze(art::Event const& evt)
  {
    // get the slices out of the event
    art::Handle<std::vector<rb::Cluster> > slicecol;
    evt.getByLabel(fSliceLabel,slicecol);

    // // get the tracks out of the event
    // art::Handle<std::vector<rb::Track> > trkcol;
    // evt.getByLabel(fTrackLabel,trkcol);

    art::PtrVector<rb::Cluster> slicelist;
    art::PtrVector<rb::CellHit> Allevthits; // PtrVector of all of the hits in the event
    
    for(unsigned int i = 0; i<slicecol->size();++i){
      art::Ptr<rb::Cluster>slice(slicecol,i);
      slicelist.push_back(slice);
      for(unsigned int ihit = 0; ihit<slice->NCell();++ihit){
        Allevthits.push_back(slice->Cell(ihit));
      }
    }

    // get  assosciations between slices and tracks
    art::FindManyP<rb::Track> fndmnytrk(slicecol,evt,fTrackLabel);

    // get the detector configuration
    unsigned int dibmask = (frh->GoodDiBlockMask()&frh->GetConfiguration());
    if(dibmask == ((1<<16)-1)){ dibmask = frh->GetConfiguration(); }
    if(dibmask == 0){ dibmask = frh->GetConfiguration(); }
    if(dibmask){
      int iD;
      iD=0;
      while(!((dibmask>>iD)&1)){
        iD++; 
        fdibfirst=iD+1;  
      }
      iD=0; 
      while(dibmask>>iD){
        iD++; 
        fdiblast=iD;
      }
    }


    art::ServiceHandle<cheat::BackTracker> bt;

    // loop over all slices and find the most effiecient slice for a neutrino interaction
    double bestEff = 0;
    unsigned int mostEffSlice = 0;  // Index of the most efficient neutrino slice
    int nuintType = 6;  // This is a index to keep track of what type of interaction this is
    // get the neutrino interaction type
    for(unsigned int iSlice = 0;iSlice<slicelist.size();++iSlice){
      if(slicelist[iSlice]->IsNoise()){ continue; }
      // get the neutrino interaction from this slice
      std::vector<cheat::NeutrinoEffPur> nuint = bt->SliceToNeutrinoInteractions(slicelist[iSlice]->AllCells(),Allevthits);
      if(nuint.size() > 0 && nuint[0].efficiency > bestEff){
        bestEff = nuint[0].efficiency;
        mostEffSlice = iSlice;
        // get the neutrino
        const simb::MCNeutrino nu = nuint[0].neutrinoInt->GetNeutrino();
        fnuE = nu.Nu().E();
        fvisE = nuint[0].energySlice;
        if(nu.CCNC() == simb::kNC){ nuintType = 5; }
        else{
          if(nu.Nu().PdgCode() == 12 || nu.Nu().PdgCode() == -12){nuintType = 4;}
          else if(nu.Mode() == 0){nuintType = 0;}
          else if(nu.Mode() == 1){nuintType = 1;}
          else if(nu.Mode() == 2){nuintType = 2;}
          else if(nu.Mode() == 3){nuintType = 3;}
        }
      }
    }

    // loop over the slices
    for(unsigned int iSlice = 0; iSlice<slicelist.size(); ++iSlice){
      // if we are not on the most efficient slice don't try to do anything
      if(iSlice != mostEffSlice){ continue; }
      // also if the efficiency is not greater than 1 or is a noise slice skip it
      if(bestEff == 0){ continue; }
      // If slice is noise, skip it
      if(slicelist[iSlice]->IsNoise()){ continue; }
      // If don't know what neutrino this slice comes from skip it
      if(nuintType == 6){ continue; }

      // get all of the tracks associated to this slice
      std::vector<art::Ptr<rb::Track> > tracks = fndmnytrk.at(iSlice);
      // if no tracks move on
      if(tracks.size() == 0){ continue; }

      frecoE = slicelist[iSlice]->TotalGeV();

      // sort the tracks by the number of hits
      std::sort(tracks.begin(),tracks.end(),SortByHits);

      //get the trackIDs for the event
      std::set<int> trackIDs = bt->GetTrackIDList();

      // don't penalize purity because of michel electron
      std::set<int> allowedDaughters;
      allowedDaughters.insert(-11);
      allowedDaughters.insert(11);
            
      // // get the id of the track to use as either signal or background
      // int trackid = -1;
      // // get the reco track id of the particle
      // int recotrackid = -1;
      
      std::set<int> primTrkIDs;
      // Get the set of primary particle track ids for this event
      for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end();++it){
        if(*it == bt->TrackIDToMotherParticle(*it)->TrackId()){ primTrkIDs.insert(*it); }
      }    
      
      // match tracks to these primary particles
      std::vector<cheat::ParticleEffPur> pep = bt->TracksToParticles(tracks, primTrkIDs, allowedDaughters, slicelist[iSlice]->AllCells());

      // get the pid that belongs to this recotrack
      art::FindOneP<remid::ReMId> foids(tracks,evt,fPidLabel);

      // loop over all of the tracks
      for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
        
        fpdg = pep[itrk].pdg;
        fintType = nuintType;
        
        // check if this track is matched to a primary particle
        if(fpdg == 0){ continue; }
        
        // Be safe, check that efficiency and purity are both greater than 0
        if(!(pep[itrk].efficiency > 0 || pep[itrk].purity > 0)){ continue; }
        
        // get the pid that belongs to this track
        art::Ptr<remid::ReMId> pid = foids.at(itrk);
                
        fnPlane = tracks[itrk]->ExtentPlane();
        
        // fill trees with variable for creating sample histograms if wanted
        if(fCreateSampleTrees){

          // Create sample histograms from contained tracks only
          
          // check if this track is contained by pid
          if(!(pid->Contained())){ continue; }
          
          // any other track cuts?
          fTrackLength = tracks[itrk]->TotalLength();
          
          // cut out strange tracklengths
          if(tracks[itrk]->TotalLength() > 6500){ continue; }

          // make sure its a 3d track
          if(!tracks[itrk]->Is3D()){ continue; }

          // get all the dedx variables filled
          fnDedx = pid->NDedx()-1;
          fnDedxUsed = pid->NMeasurementPlanes();
          // loop over all of the dE/dx measurements
          for(unsigned int idedx = 0; idedx<pid->NDedx(); ++idedx){
            fDedx = pid->Dedx(idedx);
            fVert = pid->DedxVertex(idedx);
            fUsed = pid->DedxUsed(idedx);
            // get the distance from the end of the track
            fDistDedx = (1.0 - pid->DedxLocation(idedx))*tracks[itrk]->TotalLength();
            // fill the dedx sample tree
            DedxSample->Fill();
          }
          
          // get all the scattering variables filled
          fnScat = pid->NScatters();
          for(unsigned int iscat = 0; iscat<pid->NScatters();++iscat){
            fCosScat = pid->CosScat(iscat);
            double ang = TMath::ACos(fCosScat);
            // get the distance from the end of the track and from the last scattering measurement
            TVector3 trkstart = tracks[itrk]->Start();
            double dist = 0;
            double lastdist = 0;
            for(unsigned int itraj = 0; itraj<tracks[itrk]->NTrajectoryPoints()-1; ++itraj){
              TVector3 p0 = tracks[itrk]->TrajectoryPoint(itraj);
              TVector3 p1 = tracks[itrk]->TrajectoryPoint(itraj+1);
              lastdist = (p1-p0).Mag();
              dist+=lastdist;
              if(p1 == pid->ScatLocation(iscat)){break;}
            }
            fDistLastScat = lastdist;
            fDistScat = tracks[itrk]->TotalLength()-dist;
            fScatVar = ang*ang/lastdist;
            // fill the scat sample tree
            ScatSample->Fill();
          }
          
        } // Fill sample histogram trees
        

      } // end of loop over reco tracks


      if(fCreateTrainTrees){

        fintType = nuintType;
        if(fintType < 4){
          // find the best muon track
          int bestMu = -1;
          for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
            if(abs(pep[itrk].pdg) == 13){ bestMu = itrk; }
          }
          if(bestMu >= 0){
            // get the pid that belongs to this track
            art::Ptr<remid::ReMId> pid = foids.at(bestMu);
            // check if this pid has any usable measurements of dedx and that the track is 3d
            if(pid->NMeasurementPlanes() > 0 && tracks[bestMu]->Is3D()){
              fMeasFrac = ((float)pid->NMeasurementPlanes())/((float)pid->NMeasurementPlanes()+(float)pid->NVertexPlanes());
              // check if this track is contained by pid
              // check if this event is contained
              if(ContainedEvent(slicelist[iSlice],tracks[bestMu])){
                // any other track cuts?
                fTrackLength = tracks[bestMu]->TotalLength();
                // cut out strange tracklengths
                if(fTrackLength < 6500){
                  fDedxSep = pid->DedxSeparation();
                  fScatSep = pid->ScatSeparation();
                  TreeS->Fill();
                }
              } // if contained
            }
          } // if muon matched track
        } // CC event
        else{
          // non CC event
          int bestPi = -1;
          double longestPi = -1;
          int longestTrkIdx = -1;
          double longestTrk = -1;
          // find the most muon like track
          // first look for a primary pion match
          for(unsigned int itrk = 0; itrk < tracks.size(); ++itrk){
            double trklen = tracks[itrk]->TotalLength();
            if(abs(pep[itrk].pdg) == 211 && trklen > longestPi){
              bestPi = itrk;
              longestPi = trklen;
            }
            if(trklen > longestTrk){
              longestTrkIdx = itrk;
              longestTrk = trklen;
            }
          } // end of loop over tracks
          
            // train against longest pion track if it exist
          if(bestPi >= 0){
            // get the pid that belongs to this track
            art::Ptr<remid::ReMId> pid = foids.at(bestPi);
            // check if this pid has any usable measurements of dedx and that the track is 3d
            if(pid->NMeasurementPlanes() > 0 && tracks[bestPi]->Is3D()){
              fMeasFrac = ((float)pid->NMeasurementPlanes())/((float)pid->NMeasurementPlanes()+(float)pid->NVertexPlanes());
              // check if this track is contained by pid
              if(ContainedEvent(slicelist[iSlice],tracks[bestPi])){
                fTrackLength = longestPi;
                // cut out strange tracklengths
                if(fTrackLength < 6500){
                  fpdg = pep[bestPi].pdg;
                  fDedxSep = pid->DedxSeparation();
                  fScatSep = pid->ScatSeparation();
                  TreeB->Fill();
                }
              } // if contained
            }
          } // if pion exists
          else if(longestTrkIdx >= 0){
            // get the pid that belongs to this track
            art::Ptr<remid::ReMId> pid = foids.at(longestTrkIdx);
            // check if this pid has any usable measurements of dedx and that the track is 3d
            if(pid->NMeasurementPlanes() > 0 && tracks[longestTrkIdx]->Is3D()){
              fMeasFrac = ((float)pid->NMeasurementPlanes())/((float)pid->NMeasurementPlanes()+(float)pid->NVertexPlanes());
              // check if this track is contained by pid
              if(ContainedEvent(slicelist[iSlice],tracks[longestTrkIdx])){
                fTrackLength = longestTrk;
                // cut out strange tracklengths
                if(fTrackLength < 6500){
                  fpdg = pep[longestTrkIdx].pdg;
                  fDedxSep = pid->DedxSeparation();
                  fScatSep = pid->ScatSeparation();
                  TreeB->Fill();
                }
              } // if contained
            }
          } // if long track exists
          
        } // NC event
        
      } // if fill trees for kNN training


    } // end loop over slices


  }

  //.......................................................................

  void ReMIdTrain::endJob()
  {
  }

  //.......................................................................

  bool ReMIdTrain::SortByHits(const art::Ptr<rb::Track>& a, const art::Ptr<rb::Track>& b)
  {
    return a->NCell()>b->NCell();
  }

  //.......................................................................
  bool ReMIdTrain::ContainedEvent(art::Ptr<rb::Cluster> slice, art::Ptr<rb::Track> trk)
  {
    
    bool contained = false;
    
    // check that this is a contained event
    
    // get some useful information here
    // 1) how close is the slice to the edge of the detector
    unsigned int slcncellsfromedge = 999999999;
    for(unsigned int ihit = 0; ihit < slice->NCell(); ++ihit){
      const art::Ptr<rb::CellHit>& chit = slice->Cell(ihit);
      const int plane = chit->Plane();
      unsigned int cellsfromedge = std::min((unsigned int)chit->Cell(), fgeom->Plane(plane)->Ncells() - 1 - chit->Cell());
      if(cellsfromedge < slcncellsfromedge){
        slcncellsfromedge = cellsfromedge;
       }
    }
    // 2) how close is the track to the edge of the detector
    // in planes
    unsigned int firstpl = trk->MinPlane();
    unsigned int lastpl = trk->MaxPlane();
    // in projected cells, just define directions here actually get projected cell below
    TVector3 dirbwk = -trk->Dir();
    TVector3 dirfwd = (trk->TrajectoryPoint(trk->NTrajectoryPoints()-1) - trk->TrajectoryPoint(trk->NTrajectoryPoints()-2)).Unit();
    
    if(fgeom->DetId() == novadaq::cnv::kFARDET){
      // cut the slice number of cells from edge
      if(slcncellsfromedge > 1){
        unsigned int npltofront = firstpl - (fdibfirst-1)*64;
        unsigned int npltoback = (fdiblast*64) - 1 - lastpl;
        // check that the start of track is not the front or back of the detector
        if(npltofront > 1 && npltoback > 1){
          int fwdcell = flivegeom->ProjectedLiveCellsToEdge(trk->Stop(),dirfwd);
          int bwkcell = flivegeom->ProjectedLiveCellsToEdge(trk->Start(),dirbwk);
          if(fwdcell > 10 && bwkcell > 10){
            contained = true;
          } // end if track ends projections contained
        } // end if track end planes contained
      } // end if slice contained
    } // end if far detector
    else if(fgeom->DetId() == novadaq::cnv::kNEARDET){
      // cut the slice number of cells from edge
      std::cout<<" slice cells from edge: "<<slcncellsfromedge<<std::endl;
      if(slcncellsfromedge > 1){
        // check that first diblock and last diblock give you the full detector
        std::cout<<"fdibfirst: "<<fdibfirst<<" fdiblast: "<<fdiblast<<std::endl;
        // if(fdibfirst == 1 && fdiblast == 4){
        // check that the start of track is not the front or back of the detector
        std::cout<<"firstpl: "<<firstpl<<" lastpl: "<<lastpl<<std::endl;
        if(firstpl > 1 && lastpl < 212){
          int fwdcell = flivegeom->FullNDProjectedCells(trk->Stop(), dirfwd);
          int bwkcell = flivegeom->FullNDProjectedCells(trk->Start(),dirbwk);
          // also make sure the track doesn't start too close to the end of the active region
          std::cout<<"fwdcell: "<<fwdcell<<" bwkcell: "<<bwkcell<<" start Z: "<<trk->Start().Z()<<std::endl;
          if(fwdcell > 4 && bwkcell > 8 && trk->Start().Z() < 1150){
            std::cout<<"stop Z: "<<trk->Stop().Z()<<std::endl;
            // if the track ends in the muon catcher, check that it doesn't clip the corner
            if(trk->Stop().Z() < 1275){ 
              std::cout<<"contained"<<std::endl;
              contained = true; 
            }
            else{
              
              // find the trajectory point closest in z to the transition plane
              double mindist = 99999999;
              unsigned int mintrajid = 0;
              for(unsigned int itraj = 0; itraj < trk->NTrajectoryPoints()-1; ++itraj){
                double dist = 1275 - trk->TrajectoryPoint(itraj).Z();
                if(dist > 0 && dist < mindist){
                  mindist = dist;
                  mintrajid = itraj;
                }
              }
              // get the direction at this trajectory point
              TVector3 trandir = trk->Dir();
              if(mintrajid < trk->NTrajectoryPoints() -2){
                TVector3 trandir = (trk->TrajectoryPoint(mintrajid+1) - trk->TrajectoryPoint(mintrajid)).Unit();
              }
              double ypostran = flivegeom->YPositionAtTransition(trk->TrajectoryPoint(mintrajid),trandir);
              std::cout<<"ypostran: "<<ypostran<<std::endl;
              if( ypostran < 55){ 
                std::cout<<"contained"<<std::endl;
                contained = true; 
              }
              
            } // end else in muon catcher
          } // end if track ends projections contained
        } // end if track end planes contained 
        // } // end if  expected number of diblocks
      } // end if slice contained
    } // end if near detector
    
    return contained;
    
  }
  
  
}


namespace remid{

  DEFINE_ART_MODULE(ReMIdTrain)

}