MichelETrigger_module.cc

Go to the documentation of this file.

////////////////////////////////////////////////////////////////////////
// Class:       MichelETrigger
// Module Type: filter
// File:        MichelETrigger_module.cc
//
// Generated at Wed Jul 19 10:03:20 2017 by Matthew Judah using artmod
// from cetpkgsupport v1_10_02.
////////////////////////////////////////////////////////////////////////


#include "art/Framework/Core/EDFilter.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 "DDTUtilities/AssociationUtil.h"
#include "canvas/Persistency/Common/FindManyP.h"
#include "canvas/Persistency/Common/FindOneP.h"
#include "art/Framework/Services/Optional/TFileService.h"
#include "art/Framework/Services/Optional/TFileDirectory.h"
#include "art/Framework/Services/Registry/ServiceHandle.h"

#include "MichelTriggerPrescaleTable.h"

//DDT Includes
#include "DDTBaseDataProducts/Track3D.h"
#include "DDTBaseDataProducts/BaseProducts.h"
#include "DDTBaseDataProducts/HitList.h"
#include "DDTBaseDataProducts/DAQHit.h"
#include "DDTBaseDataProducts/HitSet.h"
#include "DDTBaseDataProducts/CompareDAQHit.h"
#include "DDTBaseDataProducts/TriggerDecision.h"

//C++ Includes
#include <memory>
#include <sstream>
#include <string>
#include <TMath.h>
#include <algorithm>
#include <numeric>
#include <iomanip>
#include <iostream>
#include <list>
#include <vector>       
#include <cmath>

//ROOT Includes
#include "TH1F.h"
#include "TH2F.h"
#include "TVector3.h"


#include <memory>

namespace novaddt {
  class MichelETrigger;
}

class novaddt::MichelETrigger : public art::EDFilter {
public:

  int nEvents        = 0;
  int nTracks        = 0;
  int nStoppedTracks = 0;
  int nClusters1     = 0;
  int nMichels       = 0;
  int nTrig          = 0;

