RawDigitOverlayDUNE35tAna_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       RawDigitOverlayDUNE35tAna
// Module Type: analyzer
// File:        RawDigitOverlayDUNE35tAna_module.cc
//
// Generated at Nov 20 13:06:09 2015 by Wesley Ketchum using artmod
// from cetpkgsupport v1_08_07.
////////////////////////////////////////////////////////////////////////

#include <map>

#include "art/Framework/Core/EDAnalyzer.h"
#include "art/Framework/Core/ModuleMacros.h"
#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Principal/Run.h"
#include "art/Framework/Principal/SubRun.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "art_root_io/TFileService.h"

#include <sstream>

#include "TTree.h"

#include "larcore/Geometry/Geometry.h"

#include "DataOverlay/RawDigitAdderAna.hh"

namespace mix {
  const int kNMaxADCs = 100000000;
  class RawDigitOverlayDUNE35tAna;
  enum RawDigitTypes {
    kRawData = 0,
    kRawMC,
    kRawMixed,
    kNoType = -9999
  };
}


class mix::RawDigitOverlayDUNE35tAna : public art::EDAnalyzer {
public:
  explicit RawDigitOverlayDUNE35tAna(fhicl::ParameterSet const & p);
  // The destructor generated by the compiler is fine for classes
  // without bare pointers or other resource use.

  // Plugins should not be copied or assigned.
  RawDigitOverlayDUNE35tAna(RawDigitOverlayDUNE35tAna const &) = delete;
  RawDigitOverlayDUNE35tAna(RawDigitOverlayDUNE35tAna &&) = delete;
  RawDigitOverlayDUNE35tAna & operator = (RawDigitOverlayDUNE35tAna const &) = delete;
  RawDigitOverlayDUNE35tAna & operator = (RawDigitOverlayDUNE35tAna &&) = delete;

  // Required functions.
  void analyze(art::Event const & e) override;
  void beginJob() override;



private:

  void FillADCVariables(art::Handle<std::vector<raw::RawDigit> > digit_handle, mix::RawDigitTypes type);
  void Reset();

  std::map<mix::RawDigitTypes, std::string> fRawDigitModuleLabels;
  std::map<mix::RawDigitTypes, art::Handle<std::vector<raw::RawDigit> > > fRawDigitProducts;

  TTree *fTree;
  size_t fNADCs;
  short fADCs[mix::kNMaxADCs];
  raw::ChannelID_t fChannelIDs[mix::kNMaxADCs];
  int fTypes[mix::kNMaxADCs];
};


mix::RawDigitOverlayDUNE35tAna::RawDigitOverlayDUNE35tAna(fhicl::ParameterSet const & p)
  :
  EDAnalyzer(p)
{
  fRawDigitModuleLabels[mix::kRawData] = p.get<std::string>("DataRawDigitModuleLabel");
  fRawDigitModuleLabels[mix::kRawMC] = p.get<std::string>("MCRawDigitModuleLabel");
  fRawDigitModuleLabels[mix::kRawMixed] = p.get<std::string>("MixedRawDigitModuleLabel");
}

void mix::RawDigitOverlayDUNE35tAna::beginJob(){
  art::ServiceHandle<art::TFileService> tfs;
  fTree = tfs->make<TTree>("MRDT","MixerRawDigitTree");

  fTree->Branch("NADC",&fNADCs,"NADC/i");
  fTree->Branch("ADC",fADCs,"ADC[NADC]/S");
  fTree->Branch("ChannelID",fChannelIDs,"ChannelID[NADC]/I");
  fTree->Branch("Type",fTypes,"Type[NADC]/I");
}

void mix::RawDigitOverlayDUNE35tAna::analyze(art::Event const & e)
{

  Reset();
  art::ServiceHandle<geo::Geometry> geometry;
  std::cout<<"Num TPC channels is: " << geometry->Nchannels() << std::endl;

  for (int type_num = mix::kRawData; type_num <= mix::kRawMixed; type_num++){
    mix::RawDigitTypes type = static_cast<mix::RawDigitTypes>(type_num);
    e.getByLabel(fRawDigitModuleLabels[type],fRawDigitProducts[type]);
    if (!fRawDigitProducts[type].isValid()){
      std::cerr<<"Could not find raw digit with label: " << fRawDigitModuleLabels[type] << std::endl;
    }
    else{
      FillADCVariables(fRawDigitProducts[type],type);
    }
  }

  /*
  art::ServiceHandle<art::TFileService> tfs;

  art::Handle<std::vector<raw::RawDigit> > waveform1Handle,waveform2Handle,waveformSumHandle;
  e.getByLabel(fRawDigitModule1,waveform1Handle);
  e.getByLabel(fRawDigitModule2,waveform2Handle);
  e.getByLabel(fRawDigitModuleSum,waveformSumHandle);
  if(!waveform1Handle.isValid() || !waveform2Handle.isValid() || !waveformSumHandle.isValid()) throw std::runtime_error("some handle is not valid");
  std::vector<raw::RawDigit> const& waveform1Vector(*waveform1Handle); 
  std::vector<raw::RawDigit> const& waveform2Vector(*waveform2Handle); 
  std::vector<raw::RawDigit> const& waveformSumVector(*waveformSumHandle);

  size_t histos_to_make = fAnaAlg.CheckOverlay(waveform1Vector,waveform2Vector,waveformSumVector);

  std::vector<TH1S*> histograms(histos_to_make);
  for(size_t i_h=0; i_h<histograms.size(); i_h++){
    std::stringstream hname; hname << "h" << i_h;
    tfs->make<TH1S>(hname.str().c_str(),"GenericTitle",1,0,1);
  }
  
  fAnaAlg.CreateOutputHistograms(histograms,
                                 waveform1Vector,waveform2Vector,waveformSumVector,
                                 e.run(),e.event());
     */

  fTree->Fill();
}

void mix::RawDigitOverlayDUNE35tAna::FillADCVariables(art::Handle<std::vector<raw::RawDigit> > digit_handle, mix::RawDigitTypes type){
  //Check if the handle is valid.  If not, don't do anything :(
  if (!digit_handle.isValid()){
    std::cerr<<"Digit handle for type " << type << " not valid!  Don't attempt to fill the tree for this product"<<std::endl;
    return;
  }
  std::vector<raw::RawDigit> const& digitVector(*digit_handle);

  for (unsigned int i_raw = 0; i_raw < digitVector.size(); i_raw++){
    raw::RawDigit raw_digit = digitVector[i_raw];

    raw::ChannelID_t channel = raw_digit.Channel();

    for (unsigned int i_adc = 0; i_adc < raw_digit.ADCs().size(); i_adc++){
      fADCs[fNADCs] = raw_digit.ADCs()[i_adc];
      fChannelIDs[fNADCs] = channel;
      fTypes[fNADCs] = static_cast<int>(type);
      fNADCs++;
    }
  }

  return;
}

void mix::RawDigitOverlayDUNE35tAna::Reset(){
  fNADCs = 0;
  for (int i = 0; i < mix::kNMaxADCs; i++){
    fADCs[i] = -9999;
    fChannelIDs[i] = -9999;
    fTypes[i] = kNoType;
  }
}

DEFINE_ART_MODULE(mix::RawDigitOverlayDUNE35tAna)