TrackPlots_module.cc

Go to the documentation of this file.

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

// Root
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.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 "CMap/service/CMapService.h"
#include "DAQChannelMap/DAQChannelMap.h"
#include "DetectorPerformance/RawFebPlots.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/CellGeo.h"
#include "NovaDAQConventions/DAQConventions.h"
#include "RecoBase/Track.h"
#include "RecoBase/CellHit.h"


namespace dprf {
  
  class TrackPlots : public art::EDAnalyzer
  {
  public:
    explicit TrackPlots(fhicl::ParameterSet const &p);
    ~TrackPlots();
    
    void beginRun(const art::Run& run);
    virtual void analyze(art::Event const&);
    void endJob();
    
  private:
    
    unsigned int DcmKey(unsigned int dib,unsigned int ldc) {
      return dib*100+ldc; };
    unsigned int DiBlock(unsigned int key) {
      return key/100; };
    unsigned int LDCM(unsigned int key) {
      return key%100; };
    
    // configurables
    //const fhicl::ParameterSet &fPSet;
    const unsigned int fMinPlaneTrack,fMinCellTrack,fMinCellTrackView; ///< track cuts
    const unsigned int fMinTrackHitDCM; ///< minimum number of hits on a track in a DCM to use it for comparision to other DCMs
    const unsigned int fEventGroup; ///< Number of events per group
    const std::string fRawDigitLabel,fCellHitLabel,fTrackLabel; ///< module names
    unsigned int fCallNumber;///< keep track number of calls, for defining event groups
    typedef enum _groups { kAllG, kEvent, kNumGroup} Groups_t;
    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
    TH1F* fHTrkCellXMinusY[kNumLevel];           ///< difference in cells per view
    
    TH2F* fHTracksDCM[kNumGroup][kNumLevel];           ///< DCMs with hits on tracks
    TH2F* fHTracksDCMHits[kNumGroup][kNumLevel];       ///< Number of hits per track on DCM
    TH2F* fHTracksOtherDCMs[kNumGroup][kNumLevel];     ///< Number of other DCMs with hits on track
    TH2F* fHTracksOtherDCMHits[kNumGroup][kNumLevel];  ///< Number of hits on other DCMs with hits on track
    TH2F* fHTracksOtherDCMRatio[kNumGroup][kNumLevel]; ///< Ratio of hits on track on DCM to hits on track on all other DCMs
    TH2F* fHTrackEnd[geo::kXorY][kNumLevel];           ///< (Cell,Plane) for endpoints
    // of track in Z, in each view
    
    std::map<unsigned int, std::map<unsigned int, float> > fTracksOtherDCMRatio; // map of ratio vs. event group (first key), for each dcm (second key)
    bool TrackOK(art::Ptr<rb::Track> );
    void ResetEventGroup();
    void FinishEventGroup(unsigned int);
    void Normalize(unsigned int);
    
  };

  
}//namespace


namespace dprf {
  
  TrackPlots::TrackPlots(fhicl::ParameterSet const &pset):
    EDAnalyzer(pset),
    //fPSet(pset),
    fMinPlaneTrack(pset.get< unsigned int >("MinPlaneTrack")),
    fMinCellTrack(pset.get< unsigned int >("MinCellTrack")),
    fMinCellTrackView(pset.get< unsigned int >("MinCellTrackView")),
    fMinTrackHitDCM(pset.get< unsigned int >("MinTrackHitDCM")),
    fEventGroup(pset.get< unsigned int >("EventGroup")),
    fRawDigitLabel(pset.get< std::string >("RawDigitLabel")),
    fCellHitLabel(pset.get< std::string >("CellHitLabel")),
    fTrackLabel(pset.get< std::string >("TrackLabel")),
    fCallNumber(0),
    fHExtentPlaneTrk{0, 0} // Flag the others are uninitialized too
  { 
  }

  TrackPlots::~TrackPlots()
  {
  }
  
