RecoCheckAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief 
/// \author brebel@fnal.gov
/// \date   Jul. 15 2011 
////////////////////////////////////////////////////////////////////////

#include "TH1.h"
#include "TH2.h"

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

#include "MCCheater/BackTracker.h"
#include "RecoBase/Prong.h"
#include "RecoBase/Track.h"
#include "RecoBase/Cluster.h"
#include "RecoBase/CellHit.h"
#include "Simulation/Particle.h"
#include "Simulation/Simulation.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"

namespace cheat {
  class RecoCheckAna;
}

class TH1D;
class TH2D;

class cheat::RecoCheckAna : public art::EDAnalyzer {
public:
  explicit RecoCheckAna(fhicl::ParameterSet const &p);
  virtual ~RecoCheckAna();

  virtual void analyze(art::Event const &e);

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

  

private:
  struct Plots
  {
    Plots() : purity(0), efficiency(0), deltaLength(0),
              purityVsLength(0), efficiencyVsLength(0), deltaLengthVsLength(0),
              effvpur(0) {}
    std::set<int> pdgs;
    std::set<int> allowedDaughters; ///< To not count eg Michel electrons against purity
    TH1D *purity, *efficiency, *deltaLength;
    TH2D *purityVsLength, *efficiencyVsLength, *deltaLengthVsLength;
    TH2D* effvpur;
  };

  void CheckRecoClusters(std::vector< art::Ptr<rb::Cluster> > const& clscol,
                         art::PtrVector<rb::CellHit>          const& allhits,
                         Plots& plots);

  /// Total length in either 2D or 3D
  double TotalLengthInDetector(const sim::Particle* part,
                               geo::View_t view) const;

  std::set<int> FindVisibleProngs(const art::PtrVector<rb::CellHit>& allhits,
                                  std::set<int> ids) const;

  std::string fClusterAssnLabel;   ///< label for module clusters, showers, and tracks are associated with
  std::string fShowerModuleLabel;  ///< label for module making the showers
  std::string fTrackModuleLabel;   ///< label for module making the tracks

  bool        fCheckShowers;       ///< should we check the reconstruction of showers?
  bool        fCheckTracks;        ///< should we check the reconstruction of tracks?
  bool        fLinkByMaxEff;       ///< If set to true,  finds mc particle with max efficiency 
                                   ///< and plots the corresponding efficiency and purity. If 
                                   ///< set to false, finds mc particle with the max purity and 
                                   ///< plots the corresponding efficiency and purity.

  std::vector<Plots> fPlots;
  Plots fClusterPlots;
  Plots fShowerPlots;
};


//-------------------------------------------------------------------
cheat::RecoCheckAna::RecoCheckAna(fhicl::ParameterSet const &p)
  : EDAnalyzer(p)
{
  this->reconfigure(p);
}

//-------------------------------------------------------------------
cheat::RecoCheckAna::~RecoCheckAna() 
{
}

//-------------------------------------------------------------------
// We want to determine how well an event finding algorithm works by 
// determining the following:
// 1. How much of the actual energy for a given interaction is 
//    found in the reconstructed event
// 2. How many of the particles from a given interaction are assigned 
//    to the correct reconstructed event
// 3. How close is the reconstructed event primary vertex to the 
//    simulated interaction primary vertex
// 4. What fraction of energy in the reconstructed event is from 
//    particles not associated with the simulated interaction

