ChannelPlots_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   TODO
/// \author  shanahan@fnal.gov
////////////////////////////////////////////////////////////////////////
#include <iostream>
#include <map>
#include <vector>

// Root
#include "TH1F.h"
#include "TH2F.h"

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Principal/Handle.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

// NOvA includes
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "GeometryObjects/PlaneGeo.h"
#include "RecoBase/Track.h"
#include "RecoBase/CellHit.h"
#include "DetectorPerformance/RawFebPlots.h"
#include "CMap/service/CMapService.h"


namespace dprf {

  class RawFebPlots;
  
  class ChannelPlots : public art::EDAnalyzer
  {
  public:
    explicit ChannelPlots(fhicl::ParameterSet const &p);
    ~ChannelPlots();
    
    void beginRun(const art::Run& run);
    virtual void analyze(art::Event const&);
    void endJob();
    
  private:
    
    // configurables
    const fhicl::ParameterSet &fPSet;
    const std::vector< unsigned int > fHumanReadableFEBs;
    const unsigned int fMinPlaneTrack,fMinCellTrack,fMinCellTrackView; //track cuts
    const std::string fRawDigitLabel,fCellHitLabel,fTrackLabel; // module names
    typedef enum _cutlevel { kAll, kPass, kNumLevel} Cuts_t;
    
    struct ViewH {
      TH2F* hists[2];
    };
    std::map<unsigned int, ViewH> fHBookedChanCP;// map of booked channels in cell-plane, vs group
    TH2F* fHRawDigits[geo::kXorY][kNumLevel]; // plots of raw digits
    TH1F* fHExtentPlaneTrk[kNumLevel]; // plots of raw digits
    TH1F* fHNCellTrk[kNumLevel]; // plots of raw digits
    TH1F* fHNCellViewTrk[geo::kXorY][kNumLevel]; // plots of raw digits
    
    typedef std::map< unsigned int, RawFebPlots* > FebPlotGroup;
    typedef std::map< unsigned int, FebPlotGroup > FebPlotSet;
    
    FebPlotSet fFebRawAll;
    FebPlotSet fFebRawOnTrack;
    
    
    void FillFebRawPlots(FebPlotSet&, art::Ptr<rawdata::RawDigit> );
    bool TrackOK(art::Ptr<rb::Track> );
    std::string HRFebToStringLabel(unsigned int, 
                                   int, 
                                   unsigned int&, unsigned int&, unsigned int&, unsigned int&);
    
    void ResetEventPlots();
    void FillEventPlots();
    
  };
  
  
}//namespace

//////////////////////////////////////////////////////////////////////
namespace dprf {

  ChannelPlots::ChannelPlots(fhicl::ParameterSet const &pset):
  EDAnalyzer(pset),
  fPSet(pset),
  fHumanReadableFEBs(pset.get< std::vector< unsigned int > >("FEBs")),
  fMinPlaneTrack(pset.get< unsigned int >("MinPlaneTrack")),
  fMinCellTrack(pset.get< unsigned int >("MinCellTrack")),
  fMinCellTrackView(pset.get< unsigned int >("MinCellTrackView")),
  fRawDigitLabel(pset.get< std::string >("RawDigitLabel")),
  fCellHitLabel(pset.get< std::string >("CellHitLabel")),
  fTrackLabel(pset.get< std::string >("TrackLabel"))
  {
  }
  
  ChannelPlots::~ChannelPlots()
  {
  }