  void TrackPlots::beginRun(const art::Run& run)
  {
    if(fHExtentPlaneTrk[kAll]) return; // Don't make histograms twice
    
    art::ServiceHandle <geo::Geometry> geom;
    art::ServiceHandle <art::TFileService> tfs;
    const geo::PlaneGeo* plane;
    
    unsigned int ndcm = 12;    //   # of DCM in Far Detector 
    unsigned int ndb  = 14;    //   # of diblocks in Far Detector
    
    if (geom->DetId() == novadaq::cnv::kNDOS){
      ndcm = 3; //   # of DCM in NDOS 
      ndb  = 4;
    }
    if (geom->DetId() == novadaq::cnv::kNEARDET){
      ndcm = 4; //   # of DCM in NEAR Detector 
      ndb  = 4;
    }
    
    
    for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
      std::string cutname(cut==kAll ? "All" : "On Good Tracks");
      
      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())));
          plane=(geom->Plane(*(planes.begin())));
          
          fHTrackEnd[iv][cut]=
            tfs->make<TH2F>(Form("HTrackEnd_%s_%d",viewName.c_str(),cut),
                            Form("Track Endpoints, %s, %s",viewName.c_str(),cutname.c_str()),
                            geom->NPlanes()+1,-.5,geom->NPlanes()+.5,
                            plane->Ncells()+1,-.5,plane->Ncells()+.5);
          
        }//view 
      
      fHExtentPlaneTrk[cut]=
        //
        // *********  Hardcoded for NDOS *************
        //
        //      tfs->make<TH1F>(Form("HTrackPlanes%d",cut),
        //              Form("Track Extent (Planes), %s",cutname.c_str()),
        //              200.,-.5,199.5);
        //
        tfs->make<TH1F>(Form("HTrackPlanes%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);
      
      fHTrkCellXMinusY[cut]=
        tfs->make<TH1F>(Form("HTrkCellXMinusY_%d",cut),
                        Form("Track Cell difference, X-Y view, %s",cutname.c_str()),
                        plane->Ncells()+1,-.5,plane->Ncells()+.5);
                        //                              100, -.5, 99.5);
      
      //igrp=0 for all events, igrp=1 event-wise (or group-wise)
      for (unsigned int igrp=0; igrp<2; igrp++) {
        std::string groupname(igrp==kAll ? "(All Events)" : "Event Group");
        fHTracksDCM[igrp][cut]=
          tfs->make<TH2F>(Form("HTracksDCM_%d_%d",igrp,cut),
                          Form("Tracks in DCM (ldcm vs. db), %s %s",groupname.c_str(),cutname.c_str()),
                          ndb,+.5,(float)ndb+.5,ndcm,.5,(float)ndcm+.5);
        fHTracksDCMHits[igrp][cut]=
          tfs->make<TH2F>(Form("HTracksDCMHits_%d_%d",igrp,cut),
                          Form("Track Hits in DCM (ldcm vs. db), %s %s",groupname.c_str(),cutname.c_str()),
                          ndb,+.5,(float)ndb+.5,ndcm,.5,(float)ndcm+.5);
        
        fHTracksOtherDCMs[igrp][cut]=
          tfs->make<TH2F>(Form("HTracksOtherDCMs_%d_%d",igrp,cut),
                          Form("Number of Other DCMs with Track Hits, vs. DCM (ldcm vs. db), %s %s",groupname.c_str(),cutname.c_str()),
                          ndb,+.5,(float)ndb+.5,ndcm,.5,(float)ndcm+.5);
        
        fHTracksOtherDCMHits[igrp][cut]=
          tfs->make<TH2F>(Form("HTracksOtherDCMHits_%d_%d",igrp,cut),
                          Form("Number of Track Hits on Other DCMs, vs. DCM (ldcm vs. db), %s %s",groupname.c_str(),cutname.c_str()),
                          ndb,+.5,(float)ndb+.5,ndcm,.5,(float)ndcm+.5);
        
        fHTracksOtherDCMRatio[igrp][cut]=
          tfs->make<TH2F>(Form("HTracksOtherDCMRatio_%d_%d",igrp,cut),
                          Form("Ratio of Track Hits on Other DCMs vs. DCM (ldcm vs. db), %s %s",groupname.c_str(),cutname.c_str()),
                          ndb,+.5,(float)ndb+.5,ndcm,.5,(float)ndcm+.5);
      } //igrp
      
      
    }//Cut levels
     
  } // end beginRun

  void TrackPlots::endJob()
  {
    art::ServiceHandle <art::TFileService> tfs;
    
    // normalize "per-track" plots
    Normalize(kAllG);
    
    // make plots vs. event group
    if (!fTracksOtherDCMRatio.size()) return;
    std::map<unsigned int, std::map<unsigned int, float> >::iterator
      iter=fTracksOtherDCMRatio.begin();
    std::map<unsigned int, std::map<unsigned int, float> >::iterator
      end=fTracksOtherDCMRatio.end();
    end--;
    
    float bgroup=(float) iter->first;
    float fgroup=(float) end->first;
    
    TProfile* hratioall=
      tfs->make<TProfile>("HRatioOtherDCM_All",
                          "Ratio of Hits on other DCMs vs. event, All DCMs",
                          fTracksOtherDCMRatio.size()+1,bgroup-.5,fgroup+.5,"S");
    
    std::map<unsigned int, TProfile*> eventPlots; // per-dcm

    // loop over entries - corresponding to event groups
    for (; iter!=fTracksOtherDCMRatio.end(); iter++) {
      
      float event=(float) iter->first;
      
      // loop over DCMs within each event group
      std::map<unsigned int, float> gmap=iter->second;
      std::map<unsigned int, float>::iterator siter=gmap.begin();
      for(; siter!=gmap.end(); siter++) {
        unsigned int key=siter->first;
        unsigned int dib=DiBlock(key);
        unsigned int ldc=LDCM(key);
        float ratio=siter->second;
        
        if (eventPlots.find(key) == eventPlots.end() || !eventPlots[key]) {
          eventPlots[key]=
            tfs->make<TProfile>(Form("HRatioOtherDCM_%2.2d_%2.2d",dib,ldc),
                                Form("Ratio of Hits on other DCMs vs. event, %2.2d-%2.2d",dib,ldc),
                                fTracksOtherDCMRatio.size()+1,bgroup-.5,fgroup+.5,"S");
        }
        eventPlots[key]->Fill(event,ratio);
        hratioall->Fill(event,ratio);
      }
    }
    
  } // end endRun
  
  void TrackPlots::analyze(art::Event const& evt)
  {
    
    ResetEventGroup();
    
    // 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;
    /* -- until needed
    // 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....
       unsigned int plannum=daqchannelmap::DAQChannelMap::getPlaneGlobal(chan);
       const geo::PlaneGeo* plane=geom->Plane(plannum);
       unsigned int view=(unsigned int)plane->View();
       
       
       

    }// 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 (...) {
      return ;
    }
    
    // can it exist but have zero size?
    if (!tracks->size()) {
      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++)
      {
        
        //count DCMs on the track
        std::map<unsigned int, unsigned int>  DCMsHit;
        
        art::Ptr<rb::Track> track(tracks,tr);
        
        bool track_ok=TrackOK(track);
        
        // keep track of start and end hit in each view
        unsigned short cellx[2][2]={{9999,0},{9999,0}};
        unsigned short planex[2][2]={{9999,0},{9999,0}};
        
        // 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();
            const daqchannelmap::lchan logchan = cmap->Map()->encodeLChan(chan);
            const unsigned int dib             = cmap->Map()->getDiBlock(chan);
            const unsigned int ldc             = cmap->Map()->getDCM(chan);
            DCMsHit[DcmKey(dib,ldc)]++;
            
            
            const unsigned int plannum = cmap->Map()->getPlane(logchan);
            if (plannum!=cell->Plane()) {
              std::cout<<"Why is plane from CellHit != plane from chan???"<<std::endl;
              assert(0);
            }
            
            unsigned int cellnum=(unsigned int) cell->Cell();
         
            const geo::PlaneGeo* plane=geom->Plane(plannum);
            unsigned int view=(unsigned int)plane->View();
            
            if (plannum<planex[view][0]) {
              cellx[view][0]=cellnum;
              planex[view][0]=plannum;
            }
            if (plannum>planex[view][1]) {
              cellx[view][1]=cellnum;
              planex[view][1]=plannum;
            }
            
            
          }//cell
        
        //count, for each DCM, hits on the track in all others
        std::map<unsigned int, unsigned int>::iterator iter=DCMsHit.begin();
        for (; iter!=DCMsHit.end(); iter++) {
          unsigned int dcmkey=iter->first;
          float otherdcmhits=0;
          float otherdcm=0;
          std::map<unsigned int, unsigned int>::iterator oter=DCMsHit.begin();
          for (; oter!=DCMsHit.end(); oter++) {
            unsigned int otherdcmkey=oter->first;
            if (otherdcmkey==dcmkey) continue;
            otherdcmhits+=(float)oter->second;
            otherdcm++;
          } //other DCMs
          
          float thisdcmhits=(float) iter->second;
          float ratioOtherToThis=otherdcmhits/thisdcmhits;
          if (thisdcmhits>fMinTrackHitDCM) { 
            for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
              // plot all tracks for cut=0, only good ones for cut=1
              if (cut==1&&!track_ok) continue;
              
              unsigned int dib=DiBlock(dcmkey);
              unsigned int ldc=LDCM(dcmkey);
              
              //igrp=0 for all events, igrp=1 for group of events
              for (unsigned int igrp=0; igrp<2; igrp++) {
                fHTracksDCM[igrp][cut]->Fill(dib,ldc);
                fHTracksDCMHits[igrp][cut]->Fill(dib,ldc,thisdcmhits);
                fHTracksOtherDCMs[igrp][cut]->Fill(dib,ldc,otherdcm);
                fHTracksOtherDCMHits[igrp][cut]->Fill(dib,ldc,otherdcmhits);
                fHTracksOtherDCMRatio[igrp][cut]->Fill(dib,ldc,ratioOtherToThis);
              }
              for (unsigned int iend=0; iend<2; iend++)
                for (unsigned int iv=geo::kX; iv<geo::kXorY; iv++)
                   fHTrackEnd[iv][cut]->Fill(planex[iv][iend],cellx[iv][iend]);
            } // cut
          }//min hits on this DCM
          
        } // "this" DCM
        
      } // tracks
    
    FinishEventGroup((unsigned int) evt.id().event());
    
  }//ana
  
  bool TrackPlots::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());
    
    // and some others
    fHTrkCellXMinusY[kAll]->Fill(track->NXCell()-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());
    fHTrkCellXMinusY[kPass]->Fill(track->NXCell()-track->NYCell());
    return true;
  }

  
  void TrackPlots::ResetEventGroup()
  {
    fCallNumber++;
    if (fCallNumber) return; // don't do anything if not a new group
    
    for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
      fHTracksDCM[kEvent][cut]->Reset();
      fHTracksDCMHits[kEvent][cut]->Reset();
      fHTracksOtherDCMs[kEvent][cut]->Reset();
      fHTracksOtherDCMHits[kEvent][cut]->Reset();
      fHTracksOtherDCMRatio[kEvent][cut]->Reset();
    }
  } 

  
  void TrackPlots::FinishEventGroup(unsigned int event)
  {
    
    if (fCallNumber<fEventGroup) return; // don't do anything until group 
    // is finished...
    fCallNumber=-1; //ensure ResetEventGroup at next call
    /*
      if (fTracksOtherDCMRatio.size()==0) {
      fTracksOtherDCMRatio=std::map<unsigned int, std::map<unsigned int, float> > ();
      }
      
    */
    
    
    // first, add the group histograms to the total
    for (unsigned int cut=kPass; cut<kNumLevel; cut++) {
      fHTracksDCM[kAll][cut]->Add(fHTracksDCM[kEvent][cut]);
      fHTracksDCMHits[kAll][cut]->Add(fHTracksDCMHits[kEvent][cut]);
      fHTracksOtherDCMs[kAll][cut]->Add(fHTracksOtherDCMs[kEvent][cut]);
      fHTracksOtherDCMHits[kAll][cut]->Add(fHTracksOtherDCMHits[kEvent][cut]);
      fHTracksOtherDCMRatio[kAll][cut]->Add(fHTracksOtherDCMRatio[kEvent][cut]);
    }
    
    // normalize the group-wise plots
    Normalize(kEvent);
    
    // fill map of ratio, per event group, per dcm
    TH1* hist=fHTracksOtherDCMRatio[kEvent][kPass];
    for (unsigned int dib=1; dib<=(unsigned int) hist->GetXaxis()->GetNbins(); dib++) 
      {
        for (unsigned int ldc=1; ldc<=(unsigned int) hist->GetYaxis()->GetNbins(); ldc++) 
          {
            // not sure this is needed... was having segv's, but not sure it was due to this.
            std::map<unsigned int, std::map<unsigned int, float> >::iterator find=
              fTracksOtherDCMRatio.find(event);
            std::map<unsigned int, std::map<unsigned int, float> >::iterator end=fTracksOtherDCMRatio.end();
            if (find==end) {
              std::map<unsigned int, float> amap=std::map<unsigned int, float> ();
              fTracksOtherDCMRatio[event]=amap;
            }
            fTracksOtherDCMRatio[event][DcmKey(dib,ldc)]=
              hist->GetBinContent(dib,ldc);
          }
      }
  }

  
  void TrackPlots::Normalize(unsigned int igrp)
  {
    for (unsigned int cut=kAll; cut<kNumLevel; cut++) {
      fHTracksDCMHits[igrp][cut]->Divide(fHTracksDCM[igrp][cut]);
      fHTracksOtherDCMs[igrp][cut]->Divide(fHTracksDCM[igrp][cut]);
      fHTracksOtherDCMHits[igrp][cut]->Divide(fHTracksDCM[igrp][cut]);
      fHTracksOtherDCMRatio[igrp][cut]->Divide(fHTracksDCM[igrp][cut]);
    }
  }
} //namespace dprf

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