HighEnergyAna_module.cc

Go to the documentation of this file.

#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 "art/Framework/Services/Optional/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"

#include "DDTBaseDataProducts/CompareDAQHit.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/TriggerDecision.h"

#include <TH1.h>
#include <TH2.h>

#include <tuple>

namespace highe {
  class HighEnergyAna;
}

class highe::HighEnergyAna : public art::EDAnalyzer {
public:
  explicit HighEnergyAna(fhicl::ParameterSet const & p);
  virtual ~HighEnergyAna();

  void analyze(art::Event const & e) override;

  void beginJob() override;

private:
  typedef std::tuple<unsigned, double, double> Binning;

  void make_histo
  (std::string const name, std::string const title, Binning const& bin);

  std::string slice_label_, slice_instance_;
  std::map<std::string, TH1*> h_;
  std::map<std::string, TH2*> h2_;
  art::ServiceHandle<art::TFileService> tfs_;
};


highe::HighEnergyAna::HighEnergyAna(fhicl::ParameterSet const & p) :
  EDAnalyzer(p),
  slice_label_(p.get<std::string>("slice_label")),
  slice_instance_(p.get<std::string>("slice_instance"))
{
}

highe::HighEnergyAna::~HighEnergyAna()
{
}

void highe::HighEnergyAna::beginJob()
{
  Binning adc_bin(4096, -0.5, 4095.5);
  Binning nbin(2001, -0.5, 2000.5);
  Binning plane_bin(1001, -0.5, 1000.5);
  Binning cell_bin(401, -0.5, 400.5);

  make_histo("counter", "Counter", Binning(10, 0.5, 10.5));

  make_histo("adc_all_hits", "ADC Distribution of All Hits", adc_bin);
  make_histo("adc_all_hits_raw", "Raw ADC Distribution of All Hits", adc_bin);
  make_histo("adc_mean", "Mean ADC per Slice", adc_bin);
  make_histo("adc_max", "Maximum ADC per Slice", adc_bin);
  make_histo("adc_min", "Minimum ADC per Slice", adc_bin);
  make_histo("adc_sum", "Summed ADC per Slice", Binning(1e4, 0, 1e8));
  make_histo("adc_sum_above_200", "Summed ADC (ADC > 200) per Slice", 
             Binning(1e4, 0, 1e8));
  
  make_histo("n_hits_per_slice", "Number of Hits per Slice", nbin);
  make_histo("n_hits_per_triggered_slice", 
             "Number of Hits per Triggered Slice", nbin);
  make_histo("n_hits_above_200", "Number of Hits above 200 ADC", nbin);
  make_histo("n_hits_above_700", "Number of Hits above 700 ADC", nbin);
  make_histo("n_hits_above_1000", "Number of Hits above 1000 ADC", nbin);
  make_histo("n_hits_above_2000", "Number of Hits above 2000 ADC", nbin);

  make_histo("plane_all_hits", "Plane of All Hits", plane_bin);
  make_histo("cell_all_hits", "Cell of All Hits", cell_bin);


  h2_["trigger_rate_700"] = tfs_->make<TH2F>
    ("trigger_rate_700", "Trigger Rate",
     5, 0, 1000,
     6, 0, 30);
  h2_["trigger_rate_1000"] = tfs_->make<TH2F>
    ("trigger_rate_1000", "Trigger Rate",
     5, 0, 1000,
     6, 0, 30);
  h2_["trigger_rate_2000"] = tfs_->make<TH2F>
    ("trigger_rate_2000", "Trigger Rate",
     5, 0, 1000,
     6, 0, 30);
}

void highe::HighEnergyAna::make_histo
(std::string const name, std::string const title, Binning const& bin)
{
  h_[name] = tfs_->make<TH1F>
    (name.c_str(), title.c_str(), 
     std::get<0>(bin), std::get<1>(bin), std::get<2>(bin));
}


void highe::HighEnergyAna::analyze(art::Event const & e)
{
  art::Handle<novaddt::HitList> sorted_hits;
  e.getByLabel("tdcsort", "TDCSortedHits", sorted_hits);

  art::Handle<std::vector<novaddt::HitList> > slices;
  e.getByLabel(slice_label_, slice_instance_, slices);

  art::Handle<std::vector<novaddt::TriggerDecision> > decisions;
  e.getByLabel("highetrigger", decisions);
  
  h_["counter"]->Fill("event", 1);

  int y;
  
  if (!decisions.failedToGet())
  {
    if (decisions->size() == 1)
      h_["counter"]->Fill("HE triggered", 1);
    
    if (decisions->size() > 1)
      h_["counter"]->Fill("multiple HE triggered", 1);
  }
    
  for (auto const& hit : *sorted_hits)
    h_["adc_all_hits_raw"]->Fill(hit.ADC().val);
  
  
  for (auto const& slice : *slices)
  {
    // h_["counter"]->Fill("slice", 1);

    double adc_sum = 0.0;
    double adc_sum_above_200 = 0.0;
    std::map<unsigned, unsigned> n_hits = 
      {{0, 0}, {200, 0}, {700, 0}, {1000, 0}, {2000, 0}};

    for (auto const& hit : slice)
    {
      // h_["counter"]->Fill("hit", 1);

      h_["plane_all_hits"]->Fill(hit.Plane().val);
      h_["cell_all_hits"]->Fill(hit.Cell().val);

      auto adc = hit.ADC().val;
      h_["adc_all_hits"]->Fill(adc);
      adc_sum += adc;

      if (adc >= 200)
        adc_sum_above_200 += adc;

      // we are looking for hits that have an energy above the thresholds 
      // in n_hits
      for (auto & n : n_hits)
        if (adc >= n.first)
          ++n.second;
    }

    double adc_mean;
    unsigned adc_min, adc_max;
    if (slice.size() > 0)
    {
      adc_mean = adc_sum / double(slice.size());
      auto range = 
        std::minmax_element(slice.begin(), slice.end(), 
                            novaddt::CompareDAQHit<novaddt::ADC>());
      adc_min = range.first->ADC().val;
      adc_max = range.second->ADC().val;
    } else {
      adc_mean = 0.0;
      adc_min = 0;
      adc_max = 0;
    }

    h_["adc_mean"]->Fill(adc_mean);
    h_["adc_min"]->Fill(adc_min);
    h_["adc_max"]->Fill(adc_max);
    h_["adc_sum"]->Fill(adc_sum);
    h_["adc_sum_above_200"]->Fill(adc_sum_above_200);

    h_["n_hits_per_slice"]->Fill(n_hits[0]);
    h_["n_hits_above_200"]->Fill(n_hits[200]);
    h_["n_hits_above_700"]->Fill(n_hits[700]);
    h_["n_hits_above_1000"]->Fill(n_hits[1000]);
    h_["n_hits_above_2000"]->Fill(n_hits[2000]);

    h2_["trigger_rate_700"]->Fill(adc_mean, n_hits[700]);
    h2_["trigger_rate_1000"]->Fill(adc_mean, n_hits[1000]);
    h2_["trigger_rate_2000"]->Fill(adc_mean, n_hits[2000]);

    // if this slice is triggered
    if (n_hits[1000] >= 15 && 
        adc_mean >= 600)
      h_["n_hits_per_triggered_slice"]->Fill(n_hits[0]);
  }
}

DEFINE_ART_MODULE(highe::HighEnergyAna)