  void ChannelPlots::beginRun(const art::Run& run)
  {

    //    std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ geo::kXorY = " << "Here1" << std::endl;

    //   if(fHBookedChanCP.empty()) return; // Don't create them all twice
    
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle <art::TFileService> tfs;
    art::ServiceHandle<cmap::CMap>    cmap;


    
    for (unsigned int ichan=0; ichan<fHumanReadableFEBs.size(); ichan++) {
      const int det=geom->DetId();
      
      //      std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ geo::kXorY = " << "Here1" << std::endl;


      unsigned int group,diblock,ldcm,lfeb;
      std::string label=HRFebToStringLabel(fHumanReadableFEBs[ichan],det,
                                           group,diblock,ldcm,lfeb);
      const daqchannelmap::dchan chan = cmap->Map()->encodeDChan(det,diblock,ldcm,lfeb,0);
      
      bool detailed=true; //need to put in config...
      
      (fFebRawAll[group])[chan]=
        new RawFebPlots(label,"All",detailed,fPSet);
      (fFebRawOnTrack[group])[chan]=
        new RawFebPlots(label,"OnTracks",detailed,fPSet);
      
      //book and fill histogram showing which channels were booked
      if (fHBookedChanCP.find(group) == fHBookedChanCP.end())
        {
          ViewH hist={{0,0}};
          for (unsigned int iv=geo::kX; iv<geo::kXorY; iv++) 
            {

              geo::View_t view=(geo::View_t)iv;
              std::string viewName(view == geo::kX ? "X" : "Y");
              const std::set<unsigned int>& planes=geom->GetPlanesByView(view);
              const geo::PlaneGeo* plane=(geom->Plane(*(planes.begin())));
              //(fHBookedChanCP[group])[iv]=
              hist.hists[iv]=
                tfs->make<TH2F>(Form("HBookedChanCP_%s_%d",viewName.c_str(),group),
                                Form("Booked Channel, %s, group %d",viewName.c_str(),group),
                                geom->NPlanes()+1,-.5,geom->NPlanes()+.5,
                                plane->Ncells()+1,-.5,plane->Ncells()+.5);
            }
          fHBookedChanCP.insert(std::pair<unsigned int, ViewH>
                                (group,hist));
        }
      
      
          for (unsigned int ipix=0; 
               ipix<daqchannelmap::NUM_PIXELS; ipix++) {
            const unsigned int chan            = cmap->Map()->encodeDChan(det,diblock,ldcm,lfeb,ipix);
            const daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
            const unsigned int plannum         = cmap->Map()->getPlane(logchan);
            const geo::PlaneGeo* plane         = geom->Plane(plannum);
            const unsigned int view=(unsigned int)plane->View();
            
            (fHBookedChanCP[group]).hists[view]
              ->Fill(cmap->Map()->getPlane(logchan),
                     cmap->Map()->getCell(logchan));
          }//pixels
          
    }//booked channels
    
    // book general hit histograms
    std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ geo::kXorY = " << "Here" << std::endl;

    for (unsigned int iv=geo::kX; iv<geo::kXorY; iv++) 
      {
        std::cout<< "$$$$$$$$$$$$$$$$$$$$$$$$$$$$$$ LOOPING  $$$$$$$$ "<< std::endl;

        geo::View_t view=(geo::View_t)iv;
        std::string viewName(view == geo::kX ? "X" : "Y");
        const std::set<unsigned int>& planes=geom->GetPlanesByView(view);
        const geo::PlaneGeo* plane=(geom->Plane(*(planes.begin())));
        
        for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
          //for (Cuts_t cut=kAll; cut<kNumLevel; cut++) {
          std::string cutname(cut==kAll ? "All" : "On Good Tracks");
          fHRawDigits[iv][cut]=
            tfs->make<TH2F>(Form("HRawDigits_%s_%d",viewName.c_str(),cut),
                            Form("Raw Digits, %s, %s",viewName.c_str(),cutname.c_str()),
                            geom->NPlanes()+1,-.5,geom->NPlanes()+.5,
                            plane->Ncells()+1,-.5,plane->Ncells()+.5);
          
        }//Cut levels
        
      }//view 
    
