T0Counter_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       T0Counter
// Module Type: producer
// File:        T0Counter_module.cc
//
// Generated at Tue Mar 24 08:49:12 2015 by Matthew Thiesse using artmod
// from cetpkgsupport v1_08_04.
////////////////////////////////////////////////////////////////////////

#include "art/Framework/Core/EDProducer.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_root_io/TFileService.h"
#include "canvas/Utilities/InputTag.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"

#include "larcore/CoreUtils/ServiceUtil.h"
#include "lardataobj/RawData/ExternalTrigger.h"
#include "lardataobj/AnalysisBase/T0.h"
#include "lardata/DetectorInfoServices/DetectorClocksService.h"
#include "lardata/Utilities/AssociationUtil.h"

#include <memory>
#include <map>
#include <limits>

#include "TTree.h"

namespace dune {
  class T0Counter;
}

class dune::T0Counter : public art::EDProducer {
public:
  explicit T0Counter(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.
  T0Counter(T0Counter const &) = delete;
  T0Counter(T0Counter &&) = delete;
  T0Counter & operator = (T0Counter const &) = delete;
  T0Counter & operator = (T0Counter &&) = delete;

  // Required functions.
  void produce(art::Event & e) override;

  // Selected optional functions.
  void beginJob() override;

private:

  // data structure to hold t0 information and ExternalTrigger data products
  struct t0 {
 
    std::vector< std::pair< unsigned int, double> > idtime;
    std::vector< art::Ptr< raw::ExternalTrigger> > trigs;

    // returns true if the trigger is within tolerance of at least one other trigger in the t0
    bool testCoincidence(double tim, double tol) const {// tol is in units of ns
      for(auto const &i : idtime) {
        if (tim >= i.second-tol && tim <= i.second+tol) return true;
      }
      return false;
    }

    // calculate average time of hit
    double avgTime() const {
      double t = 0;
      int n = 0;
      for(auto const &i : idtime) t += i.second, ++n;
      if (n != 0) return t/n;
      return std::numeric_limits<double>::min();
    }

    // insert another trigger into the t0 object
    void insert(double tim, unsigned int id, art::Ptr< raw::ExternalTrigger> trig) {
      idtime.push_back(std::pair<unsigned int,double>(id,tim));
      trigs.push_back(trig);
      return;
    }

    // constructor
    t0(double tim, unsigned int id, art::Ptr< raw::ExternalTrigger> trig) {
      insert(tim,id,trig);
    }
  };

  double tick2Time(unsigned int t);

  // fhicl parameters
  std::string fTriggerModuleLabel;
  double fClockSpeedCounter;
  double fCombinedTimeDelay;
  int fCoincidenceTolerance;
  bool fVerbose;
  bool fMakeTree;

  // calculated parameters
  double fSampleTimeCounter;
  double fTriggerOffsetTPC;

  // anab::T0 data members to fill
  double t0time;
  unsigned int trigtype;
  int trigbits;
  //int t0id; // unused

  // output TTree
  TTree *fTree;
  int run;
  int subrun;
  int event;
  int ntrigs;
  unsigned int auxdetid;
  unsigned int tick;
  double time;


};

/////////////////////////////////////////////////////////////////////////////////////

dune::T0Counter::T0Counter(fhicl::ParameterSet const & p)
  : EDProducer{p}, fTriggerModuleLabel(p.get<std::string>("TriggerModuleLabel")),
    fClockSpeedCounter(p.get<double>("ClockSpeedCounter")), // MHz
    fCombinedTimeDelay(p.get<double>("CombinedTimeDelay")), // ns
    fCoincidenceTolerance(p.get<int>("CoincidenceTolerance")), // ticks
    fVerbose(p.get<bool>("Verbose",true)),
    fMakeTree(p.get<bool>("MakeTree",false))
{
  fSampleTimeCounter = 1.e3/fClockSpeedCounter;//ns

  auto const clockData = art::ServiceHandle<detinfo::DetectorClocksService const>()->DataForJob();

  fTriggerOffsetTPC = clockData.TriggerOffsetTPC()*1.e3; // ns

  produces< std::vector< anab::T0 > >();
  produces< art::Assns< anab::T0, raw::ExternalTrigger> >();
}

/////////////////////////////////////////////////////////////////////////////////////

void dune::T0Counter::produce(art::Event & e)
{
  // make collection of T0 objects, and associations
  std::unique_ptr< std::vector< anab::T0> > t0col(new std::vector<anab::T0>);
  std::unique_ptr< art::Assns< anab::T0, raw::ExternalTrigger> > assn( new art::Assns< anab::T0, raw::ExternalTrigger>);

  // get basic info
  run = e.run();
  subrun = e.subRun();
  event = e.id().event();
  
  // get raw::ExternalTriggers
  auto externalTriggerListHandle = e.getHandle< std::vector< raw::ExternalTrigger> >(fTriggerModuleLabel);
  if (externalTriggerListHandle)
    {

      std::vector< art::Ptr< raw::ExternalTrigger> > trigs;
      art::fill_ptr_vector(trigs,externalTriggerListHandle);

      ntrigs = externalTriggerListHandle->size();

      // this vector will contain t0 objects in a nice, convenient, organised package
      std::vector<t0> t0vect;

      // fill t0vect with triggers, group if times are within fCoincidenceTolerance
      for (auto const& trig : trigs) {
        auxdetid = trig->GetTrigID();
        tick = trig->GetTrigTime();
        time = tick2Time(tick);
        if (fMakeTree) fTree->Fill();
    
        // convert # PENN ticks tolerance to nsec
        double tol = fCoincidenceTolerance*fSampleTimeCounter;

        // loop over t0vect
        std::vector<t0>::iterator t0it;
        for (t0it = t0vect.begin(); t0it != t0vect.end(); ++t0it) {      
          // get t0 object
          t0* t = &(*t0it);
          // test for coincidence with all triggers associated with t0 object
          if (t->testCoincidence(time,tol)) {
            // insert trigger into t0
            t->insert(time,auxdetid,trig);
            // stop loop after one hit to prevent multiple insertions
            break;
          }      
        }
        // if the above loop didn't find a coincidence, create a new t0 object
        if (t0it == t0vect.end()) t0vect.push_back(t0(time,auxdetid,trig));
      }
  
      // print information, if interested
      if (fVerbose) {
        int num = 0;
        for (auto const &i : t0vect) {
          std::cout << "Trigger " << num << " has " << i.idtime.size() << " triggers associated." << std::endl;
          for (auto const &tit : i.idtime) {
            std::cout << "       time=" << std::fixed << std::setprecision(1) << tit.second << "  ID=" << tit.first << std::endl;
          }
          ++num;
        }
      }


      // make anab::T0 and Assns<anab::T0, raw::ExternalTrigger>
      for (auto const &i : t0vect) {
        t0time = i.avgTime(); // calculate average time of triggers, use as T0 time
        trigtype = 1; // use 1 for counters
        trigbits = i.idtime.size(); // number of associated triggers
        t0col->push_back(anab::T0(t0time,
                                  trigtype,
                                  trigbits,
                                  (*t0col).size()
                                  ));
        util::CreateAssn(*this, e, *t0col, i.trigs, *assn);
      }                       

    }
  e.put(std::move(t0col));
  e.put(std::move(assn));

}

/////////////////////////////////////////////////////////////////////////////////////

void dune::T0Counter::beginJob()
{
  art::ServiceHandle<art::TFileService> tfs;

  if (fMakeTree) {
    fTree = tfs->make<TTree>("T0Counter","T0Counter");
    fTree->Branch("run",&run,"run/I");
    fTree->Branch("subrun",&subrun,"subrun/I");
    fTree->Branch("event",&event,"event/I");
    fTree->Branch("ntrigs",&ntrigs,"ntrigs/I");
    fTree->Branch("auxdetid",&auxdetid,"auxdetid/I");
    fTree->Branch("tick",&tick,"tick/I");
    fTree->Branch("time",&time,"time/D");
  }
}

/////////////////////////////////////////////////////////////////////////////////////

double dune::T0Counter::tick2Time(unsigned int t)
{
  double calctime = (double)t+0.5;//get central value of tick bin
  calctime *= fSampleTimeCounter;//convert tick to time in ns
  calctime += fTriggerOffsetTPC;//adjust for possible TPC window offset
  calctime -= fCombinedTimeDelay;//correct for delay in signal propagation in cables and processing
  return calctime;
}

DEFINE_ART_MODULE(dune::T0Counter)