void cheat::RecoCheckAna::analyze(art::Event const &e) 
{
  // check that this is MC, stop if it isn't
  if(e.isRealData()){
    mf::LogWarning("RecoCheckAna") << "attempting to run MC truth check on "
                                   << "real data, bail";
    return;
  }


  // get a handle on the clusters that the clusters/tracks/showers are associated with (slices)
  art::Handle<std::vector<rb::Cluster> > hslices;
  e.getByLabel(fClusterAssnLabel, hslices);
  if(hslices.failedToGet()){
    mf::LogWarning("RecoCheckAna") << "could not recover slices";
  }


  // get all the association finder objects
  art::FindManyP<rb::Track>   fmtrack(hslices,e,fTrackModuleLabel);
  art::FindManyP<rb::Prong>   fmprong(hslices,e,fShowerModuleLabel);

  // loop over the slices
  for(unsigned int islc = 0; islc < hslices->size(); ++islc){
    art::PtrVector<rb::CellHit> allhits = hslices->at(islc).AllCells();    
    if(fCheckTracks){
      if(fmtrack.isValid()){
        std::vector< art::Ptr<rb::Track> > trkcol = fmtrack.at(islc);
        std::vector< art::Ptr<rb::Cluster> > tracks;
        for(unsigned int itrk = 0; itrk < trkcol.size(); ++itrk){
          tracks.push_back(trkcol[itrk]);
        }
        for(unsigned int n = 0; n < fPlots.size(); ++n){
          CheckRecoClusters(tracks, allhits, fPlots[n]);
        }
      }
      else{
        mf::LogWarning("RecoCheckAna") << "could not recover tracks";
      }
    }
    
    if(fCheckShowers){
      if(fmprong.isValid()){
        std::vector< art::Ptr<rb::Prong> > shwcol = fmprong.at(islc);
        std::vector< art::Ptr<rb::Cluster> > shws;
        for(unsigned int ishw = 0; ishw < shwcol.size(); ++ishw){
          shws.push_back(shwcol[ishw]);
        }
        this->CheckRecoClusters(shws, allhits, fShowerPlots);
      }
      else{
        mf::LogWarning("RecoCheckAna") << "could not recover showers";
      }
    }

  }

}

//-------------------------------------------------------------------
void cheat::RecoCheckAna::reconfigure(fhicl::ParameterSet const & p) 
{
  fClusterAssnLabel   = p.get< std::string >("ClusterAssnLabel"  );
  fShowerModuleLabel  = p.get< std::string >("ShowerModuleLabel" );
  fTrackModuleLabel   = p.get< std::string >("TrackModuleLabel"  );

  fCheckShowers       = p.get< bool        >("CheckShowers" );
  fCheckTracks        = p.get< bool        >("CheckTracks"  );
  fLinkByMaxEff       = p.get< bool        >("LinkByMaxEff" );
}

//-------------------------------------------------------------------
void cheat::RecoCheckAna::beginJob() 
{
  art::ServiceHandle<art::TFileService> tfs;

  if(fCheckTracks){
    for(int n = 0; n < 4; ++n){
      std::string name;
      std::set<int> pdgs;
      std::set<int> allowedDaughters;

      switch(n){
      case 0:
        name = "all";
        break;
      case 1:
        name = "muon";
        pdgs.insert(-13); pdgs.insert(+13);
        allowedDaughters.insert(-11); allowedDaughters.insert(+11);
        break;
      case 2:
        name = "proton";
        pdgs.insert(-2212); pdgs.insert(+2212);
        break;
      case 3:
        name = "piPlusMinus";
        pdgs.insert(-211); pdgs.insert(+211);
        break;
      }

      Plots p;
      p.purity = tfs->make<TH1D>((name+"Purity").c_str(), ";Purity;Tracks", 110, 0, 1.1);
      p.efficiency = tfs->make<TH1D>((name+"Efficiency").c_str(), ";Efficiency;Tracks", 110, 0, 1.1);
      p.deltaLength = tfs->make<TH1D>((name+"DeltaLength").c_str(), ";#DeltaL (cm);Tracks", 100, -50, +50);
      p.purityVsLength = tfs->make<TH2D>((name+"PurityVsLength").c_str(), ";True length (cm);Purity;Tracks", 300, 0, 1500, 110, 0, 1.1);
      p.efficiencyVsLength = tfs->make<TH2D>((name+"EfficiencyVsLength").c_str(), ";True length (cm);Efficiency;Tracks", 300, 0, 1500, 110, 0, 1.1);
      p.deltaLengthVsLength = tfs->make<TH2D>((name+"DeltaLengthVsLength").c_str(), ";True length (cm);#DeltaL (cm);Tracks", 300, 0, 1500, 100, -50, +50);
      p.effvpur = tfs->make<TH2D>((name+"EffvPurity").c_str(), ";Efficiency;Purity", 110, 0, 1.1, 110, 0, 1.1);
      p.pdgs = pdgs;
      p.allowedDaughters = allowedDaughters;

      fPlots.push_back(p);
    }
  }

  if(fCheckShowers){
    fShowerPlots.purity      = tfs->make<TH1D>("showerPurity",      ";Purity;Showers",          110, 0., 1.1);
    fShowerPlots.efficiency  = tfs->make<TH1D>("showerEfficiency",  ";Efficiency;Showers",      110, 0., 1.1);
  }

}