    // book track cut histograms
    for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
      std::string cutname(cut==kAll ? "All Tracks" : "Good Tracks");
      fHExtentPlaneTrk[cut]=
        tfs->make<TH1F>(Form("HExtentPlaneTrk_%d",cut),
                        Form("Track Extent, Planes, %s",cutname.c_str()),
                        geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
      fHNCellTrk[cut]=
        tfs->make<TH1F>(Form("HNCellTrk_%d",cut),
                        Form("Track Number of Cells, %s",cutname.c_str()),
                        geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
      fHNCellViewTrk[geo::kX][cut]=
        tfs->make<TH1F>(Form("HNCellViewTrk_X_%d",cut),
                        Form("Track Number of Cells, X, %s",cutname.c_str()),
                        geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
      fHNCellViewTrk[geo::kY][cut]=
        tfs->make<TH1F>(Form("HNCellViewTrk_Y_%d",cut),
                        Form("Track Number of Cells, Y, %s",cutname.c_str()),
                        geom->NPlanes()+1,-.5,geom->NPlanes()+.5);
    }// track cuts
    
  }
  void ChannelPlots::endJob()
  {
  }
  
  
  std::string ChannelPlots::HRFebToStringLabel(
                                               unsigned int hrfeb, int det, 
                                               unsigned int &group, unsigned int &diblock, unsigned int &ldcm,
                                               unsigned int &lfeb)
  {
    // return a useful label from a human-readable FEB id.  And return
    // other useful info
    group= (hrfeb/1000000)%100;
    diblock=(hrfeb/  10000)%100;
    ldcm=   (hrfeb/    100)%100;
    lfeb=   (hrfeb/      1)%100;
    
    std::string label(Form("g%d_%d_%2.2d_%2.2d_%2.2d",group,det,diblock,ldcm,lfeb));
    return label;
    
  }  
  
  void ChannelPlots::analyze(art::Event const& evt)
  {
    
    ResetEventPlots();
   
    // look for RawDigits.  Don't continue if none
    art::Handle< std::vector<rawdata::RawDigit> > digs;
    try {
      evt.getByLabel(fRawDigitLabel,digs);
    }
    catch (...) {
      mf::LogWarning("CPLTWarn") << "Failed to find a RawDigit, bail";
      return ;
    }
    
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<cmap::CMap>    cmap;
    
    // make "all hits" plots
    for (unsigned int id=0; id<digs->size(); id++)
      {

        art::Ptr<rawdata::RawDigit> dig(digs,id);
        unsigned int chan=dig->DaqChannel();
        
        // fill general plots
        // need to use "offline" view, not online....
        daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
        const unsigned int plannum   = cmap->Map()->getPlane(logchan);
        const geo::PlaneGeo* plane   = geom->Plane(plannum);
        const unsigned int view      = (unsigned int)plane->View();
        
        
        fHRawDigits[view][kAll]->Fill(plannum, cmap->Map()->getCell(logchan));
        
        FillFebRawPlots(fFebRawAll,dig);

        
      }// digit id
    
    // look for Tracks. OK to return w/out error if none
    art::Handle< std::vector<rb::Track> > tracks;
    try {
      evt.getByLabel(fTrackLabel,tracks);
    }
    catch (...) {
      FillEventPlots();
      return ;
    }
    
    // can it exist but have zero size?
    if (!tracks->size()) {
      FillEventPlots();
      return;
    }
    
    // get the CellHits, so we can find RawHits on tracks
    art::Handle< std::vector<rb::CellHit> > cells;
    try {
      evt.getByLabel(fCellHitLabel,cells);
    }
    catch (...) {
      mf::LogError("CPLTWarn") << "\nFound Tracks but no CellHits!\n";
      assert(0);
    }
    
    for(unsigned int tr = 0; tr < tracks->size(); tr++)
      {
        art::Ptr<rb::Track> track(tracks,tr);
        
        // does the track have a decent number of hits and planes in each view?
        
        if (TrackOK(track)) {
          // loop over CellHits on tracks
          for (unsigned int ic=0; ic<track->NCell(); ic++) 
            {
              art::Ptr<rb::CellHit> cell   = track->Cell(ic);
              const unsigned int chan      = cell->DaqChannel();
              daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
              const unsigned int plannum   = cmap->Map()->getPlane(logchan);
              const geo::PlaneGeo* plane   = geom->Plane(plannum);
              unsigned int view=(unsigned int)plane->View();
              
              fHRawDigits[view][kPass]->Fill(cmap->Map()->getPlane(logchan),
                                             cmap->Map()->getCell(logchan));
              
              FillFebRawPlots(fFebRawOnTrack,cell);
            }//cell
        } // track is OK
        
      } // tracks
    
    FillEventPlots();
    
  } //ana
  
  void ChannelPlots::FillFebRawPlots(FebPlotSet& plots, 
                                     art::Ptr<rawdata::RawDigit> dig)
  {
    
    //Fill all FebRawPlots for this channel, in all groups
    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<cmap::CMap>    cmap;
    
    unsigned int chan=dig->DaqChannel();
    
    //key for the FEB is pixel 0 dchan for that FEB.
    unsigned int lookupchan=chan & 
      ~(daqchannelmap::MASK_PIXELID);
    
    //loop over groups, and fill histograms if this channel is booked
    for (FebPlotSet::iterator iter=plots.begin(); iter!=plots.end(); iter++) 
      {
        
        //unsigned int group=iter->first;
        if (iter->second.find(lookupchan) != iter->second.end()) {
          iter->second[lookupchan]->AddHit(dig->ADC(),
                                           cmap->Map()->getPixel(chan));
        }
      } // groups
  } // FillFebRawPlots

  bool ChannelPlots::TrackOK(art::Ptr<rb::Track> track)
  {
    //some histos to monitor the cuts
    fHExtentPlaneTrk[kAll]->Fill(track->ExtentPlane());
    fHNCellTrk[kAll]->Fill(track->NCell());
    fHNCellViewTrk[geo::kX][kAll]->Fill(track->NXCell());
    fHNCellViewTrk[geo::kY][kAll]->Fill(track->NYCell());
    
    // enough planes (extent)
    if (track->ExtentPlane() < fMinPlaneTrack)     return false;
    
    // enough cells
    if (track->NCell() < fMinCellTrack)            return false;
    
    // enough cells in each view
    if (track->NXCell() < fMinCellTrackView ||
        track->NYCell() < fMinCellTrackView )       return false;
    
    fHExtentPlaneTrk[kPass]->Fill(track->ExtentPlane());
    fHNCellTrk[kPass]->Fill(track->NCell());
    fHNCellViewTrk[geo::kX][kPass]->Fill(track->NXCell());
    fHNCellViewTrk[geo::kY][kPass]->Fill(track->NYCell());
    return true;
  }
  void ChannelPlots::FillEventPlots()
  {
    // Fill the "per-event" plots
    FebPlotSet::iterator siter=fFebRawAll.begin();
    for (;siter!=fFebRawAll.end();siter++) {
      FebPlotGroup::iterator giter=(siter->second).begin();
      for (;giter!=(siter->second).end();giter++) 
        (giter->second)->FinishEvent();
    }
    
    siter=fFebRawOnTrack.begin();
    for (;siter!=fFebRawOnTrack.end();siter++) {
      FebPlotGroup::iterator giter=(siter->second).begin();
      for (;giter!=(siter->second).end();giter++) 
        (giter->second)->FinishEvent();
    }
  }
  void ChannelPlots::ResetEventPlots()
  {
    // Reset the "per-event" counters
    FebPlotSet::iterator siter=fFebRawAll.begin();
    for (;siter!=fFebRawAll.end();siter++) {
      FebPlotGroup::iterator giter=(siter->second).begin();
      for (;giter!=(siter->second).end();giter++) 
        (giter->second)->Reset();
    }
    
    siter=fFebRawOnTrack.begin();
    for (;siter!=fFebRawOnTrack.end();siter++) {
      FebPlotGroup::iterator giter=(siter->second).begin();
      for (;giter!=(siter->second).end();giter++) 
        (giter->second)->Reset();
    }
  }
} //namespace dprf

///////////////////////////////////////////////////////////////
namespace dprf
{
  DEFINE_ART_MODULE(ChannelPlots)
}