RecordNoiseSpectrumFile_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
/// \brief   TODO
/// \author Christopher Backhouse - bckhouse@caltech.edu
////////////////////////////////////////////////////////////////////////

// Framework includes
#include "art/Framework/Core/EDAnalyzer.h"
#include "canvas/Persistency/Common/Ptr.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/SubRun.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "fhiclcpp/ParameterSet.h"

// NOvA includes
#include "ChannelInfo/BadChanList.h"
#include "Geometry/Geometry.h"
#include "RecoBase/CellHit.h"
#include "RecoBase/Cluster.h"
#include "SummaryData/CosmicExposure.h"

#include "TH3.h"

namespace rsim
{
  class RecordNoiseSpectrumFile: public art::EDAnalyzer
  {
  public:
    explicit RecordNoiseSpectrumFile(const fhicl::ParameterSet& pset);
    virtual ~RecordNoiseSpectrumFile();

    virtual void beginJob();
    virtual void analyze(const art::Event& evt);
    virtual void endSubRun(const art::SubRun& sr);
    virtual void endJob();

  protected:
    int fSubrunCount; ///< We can only handle input of precisely one subrun
    TH3* fNoiseMap;
  };

  //----------------------------------------------------------------------
  RecordNoiseSpectrumFile::RecordNoiseSpectrumFile(const fhicl::ParameterSet& pset) :
    EDAnalyzer(pset),
    fSubrunCount(0)
  {
  }

  //----------------------------------------------------------------------
  RecordNoiseSpectrumFile::~RecordNoiseSpectrumFile()
  {
  }

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

    fNoiseMap = tfs->make<TH3D>("h3_adc0adc1adc2", ";ADC0;ADC1;ADC2",
                                150, 0, 150,
                                150, 0, 150,
                                150, 0, 150);
  }

  //----------------------------------------------------------------------
  void RecordNoiseSpectrumFile::analyze(const art::Event& evt)
  {
    art::Handle<std::vector<rb::Cluster> > slices;
    evt.getByLabel("slicer", slices);

    const unsigned int slcMax = slices->size();
    for(unsigned int slcIdx = 0; slcIdx < slcMax; ++slcIdx){
      art::Ptr<rb::Cluster> slice(slices, slcIdx);

      if(!slice->IsNoise()) continue;

      const unsigned int cellMax = slice->NCell();
      for(unsigned int cellIdx = 0; cellIdx < cellMax; ++cellIdx){
        art::Ptr<rb::CellHit> chit = slice->Cell(cellIdx);

        fNoiseMap->Fill(chit->ADC(0)+0.5, 0.0, 0.0);
      } // end for cellIdx
    } // end for slcIdx
  }

  //----------------------------------------------------------------------
  void RecordNoiseSpectrumFile::endSubRun(const art::SubRun& sr)
  {
    ++fSubrunCount;

    art::Handle<sumdata::CosmicExposure> ce;
    sr.getByLabel("expo", ce);

    mf::LogInfo("RecordNoiseSpectrumFile") << "Total live time is " << 1e9*ce->totlivetime << " ns";

    art::ServiceHandle<geo::Geometry> geom;
    art::ServiceHandle<chaninfo::BadChanList> badChans;

    unsigned int nLiveChans = 0;

    const std::set<unsigned int>& planes = geom->GetPlanesByView(geo::kXorY);
    for(std::set<unsigned int>::const_iterator it = planes.begin();
        it != planes.end(); ++it){
      const unsigned int plane = *it;
      const unsigned int maxCell = geom->Plane(plane)->Ncells();
      for(unsigned int cell = 0; cell < maxCell; ++cell){
        if(!badChans->IsBad(plane, cell)) ++nLiveChans;
      } // end for cell
    } // end for it
    mf::LogInfo("RecordNoiseSpectrumFile") << "There are " << nLiveChans << " live channels";

    // Normalization of noise map is probability /cell /ns
    fNoiseMap->Scale(1.0/(nLiveChans*1e9*ce->totlivetime));
  }

  //----------------------------------------------------------------------
  void RecordNoiseSpectrumFile::endJob()
  {
    assert(fSubrunCount == 1 && "Must run over precisely one subrun. "
           "We're not smart enough to do the normalization in other cases");
  }

  DEFINE_ART_MODULE(RecordNoiseSpectrumFile)

}//namespace