//-------------------------------------------------------------------
void cheat::RecoCheckAna::CheckRecoClusters(std::vector< art::Ptr<rb::Cluster> > const& clscol,
                                            art::PtrVector<rb::CellHit>          const& allhits,
                                            Plots& plots)
{
  // grab the set of eve IDs for this plist
  art::ServiceHandle<cheat::BackTracker> bt;
  std::set<int> trackIDs = bt->GetTrackIDList();

  // Only keep the ones associated with not completely unreasonable tracks
  trackIDs = FindVisibleProngs(allhits, trackIDs);

  // Empty set means accept all PDGs. Otherwise, filter out the unwanted ones
  if(!plots.pdgs.empty()){
    std::set<int> newTrackIDs;
    for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
      const int pdg = bt->TrackIDToParticle(*it)->PdgCode();
      if(plots.pdgs.find(pdg) != plots.pdgs.end()) newTrackIDs.insert(*it);
    } // end for it
    trackIDs = newTrackIDs;
  }

  std::map<int, int> parents;
  if(!plots.allowedDaughters.empty()){
    for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
      const sim::Particle* mother = bt->TrackIDToParticle(*it);
      if(!mother) continue; // Can this happen?

      for(int i = 0; i < mother->NumberDaughters(); ++i){
        const int id = mother->Daughter(i);
        assert(bt->TrackIDToParticle(id));
        if(plots.allowedDaughters.count(bt->TrackIDToParticle(id)->PdgCode())){
          parents[id] = *it;
        }
      } // end for i
    } // end for it
  } // end if

  // Indexed by what view they're in. X, Y, or 3D
  assert(geo::kXorY == 2); // Ensure enum still works how we think
  std::map<int, double> bestPur[3], bestEff[3], bestLength[3];

  // Go through each reconstructed object, accumulating the best scores seen
  // for each trackID we're interested in.
  for(size_t c = 0; c < clscol.size(); ++c){
    const rb::Cluster* cls = clscol[c].get();
    const geo::View_t view = cls->View();

    // get the hits associated with this object
    const art::PtrVector<rb::CellHit> hits = cls->AllCells();

    // use the cheat::BackTracker to find purity and efficiency for these
    // hits for each particle ID
    std::map<int, double> purMap, effMap;
    bt->HitCollectionPurity(trackIDs, hits, &purMap, &parents);
    bt->HitCollectionEfficiency(trackIDs, hits, allhits, view, &effMap);

    for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
      const double p = purMap[*it];
      const double e = effMap[*it];

      if(( fLinkByMaxEff && e > bestEff[view][*it]) ||
         (!fLinkByMaxEff && p > bestPur[view][*it])) {
        bestEff[view][*it] = e;
        bestPur[view][*it] = p;

        // If we were told to fill a track length histogram then downcast
        // should be safe.
        if(plots.deltaLength) bestLength[view][*it] = ((rb::Track*)cls)->TotalLength();
      }
    } // end for it
  } // end for c

  // For every track. If it was matched in 3D, fill the histograms.
  // If there is no 3D match look for 2 2D matches and fill with weight 0.5
  // If required we can add a flag to only try 3D or only try 2D.
  for(std::set<int>::iterator it = trackIDs.begin(); it != trackIDs.end(); ++it){
    // Try 3D first
    for(int view = geo::kXorY; view >= geo::kX; --view){
      // If no 3D match, continue to two 2D options
      if(view == geo::kXorY && bestEff[view].find(*it) == bestEff[view].end()) continue;

      const double weight = (view == geo::kXorY) ? 1 : 0.5;

      if(plots.purity) plots.purity->Fill(bestPur[view][*it], weight);
      if(plots.efficiency) plots.efficiency->Fill(bestEff[view][*it], weight);
      if(plots.effvpur) plots.effvpur->Fill(bestEff[view][*it], bestPur[view][*it], weight);

      const double trueLength = TotalLengthInDetector(bt->TrackIDToParticle(*it), geo::View_t(view));

      if(plots.purityVsLength) plots.purityVsLength->Fill(trueLength, bestPur[view][*it], weight);
      if(plots.efficiencyVsLength) plots.efficiencyVsLength->Fill(trueLength, bestEff[view][*it], weight);

      if(plots.deltaLength) plots.deltaLength->Fill(bestLength[view][*it]-trueLength, weight);
      if(plots.deltaLengthVsLength) plots.deltaLengthVsLength->Fill(trueLength, bestLength[view][*it]-trueLength, weight);

      // If we did find a 3D match, don't try 2D
      if(bestEff[geo::kXorY].find(*it) != bestEff[geo::kXorY].end()) break;
    } // end for view
  } // end for it
}

