MergeDaqCollections_module.cc

Go to the documentation of this file.

///////////////////////////////////////////////////////////////////////////
// \file    MergeDaqCollections_module.cxx
// \brief   Merge collections.
// At this step collection from primary and secondary files are merged under the same label. You can set RawDigit::SetMC(true or false). Which will let the calibrator use mc or real calibration constants when you analyze the merged file. 
// I learned from MergeCollections_module.cc
// \author Enhao Song es4re@virginia.edu
/////////////////////////////////////////////////////////////////////////// 

#include <math.h>
#include <string>
#include <algorithm>
#include <limits>

#include "art/Framework/Principal/Event.h"
#include "art/Framework/Principal/Handle.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"
#include "fhiclcpp/ParameterSet.h"
#include "messagefacility/MessageLogger/MessageLogger.h"
#include "art/Framework/Core/EDProducer.h"
#include "art/Framework/Core/ModuleMacros.h"

#include "DAQChannelMap/DAQChannelMap.h"
#include "RawData/RawDigit.h"
#include "RawData/RawTrigger.h"
#include "RawData/DAQHeader.h"
#include "RecoBase/CellHit.h"
#include "Calibrator/Calibrator.h"
#include "MCCheater/BackTracker.h"

#include "Utilities/AssociationUtil.h"

///group hits according to time and space
namespace MergeDaqCollections {

  class MergeDaqCollections : public art::EDProducer {
    public:
      explicit MergeDaqCollections(fhicl::ParameterSet const &pset);
      virtual ~MergeDaqCollections();

      void produce(art::Event& evt);

    protected:
      // internal methods

      // configuration settings
      std::string  fFirstDaqCollection;     /// Label on the main collection
      std::string  fSecondCollection;       /// Label on the main collection

      bool fFirstCollectionIsMC;            /// First collection is MC
      bool fSecondCollectionIsMC;           /// Second collection is MC
      bool fUseSecondaryTrigger;            /// Want trigger from 2nd collection rather than 1st

      bool fRemoveFirstCollectionNoise;     /// Don't mix in noise hits from first collection, if it is MC.
      bool fRemoveSecondCollectionNoise;    /// Don't mix in noise hits from second collection, if it is MC.
      bool fOverlayWithTDCCut;              /// Only overlay secondary hits within a TDC range of the main event hits. Since we are using the secondary hits as background noise hits. We only need them in the time window close to the signal hits. We can still use the cosmic files before overlay for as a larger data set for background study. This way we can save a lot of disk space for the overlaid files and computing time in the following analysis chain.
      int fTDCExtendValue;

      std::string fInstanceLabel;           /// Instance label
  };


  //----------------------------------------------------------------------
  MergeDaqCollections::MergeDaqCollections(fhicl::ParameterSet const &pset) :
    fFirstDaqCollection         (pset.get< std::string >("FirstDaqCollection")),
    fSecondCollection           (pset.get< std::string >("SecondCollection")),
    fFirstCollectionIsMC        (pset.get< bool >("FirstCollectionIsMC")),
    fSecondCollectionIsMC       (pset.get< bool >("SecondCollectionIsMC")),
    fUseSecondaryTrigger        (pset.get< bool >("UseSecondaryTrigger")),
    fRemoveFirstCollectionNoise (pset.get< bool >("RemoveFirstCollectionNoise")),
    fRemoveSecondCollectionNoise(pset.get< bool >("RemoveSecondCollectionNoise")),
    fOverlayWithTDCCut(pset.get< bool >("OverlayWithTDCCut")),
    fTDCExtendValue(pset.get< int >("TDCExtendValue")),
    fInstanceLabel              (pset.get< std::string >("InstanceLabel"))
  {
    produces< std::vector<rawdata::RawDigit>                   >(fInstanceLabel);
    produces< std::vector<rawdata::RawTrigger>                 >(fInstanceLabel);
  }

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

  //----------------------------------------------------------------------
  void MergeDaqCollections::produce(art::Event& evt)
  {
    art::ServiceHandle<calib::Calibrator> cal;
    art::ServiceHandle<cheat::BackTracker> bt;

    // BackTracker has not yet seen MC mixed from secondary file
    if(fSecondCollectionIsMC) bt->Rebuild(evt);

    std::unique_ptr< std::vector<rawdata::RawDigit> > digitcol (new std::vector<rawdata::RawDigit>);
    std::unique_ptr< std::vector<rawdata::RawTrigger> > triggercol (new std::vector<rawdata::RawTrigger>);

    art::Handle< std::vector<rawdata::RawDigit> > digit1list;
    evt.getByLabel(fFirstDaqCollection, digit1list);

    art::Handle< std::vector<rawdata::RawDigit> > digit2list;
    evt.getByLabel(fSecondCollection, digit2list);   


    // First Collection Digits
    auto min_TDC = std::numeric_limits<int32_t>::max();
    auto max_TDC = std::numeric_limits<int32_t>::min();
    for(unsigned int i = 0; i < digit1list->size(); ++i){     
      rawdata::RawDigit digit ((*digit1list)[i]);
      digit.SetMC(fFirstCollectionIsMC);
      if(fFirstCollectionIsMC && fRemoveFirstCollectionNoise){
        // Make cell hit so backtracker can 
        // attempt to match to particle
        if(!bt->IsNoise(cal->MakeCellHit(&digit)))
          digitcol->push_back(digit);
      }
      else
        digitcol->push_back(digit);
    }
    //find max and min TDC of the digitcol
    if (fOverlayWithTDCCut){
      for(auto digit: *digitcol){
        min_TDC = std::min(min_TDC, digit.TDC());
        max_TDC = std::max(max_TDC, digit.TDC());
      }
    }

    //extend time window on each end with fTDCExtendValue
    min_TDC = std::max(min_TDC-fTDCExtendValue, -3200); // data always has TDC > -3200 TDC which is 50 us, so we only overlay data beyond this time point.
    max_TDC = max_TDC+fTDCExtendValue;

    // Second Collection Digits   
    for(unsigned int i = 0; i < digit2list->size(); ++i){
      rawdata::RawDigit digit ((*digit2list)[i]);
      if (fOverlayWithTDCCut){
        if (digit.TDC()>max_TDC || digit.TDC()<min_TDC)
          continue;
      }
      digit.SetMC(fSecondCollectionIsMC);     
      if(fSecondCollectionIsMC && fRemoveSecondCollectionNoise){
        // Make cell hit so backtracker can 
        // attempt to match to particle
        if(!bt->IsNoise(cal->MakeCellHit(&digit)))
          digitcol->push_back(digit);
      }
      else
        digitcol->push_back(digit);
    }


    // It can't be assumed which trigger is appropriate
    art::Handle< std::vector<rawdata::RawTrigger> > triggerlist;

    if(fUseSecondaryTrigger){
      evt.getByLabel(fSecondCollection, triggerlist);

      for (unsigned int i = 0; i < triggerlist->size(); ++i) {
        rawdata::RawTrigger trigger ((*triggerlist)[i]);
        triggercol->push_back(trigger);
      }
    } else {
      evt.getByLabel(fFirstDaqCollection, triggerlist);

      for (unsigned int i = 0; i < triggerlist->size(); ++i) {
        rawdata::RawTrigger trigger ((*triggerlist)[i]);
        triggercol->push_back(trigger);
      }
    }

    // put the collections in the event
    evt.put(std::move(digitcol),   fInstanceLabel);
    evt.put(std::move(triggercol), fInstanceLabel);

  }

  DEFINE_ART_MODULE(MergeDaqCollections)
} // end namespace MergeDaqCollections
/////////////////////////////////////////////////////////////////////////