  explicit MichelETrigger(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.
  MichelETrigger(MichelETrigger const &) = delete;
  MichelETrigger(MichelETrigger &&) = delete;
  MichelETrigger & operator = (MichelETrigger const &) = delete;
  MichelETrigger & operator = (MichelETrigger &&) = delete;

  // Required functions.
  bool filter(art::Event & e) override;

  float GetSmartPrescale(int xcell, int ycell);
  bool GetSmartPrescaleDecision(double prescale);
  // Selected optional functions.

  void endJob() override;
  void reconfigure(fhicl::ParameterSet const & p);

private:
  const int kTDC_PER_USEC = 64;
  const double cw = 3.97;
  const double pw = 6.654;

  // Declare member data here.
  unsigned _prescale;
  std::string _tracklabel;
  std::string _hitslabel;
  double fTimeUpperTrack;
  double fTimeLowerTrack;
  double fTimeUpperCluster;
  double fTimeLowerCluster;
  double fDistMuon;
  double fDistMichel;
  double fMaxE;
  int fFiducial_Z_min;
  int fFiducial_Z_max;
  int fFiducial_X_min;
  int fFiducial_X_max;
  int fFiducial_Y_min;
  int fFiducial_Y_max;
  bool fapplySmartPrescale;
  double fadditionalPrescale;
  int _trigger_counts = 0;
};


novaddt::MichelETrigger::MichelETrigger(fhicl::ParameterSet const & p)
:
  _prescale         (p.get<unsigned>("prescale")           ),
  _tracklabel       (p.get<std::string>("TrackLabel")      ),
  _hitslabel        (p.get<std::string>("HitsTag")         ),
  fTimeUpperTrack   (p.get<double>("time1_u")*kTDC_PER_USEC),
  fTimeLowerTrack   (p.get<double>("time1_l")*kTDC_PER_USEC),
  fTimeUpperCluster (p.get<double>("time2_u")*kTDC_PER_USEC),
  fTimeLowerCluster (p.get<double>("time2_l")*kTDC_PER_USEC),
  fDistMuon         (p.get<double>("MaxDistance")          ),
  fDistMichel       (p.get<double>("DistanceMichel")       ),
  fMaxE             (p.get<double>("EnergyCutOff")         ),
  fFiducial_Z_min   (p.get<int>("fiducial.ZetRange.min")   ),
  fFiducial_Z_max   (p.get<int>("fiducial.ZetRange.max")   ),
  fFiducial_X_min   (p.get<int>("fiducial.XRange.min")     ),
  fFiducial_X_max   (p.get<int>("fiducial.XRange.max")     ), 
  fFiducial_Y_min   (p.get<int>("fiducial.YRange.min")     ),            
  fFiducial_Y_max   (p.get<int>("fiducial.YRange.max")     ),
  fapplySmartPrescale(p.get<bool>("applySmartPrescale")     ),
  fadditionalPrescale(p.get<double>("additionalPrescale")   )
{
  produces<std::vector<TriggerDecision>>();
}

bool cmp_tdc(const novaddt::DAQHit& h, const novaddt::TDC& tdc){
  return h.TDC() < tdc;
}

float novaddt::MichelETrigger::GetSmartPrescale(int xcell, int ycell)
{
  //m and b come from the conversion of 383 cells (X/Y dimension of FD)
  //to 95 cells
  float m = 62./246.;
  float b = 14 - 54*m;
  float xcell_convert = m*xcell + b;
  float ycell_convert = m*ycell + b;
  if(xcell_convert > 94) xcell_convert = 94;
  if(ycell_convert > 94) ycell_convert = 94;
  return prescale_table[floor(xcell_convert)][floor(ycell_convert)]
    * fadditionalPrescale; 
}


bool novaddt::MichelETrigger::GetSmartPrescaleDecision(double prescale) 
  {
    // Is it good enough for small prescale factors ~1/100?     
    double r = (double)rand() / (double)RAND_MAX;
    return r < prescale;
  }


bool novaddt::MichelETrigger::filter(art::Event & e)
{
  nEvents++;

  std::unique_ptr<std::vector<TriggerDecision>>
    trigger_decisions(new std::vector<TriggerDecision>);
  bool isTriggering = false;

  art::Handle<std::vector<Track3D>> Tracks;
  e.getByLabel(_tracklabel, Tracks);

  //get the hit list for each track, use this to calculate mean track time
  art::FindOneP<novaddt::HitList> fohl(Tracks,e,_tracklabel);
  assert(fohl.isValid());
  assert(fohl.size() == Tracks->size());

  //Grab all time sorted hits
  auto cellHits = e.getValidHandle<novaddt::HitList>(_hitslabel);

  //This will store the Michel Clusters id'd at the end of a track
  std::vector<std::vector<novaddt::DAQHit>> result;

  //Loop over all the tracks in this event
  for(uint iii = 0; iii < Tracks->size(); ++iii){
    nTracks++;
    //Grab track
    novaddt::Track3D track = Tracks->at(iii);

    //Make Sure track is 3D
    if(!track.Is3D()) continue;

    TVector3 trk_end[2];
    trk_end[0] = track.Start();
    trk_end[1] = track.End();

    for(int which_end = 0; which_end < 2; which_end++){
      if( (trk_end[which_end].X() < fFiducial_X_min) ||
          (trk_end[which_end].X() > fFiducial_X_max)) continue;
      if( (trk_end[which_end].Y() < fFiducial_Y_min) ||
          (trk_end[which_end].Y() > fFiducial_Y_max)) continue; 
      if( (trk_end[which_end].Z() < fFiducial_Z_min) ||
          (trk_end[which_end].Z() > fFiducial_Z_max)) continue;

      nStoppedTracks++;

      //Track stops in the detector. Start looking for the Michel Cluster
      std::vector<novaddt::DAQHit> meCluster;

      //Get Average time of the stopped track
      //get the hit list for this track
      art::Ptr<novaddt::HitList> this_hit_list = fohl.at(iii);

      double avg_t = 0;
      double totW  = 0;
      double min_TDC = 0;
      double max_TDC = 0;

      for(uint trackHits = 0; trackHits < this_hit_list->size(); trackHits++){
        if(trackHits == 0)
          min_TDC = this_hit_list->at(trackHits).TDC().val;
        if( this_hit_list->at(trackHits).ADC().val < 50) continue;
        if( this_hit_list->at(trackHits).TDC().val < min_TDC)
          min_TDC = this_hit_list->at(trackHits).TDC().val;
        double w = this_hit_list->at(trackHits).ADC().val;
        avg_t += this_hit_list->at(trackHits).TDC().val * w;
        totW += w;
      }// loop over track hits

      avg_t = avg_t/totW;


      //Get start and end cell hits that fall in the time window
      const double firstT = (avg_t+fTimeLowerTrack);
      const double lastT  = (avg_t+fTimeUpperTrack);

      std::vector<novaddt::DAQHit> daqHits = *(cellHits);


      const auto cellHitsBegin = std::lower_bound(daqHits.begin(),
                                                   daqHits.end(), firstT,
                                                   cmp_tdc);
      const auto cellHitsEnd   = std::lower_bound(cellHitsBegin, daqHits.end(),
                                                  lastT, cmp_tdc);

      //std::cout << avg_t << "\t" << firstT << "\t" << lastT << std::endl;

      //Hold onto potential Michel Cluster hits
      std::vector<novaddt::DAQHit> meList;
      double michel_pos[3] = {0,0,0};
      int seed_plane_pos = -5;
      double seed_time   = 0;
      double maxE = 0;
      double distance_min = 1000;
      for(auto cellHitsIter = cellHitsBegin; cellHitsIter != cellHitsEnd; ++cellHitsIter)
        {
          novaddt::DAQHit c = *cellHitsIter;
          if(c.ADC() < 50) continue;

          //      std::cout << c.TDC().val << std::endl;

          double cellpos[3] = {0,0,0};
          cellpos[0] = c.Cell().val;
          cellpos[1] = c.Cell().val;
          cellpos[2] = c.Plane().val;

          if(c.View().val == daqchannelmap::X_VIEW)
            cellpos[1] = trk_end[which_end].Y();
          else
            cellpos[0] = trk_end[which_end].X();

          TVector3 distance;
          distance.SetXYZ((cellpos[0] - trk_end[which_end].X())*cw,
                          (cellpos[1] - trk_end[which_end].Y())*cw,
                          (cellpos[2] - trk_end[which_end].Z())*pw);

          double dist = distance.Mag();

          if(dist < fDistMuon && c.ADC().val < fMaxE){
            novaddt::DAQHit hit = c;
            meList.push_back(hit);
            //std::cout << (hit.TDC().val - avg_t)/kTDC_PER_USEC << "\t" <<
            //  dist << "\t" << hit.ADC().val << std::endl;
            if(hit.ADC() > maxE && dist < distance_min){
              seed_time = hit.TDC().val;
              seed_plane_pos = cellpos[2];
              maxE = hit.ADC().val;
              michel_pos[0] = cellpos[0];
              michel_pos[1] = cellpos[1];
              michel_pos[2] = cellpos[2];
            }
          }
        }// First loop over cell hits

      //Didn't find any me candidates don't bother looking any further
      if(seed_plane_pos < 0) continue;
      nClusters1++;

      max_TDC = seed_time;// fTimeUpperCluster;
      double dist2 = 1000;
      auto iter = meList.begin();
      while( iter != meList.end()){
        novaddt::DAQHit c = *iter;

        double plane = c.Plane().val;

        dist2 = abs( plane - seed_plane_pos)*pw;

        //Make second timing cut
        //Select only cell hits that happen
        //within a certain time of the michel e candidate hit

        if( ((seed_time - fTimeLowerCluster) < (*iter).TDC().val ) &&
            ((seed_time + fTimeUpperCluster) > (*iter).TDC().val ) &&
            (dist2 < fDistMichel)){
          meCluster.push_back(*iter);
          if ((*iter).TDC().val > max_TDC)
            max_TDC = (*iter).TDC().val;
        }
        iter++;
      }//end second loop on michel hits

      //Number Cut: Only Keep Michel Cluster at least one cell
      //hit and has cell hits in both views
      int x_hits = 0;
      int y_hits = 0;
      double TotalADC = 0;
      if(meCluster.size() > 0){
        for(uint ihit = 0; ihit < meCluster.size(); ihit++){
          if(meCluster[ihit].View().val == daqchannelmap::X_VIEW)
            x_hits++;
          else
            y_hits++;
          TotalADC += meCluster[ihit].ADC().val;
        }
      }

      if( x_hits > 0 && y_hits > 0){
        result.push_back(meCluster);
        _trigger_counts++;
        nMichels++;
        //std::cout << michel_pos[0] << "\t" << michel_pos[1] << "\t" <<
        //michel_pos[2] << "\t" << std::endl;
      }
      if((_trigger_counts%_prescale == _prescale-1) && (x_hits > 0) && (y_hits > 0)) {
        if(fapplySmartPrescale){
          double prescaleValue = GetSmartPrescale(michel_pos[0],michel_pos[1]);
          bool shouldTrigger = GetSmartPrescaleDecision(prescaleValue);
          if(shouldTrigger){
            trigger_decisions->emplace_back(min_TDC, max_TDC - min_TDC, 
                                daqdataformats::trigID::TRIG_ID_DATA_MICHEL,
                                _prescale);
            nTrig++;
          }
        }
        else{
          trigger_decisions->emplace_back(min_TDC, max_TDC - min_TDC, daqdataformats::trigID::TRIG_ID_DATA_MICHEL, _prescale);
          nTrig++;
        }
      }
    }// loop over ends of track   
  }//end loop over tracks
  
  if(result.size() > 0 && trigger_decisions->size() > 0) isTriggering = true;
  e.put(std::move(trigger_decisions));
  return isTriggering;
} //end Filter()



void novaddt::MichelETrigger::endJob()
{
  std::cout << "=== novaddt::MichelETrigger endJob" << std::endl;
  std::cout << "\tNumber of events:               " << nEvents    << std::endl; 
  std::cout << "\tNumber of tracks:               " << nTracks    << std::endl; 
  std::cout << "\tNumber of stopped tracks:       " << nStoppedTracks 
            << std::endl; 
  std::cout << "\tNumber of cluster seeds:        " << nClusters1 << std::endl; 
  std::cout << "\tNumber of Michels tagged:       " << nMichels   << std::endl;
  std::cout << "\tNumber of Michels Triggered:    " << nTrig      << std::endl;
}



void novaddt::MichelETrigger::reconfigure(fhicl::ParameterSet const & p)
{
  // Implementation of optional member function here.
}

DEFINE_ART_MODULE(novaddt::MichelETrigger)