//-------------------------------------------------------------------
double cheat::RecoCheckAna::TotalLengthInDetector(const sim::Particle* p,
                                                  geo::View_t view) const
{
  const unsigned int N = p->NumberTrajectoryPoints();

  if(N < 2){
    mf::LogWarning("RecoCheckAna") << "Particle has no trajectory points";
    return 0;
  }

  art::ServiceHandle<geo::Geometry> geom;

  double dist = 0;

  for(unsigned int n = 0; n < N-1; ++n){
    TVector3 p0 = p->Position(n).Vect();
    TVector3 p1 = p->Position(n+1).Vect();
    if(p0.Z() < 0 || p1.Z() < 0 ||
       p0.Z() > geom->DetLength() || p1.Z() > geom->DetLength() ||
       fabs(p0.X()) > geom->DetHalfWidth() || fabs(p1.X()) > geom->DetHalfWidth() ||
       fabs(p0.Y()) > geom->DetHalfHeight() || fabs(p1.Y()) > geom->DetHalfHeight()) continue;

    if(view == geo::kX){
      p0.SetY(0);
      p1.SetY(0);
    }
    if(view == geo::kY){
      p0.SetX(0);
      p1.SetX(0);
    }
    dist += (p1-p0).Mag();
  }

  return dist;
}

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
};

//-------------------------------------------------------------------
std::set<int> cheat::RecoCheckAna::
FindVisibleProngs(const art::PtrVector<rb::CellHit>& allhits,
                  std::set<int> ids) const
{
  std::map<int, Counts> counts;

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

  // loop over all hits and check if this track id contributed to them
  for(unsigned int i = 0; i < allhits.size(); ++i){
    //get vector of track id's that contribute to this hit
    const art::Ptr<rb::CellHit> hit = allhits[i];
    const std::vector<TrackIDE> trackIDs = bt->HitToTrackIDE(hit);
    // loop over all the returned track ids
    int nTracks = 0; // number of tracks that contributed energy to the hit
    int loneID = -1; // Which ID it was if there was only one
    for(unsigned int j = 0; j < trackIDs.size(); ++j){
      const int trackID = trackIDs[j].trackID;
      if(trackID != sim::kNoiseId){
        ++nTracks;
        loneID = trackID;
      }

      // check if the current id matches the input track id
      if(ids.find(trackID) != ids.end()){
        ++counts[trackID].nHits[hit->View()];
      }
    } // end for j
    if(nTracks == 1){
      ++counts[loneID].nOnlyHits[hit->View()];
    }
  } // end for i

  const int kMinHits = 2; // minimum number of hits to be considered visible
  const int kMinOnlyHit = 1; // minimum number of hits in which it is only particle contributing energy to be visible

  std::set<int> ret;
  for(std::set<int>::iterator it = ids.begin(); it != ids.end(); ++it){
    const Counts c = counts[*it];
    // Check if this track passes the visible prong requirements
    if(c.nOnlyHits[geo::kX] >= kMinOnlyHit &&
       c.nOnlyHits[geo::kY] >= kMinOnlyHit){

      if(c.nHits[geo::kX] >= kMinHits && c.nHits[geo::kY] >= 1)
        ret.insert(*it);

      if(c.nHits[geo::kY] >= kMinHits && c.nHits[geo::kX] >= 1)
        ret.insert(*it);
    }
  }
  return ret;
}

DEFINE_ART_MODULE(cheat::RecoCheckAna)