ChanOcc_module.cc

Go to the documentation of this file.

///
/// \file    ChanOcc.cxx
/// \brief   Make ntuples of channel occupancies
/// \version $Id: ChanOcc.cxx,v 1.7 2012-09-25 16:06:22 gsdavies Exp $
/// \author  messier@indiana.edu
///
#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Core/ModuleMacros.h"

#include "CMap/service/CMapService.h"
#include "GeometryObjects/CellGeo.h"
#include "Geometry/Geometry.h"
#include "GeometryObjects/PlaneGeo.h"
#include "RawData/RawDigit.h"

#include "TNtuple.h"

class TNtuple;

namespace comi {
  /// Module to produce ntuples  of channel occupancies
  class ChanOcc : public art::EDAnalyzer {
  public:
    explicit ChanOcc(fhicl::ParameterSet const& pset);
    virtual ~ChanOcc();
    
    void beginJob();
    void endSubRun(art::SubRun const& sr);
    
    void analyze(art::Event const& evt);
    void reconfigure(const fhicl::ParameterSet& p);

  private:
    static const int kPLANE_MAX = 2048;
    static const int kCELL_MAX  = 512;
    double fNhit[kPLANE_MAX][kCELL_MAX]; ///< Number of hits by plane/cell
    double fNevt;                        ///< Number of events seen
    
    TNtuple* fNt; ///< Ntuple summarizing channel occupancies
  private:
    float fADCThresh; ///< Filter hits below this ADC threshold
  };

  //......................................................................
  
  ChanOcc::ChanOcc(fhicl::ParameterSet const& p)
    : EDAnalyzer(p)
  {
    memset(fNhit, 0, kPLANE_MAX*kCELL_MAX*sizeof(float));
    this->reconfigure(p);
  }
  
  //......................................................................
  
  void ChanOcc::beginJob() 
  {
    art::ServiceHandle<art::TFileService> f;
    fNt = f->make<TNtuple>("occ",
                           "Channel Occupancies",
                           "run:subrun:view:ip:ic:occ");
  }
  
  //......................................................................
  
  ChanOcc::~ChanOcc() { }
  
  //......................................................................
  
  void ChanOcc::analyze(art::Event const& evt)
  {
    unsigned int i;
    
    fNevt += 1.0;
    
    //
    // Pull out the raw digits from the event
    //
    art::Handle< std::vector<rawdata::RawDigit> > digidummy;
    evt.getByLabel("daq",digidummy);
    
    //
    // Shuffle the data into a vector we can manipulate
    //
    std::vector< art::Ptr<rawdata::RawDigit> > digi;
    for (i=0; i<digidummy->size(); ++i) {
      digi.push_back(art::Ptr<rawdata::RawDigit>(digidummy, i));
    }
    
    //
    // Record the number of hits
    //
    art::ServiceHandle<cmap::CMap> cmap;
    for (i=0; i<digi.size(); ++i) {
      if (digi[i].get()->ADC()<fADCThresh) continue;
      
      int p = cmap->GetPlane(digi[i].get());
      int c = cmap->GetCell(digi[i].get());
      
      if (p>=0 && p<kPLANE_MAX &&	c>=0 && c<kCELL_MAX) {
        fNhit[p][c] += 1.0;
      }
      else abort();
    }
  }
  
  //......................................................................
  
  void ChanOcc::endSubRun(art::SubRun const& sr) 
  {
    if (fNevt<=10000) return;
    
    int run    = sr.run();
    int subrun = sr.subRun();
    
    art::ServiceHandle<geo::Geometry> geom;
    unsigned int i, j;
    for (i=0; i<geom->NPlanes(); ++i) {
      int view = geom->Plane(i)->View();
      for (j=0; j<geom->Plane(i)->Ncells(); ++j) {
        fNt->Fill(run, subrun, view, i, j, fNhit[i][j]/fNevt);
      }
    }
    memset(fNhit, 0, kPLANE_MAX*kCELL_MAX*sizeof(float));
    fNevt = 0;
  }
  
  //......................................................................
  
  void ChanOcc::reconfigure(const fhicl::ParameterSet& p) 
  {
    fADCThresh = p.get<float>("ADCThresh");
  }

} // end namespace comi  
////////////////////////////////////////////////////////////////////////
DEFINE_ART_MODULE(comi::ChanOcc)
////////////////////////////////////////////////////////////////////////