diff --git a/TrigScint/include/TrigScint/AnalyticalRecHitProducer.h b/TrigScint/include/TrigScint/AnalyticalRecHitProducer.h
new file mode 100644
index 000000000..cdae47907
--- /dev/null
+++ b/TrigScint/include/TrigScint/AnalyticalRecHitProducer.h
@@ -0,0 +1,118 @@
+/**
+ * @file AnalyticalRecHitProducer.h
+ * @brief Class that builds recHits
+ * @author Andrew Whitbeck, TTU
+ */
+
+#ifndef TRIGSCINT_TRIGSCINTDIGIPRODUCER_H
+#define TRIGSCINT_TRIGSCINTDIGIPRODUCER_H
+
+/*~~~~~~~~~~*/
+/*   ROOT   */
+/*~~~~~~~~~~*/
+#include "TRandom3.h"
+#include "TVectorD.h"
+
+// LDMX
+#include "DetDescr/TrigScintID.h"
+#include "Recon/Event/EventConstants.h"
+#include "Tools/NoiseGenerator.h"
+#include "TrigScint/Event/TrigScintHit.h"
+#include "TrigScint/Event/TrigScintQIEDigis.h"
+
+/*~~~~~~~~~~~~~~~*/
+/*   Framework   */
+/*~~~~~~~~~~~~~~~*/
+#include "Framework/Configure/Parameters.h"
+#include "Framework/EventProcessor.h"
+
+/*~~~~~~~~~~~*/
+/* TrigScint */
+/*~~~~~~~~~~~*/
+#include "TrigScint/SimQIE.h"
+
+namespace trigscint {
+
+/**
+ * @class AnalyticalRecHitProducer
+ * @brief Organizes digis into TrigScintHits, linearizes TDC
+ * and ADC info, and converts amplitudes to PEs
+
+ */
+class AnalyticalRecHitProducer : public framework::Producer {
+ public:
+  AnalyticalRecHitProducer(const std::string& name, framework::Process& process);
+
+  ~AnalyticalRecHitProducer();
+
+  /**
+   * Callback for the processor to configure itself from the given set
+   * of parameters.
+   *
+   * @param parameters ParameterSet for configuration.
+   */
+  void configure(framework::config::Parameters& parameters) final override;
+
+  void produce(framework::Event& event);
+
+ private:
+  /**
+   * Reconstruct true charge deposited in each time sample
+   * @param adc array of adcs for give event, cell
+   * @param tdc array of tdcs for give event, cell
+   * @param sample sample of interest
+   */
+  std::vector<double>  ChargeReconstruction(std::vector<int>adc
+                                ,std::vector<int>tdc
+                                ,int sample=2);
+
+
+  /// QIE TDC Current threshold
+  float tdc_thr_;
+
+  /// QIE Sampling frequency (in MHz)
+  float qie_sf_{40.};
+
+  /// Class to set the verbosity level.
+  // TODO: Make use of the global verbose parameter.
+  bool verbose_{false};
+
+  /// Name of the input collection containing the sim hits
+  std::string inputCollection_;
+
+  /// Name of the pass that the input collection is on (empty string means take
+  /// any pass)
+  std::string inputPassName_;
+
+  /// Name of the output collection that will be used to stored the
+  /// digitized trigger scintillator hits
+  std::string outputCollection_;
+
+  /// SiPM gain
+  double gain_{1e6};
+
+  /// QIE pedestal
+  double pedestal_{6.0};
+
+  /// QIE pedestal
+  double noise_{1.5};
+
+  /// Total MeV per MIP
+  double mevPerMip_{1.40};
+
+  /// Total number of photoelectrons per MIP
+  double pePerMip_{13.5};
+
+  /// Sample of interest
+  int sample_of_interest_{2};
+
+  /// Input pulse shape for fitting
+  std::string input_pulse_shape_;
+
+  /// Input pulse parameters for fitting
+  std::vector<float> pulse_params_;
+};
+
+}  // namespace trigscint
+
+#endif
\ No newline at end of file
diff --git a/TrigScint/include/TrigScint/NumericalRecHitProducer.h b/TrigScint/include/TrigScint/NumericalRecHitProducer.h
new file mode 100644
index 000000000..1e7ae8f94
--- /dev/null
+++ b/TrigScint/include/TrigScint/NumericalRecHitProducer.h
@@ -0,0 +1,130 @@
+/**
+ * @file NumericalRecHitProducer.h
+ * @brief Class that builds recHits
+ * @author Andrew Whitbeck, TTU
+ */
+
+#ifndef TRIGSCINT_TRIGSCINTDIGIPRODUCER_H
+#define TRIGSCINT_TRIGSCINTDIGIPRODUCER_H
+
+/*~~~~~~~~~~*/
+/*   ROOT   */
+/*~~~~~~~~~~*/
+#include "TRandom3.h"
+#include "TVectorD.h"
+
+// LDMX
+#include "DetDescr/TrigScintID.h"
+#include "Recon/Event/EventConstants.h"
+#include "Tools/NoiseGenerator.h"
+#include "TrigScint/Event/TrigScintHit.h"
+#include "TrigScint/Event/TrigScintQIEDigis.h"
+
+/*~~~~~~~~~~~~~~~*/
+/*   Framework   */
+/*~~~~~~~~~~~~~~~*/
+#include "Framework/Configure/Parameters.h"
+#include "Framework/EventProcessor.h"
+
+/*~~~~~~~~~~~*/
+/* TrigScint */
+/*~~~~~~~~~~~*/
+#include "TrigScint/SimQIE.h"
+
+namespace trigscint {
+
+/**
+ * @class NumericalRecHitProducer
+ * @brief Organizes digis into TrigScintHits, linearizes TDC
+ * and ADC info, and converts amplitudes to PEs
+
+ */
+class NumericalRecHitProducer : public framework::Producer {
+ public:
+  NumericalRecHitProducer(const std::string& name, framework::Process& process);
+
+  ~NumericalRecHitProducer();
+
+  /**
+   * Callback for the processor to configure itself from the given set
+   * of parameters.
+   *
+   * @param parameters ParameterSet for configuration.
+   */
+  void configure(framework::config::Parameters& parameters) final override;
+
+  void produce(framework::Event& event);
+
+  /**
+   * Const function for pulse fitting
+   * @param params an array of 2 elements specifying
+   * pulse arrival time and pulse amplitude (Total integral)
+   */
+  double CostFunction(const double* params);
+
+  /// QIE Sampling frequency (in MHz)
+  float qie_sf_{40.};
+
+ private:
+  /**
+   * Reconstruct true charge deposited in each time sample
+   * @param adc array of adcs for give event, cell
+   * @param tdc array of tdcs for give event, cell
+   * @param sample sample of interest
+   */
+  Double_t ChargeReconstruction(std::vector<int>adc
+                                ,std::vector<int>tdc
+                                ,int sample=2);
+
+  /// Linearized charge. (Will be updated every time sample)
+  double Qm{0};
+  
+  /// Time of crossing tdc threshold (Will be updated every time sample)
+  double tm{0};
+
+  /// QIE TDC Current threshold
+  float tdc_thr_;
+
+  /// Class to set the verbosity level.
+  // TODO: Make use of the global verbose parameter.
+  bool verbose_{false};
+
+  /// Name of the input collection containing the sim hits
+  std::string inputCollection_;
+
+  /// Name of the pass that the input collection is on (empty string means take
+  /// any pass)
+  std::string inputPassName_;
+
+  /// Name of the output collection that will be used to stored the
+  /// digitized trigger scintillator hits
+  std::string outputCollection_;
+
+  /// SiPM gain
+  double gain_{1e6};
+
+  /// QIE pedestal
+  double pedestal_{6.0};
+
+  /// QIE pedestal
+  double noise_{1.5};
+
+  /// Total MeV per MIP
+  double mevPerMip_{1.40};
+
+  /// Total number of photoelectrons per MIP
+  double pePerMip_{13.5};
+
+  /// Sample of interest
+  int sample_of_interest_{2};
+
+  /// Input pulse shape for fitting
+  std::string input_pulse_shape_;
+
+  /// Input pulse parameters for fitting
+  std::vector<float> pulse_params_;
+};
+
+}  // namespace trigscint
+
+#endif
\ No newline at end of file
diff --git a/TrigScint/python/trigScint.py b/TrigScint/python/trigScint.py
index f8fcbec77..b6ea1b512 100644
--- a/TrigScint/python/trigScint.py
+++ b/TrigScint/python/trigScint.py
@@ -218,6 +218,53 @@ def pad3() :
         rechit.input_collection  = 'trigScintQIEDigisPad3'
         rechit.output_collection = 'trigScintRecHitsPad3'
         return rechit
+    
+
+class AnalyticalRecHitProducer(ldmxcfg.Producer) :
+    """Configuration for rechit producer for Trigger Scintillators"""
+
+    def __init__(self,name) :
+        super().__init__(name,'trigscint::AnalyticalRecHitProducer','TrigScint')
+
+        self .mev_per_mip = 0.4   #\
+                                  # >>>both are for converting edep to PEs 
+        self.pe_per_mip = 100.    #/
+        self.pedestal= 6.0        # QIE pedestal value (in fC)
+        self.elec_noise = 1.5     # QIE Electronic noise (in fC)
+        self.gain = 1.e6      # SiPM Gain
+        self.input_collection="trigScintQIEDigisUp"
+        self.input_pass_name=""   #take any pass
+        self.output_collection="trigScintRecHitsUp"
+        self.verbose = False
+        self.sample_of_interest=2 # Sample of interest. Range 0 to 3
+
+        self.input_pulse_shape="Expo" # Name of the input pulse class
+        self.expo_k=0.1          # Inverse of decay time of piece-wise exponential
+        self.expo_tmax=5.0       # Time at which piece-wise exponential peaks
+        self.tdc_thr = 3.4       # Threshold current in uA for TDC latch (as used by QIE)
+        self.qie_sf = 40.        # QIE sampling frequency in MHz
+
+    def pad1() : 
+        """Get the rechit producer for upstream pad"""
+        rechit = AnalyticalRecHitProducer( 'trigScintRecHitsPad1' )
+        rechit.input_collection  = 'trigScintQIEDigisPad1'
+        rechit.output_collection = 'trigScintRecHitsPad1'
+        return rechit
+
+    def pad2() : 
+        """Get the rechit producer for downstream pad"""
+        rechit = AnalyticalRecHitProducer( 'trigScintRecHitsPad2' )
+        rechit.input_collection  = 'trigScintQIEDigisPad2'
+        rechit.output_collection = 'trigScintRecHitsPad2'
+        return rechit
+
+    def pad3() :
+        """Get the rechit producer for tagger pad"""
+        rechit = AnalyticalRecHitProducer( 'trigScintRecHitsPad3' )
+        rechit.input_collection  = 'trigScintQIEDigisPad3'
+        rechit.output_collection = 'trigScintRecHitsPad3'
+        return rechit
+
 
 class TrigScintClusterProducer(ldmxcfg.Producer) :
     """Configuration for cluster producer for Trigger Scintillators"""
diff --git a/TrigScint/src/TrigScint/AnalyticalRecHitProducer.cxx b/TrigScint/src/TrigScint/AnalyticalRecHitProducer.cxx
new file mode 100644
index 000000000..10fe033b0
--- /dev/null
+++ b/TrigScint/src/TrigScint/AnalyticalRecHitProducer.cxx
@@ -0,0 +1,203 @@
+#include "TrigScint/AnalyticalRecHitProducer.h"
+
+#include <iomanip>
+#include <iostream>
+
+#include "Framework/Exception/Exception.h"
+#include "Framework/RandomNumberSeedService.h"
+#include "TLinearFitter.h"
+#include "TMath.h"
+
+namespace trigscint {
+
+AnalyticalRecHitProducer::AnalyticalRecHitProducer(const std::string &name,
+                                                   framework::Process &process)
+    : Producer(name, process) {}
+
+AnalyticalRecHitProducer::~AnalyticalRecHitProducer() {}
+
+void AnalyticalRecHitProducer::configure(
+    framework::config::Parameters &parameters) {
+  // Configure this instance of the producer
+  pedestal_ = parameters.getParameter<double>("pedestal");
+  noise_ = parameters.getParameter<double>("elec_noise");
+  gain_ = parameters.getParameter<double>("gain");
+  mevPerMip_ = parameters.getParameter<double>("mev_per_mip");
+  pePerMip_ = parameters.getParameter<double>("pe_per_mip");
+  inputCollection_ = parameters.getParameter<std::string>("input_collection");
+  inputPassName_ = parameters.getParameter<std::string>("input_pass_name");
+  outputCollection_ = parameters.getParameter<std::string>("output_collection");
+  verbose_ = parameters.getParameter<bool>("verbose");
+  sample_of_interest_ = parameters.getParameter<int>("sample_of_interest");
+  tdc_thr_ = parameters.getParameter<double>("tdc_thr");
+  qie_sf_ = parameters.getParameter<double>("qie_sf");
+
+  input_pulse_shape_ =
+      parameters.getParameter<std::string>("input_pulse_shape");
+  if (input_pulse_shape_ == "Expo") {
+    pulse_params_.clear();
+    pulse_params_.push_back(parameters.getParameter<double>("expo_k"));
+    pulse_params_.push_back(parameters.getParameter<double>("expo_tmax"));
+
+    ldmx_log(debug) << "expo_k =" << pulse_params_[0];
+    ldmx_log(debug) << "expo_tmax =" << pulse_params_[1];
+  }
+}
+
+void AnalyticalRecHitProducer::produce(framework::Event &event) {
+  // initialize QIE object for linearizing ADCs
+  SimQIE qie;
+
+  // Ensure the sample of interest <4
+  if (sample_of_interest_ > 3) {
+    ldmx_log(error) << "sample_of_interest_ should be one of 0,1,2,3\n"
+                    << "Currently, sample_of_interest = " << sample_of_interest_
+                    << "\n";
+    return;
+  }
+
+  // looper over sim hits and aggregate energy depositions
+  // for each detID
+  const auto digis{event.getCollection<trigscint::TrigScintQIEDigis>(
+      inputCollection_, inputPassName_)};
+
+  std::vector<ldmx::TrigScintHit> trigScintHits;
+  for (const auto &digi : digis) {
+    ldmx::TrigScintHit hit;
+    auto adc{digi.getADC()};
+    auto tdc{digi.getTDC()};
+
+    hit.setModuleID(0);
+    hit.setBarID(digi.getChanID());
+    hit.setBeamEfrac(-1.);
+
+    if (tdc[sample_of_interest_] > 49)
+      hit.setTime(-999.);
+    else
+      hit.setTime(tdc[sample_of_interest_] * 0.5);
+
+    auto Charge_time = ChargeReconstruction(adc, tdc, sample_of_interest_);
+    auto Charge = Charge_time[0];
+    auto time_est = Charge_time[1];
+
+    hit.setAmplitude(Charge);
+    hit.setEnergy(Charge * 6250. / gain_ * mevPerMip_ / pePerMip_);  // MeV
+    hit.setPE(Charge * 6250. / gain_);
+    if (time_est != -1) {hit.setTime(time_est);}
+    if (hit.getTime() > 0.) {
+      //std::cout << "adding hit " << digi.getChanID() << " " << time_est << " " <<  hit.getEnergy() << std::endl;
+      trigScintHits.push_back(hit);
+    }
+  }
+  // Create the container to hold the
+  // digitized trigger scintillator hits.
+
+  event.add(outputCollection_, trigScintHits);
+}
+std::vector<double>  AnalyticalRecHitProducer::ChargeReconstruction(std::vector<int> adc,
+                                                        std::vector<int> tdc,
+                                                        int sample) {
+  int npulses = 0;               // No. of true pulses
+  int poi = -1;                   // The pulse of interest
+  std::vector<float> Charge_;    // stores pulse amplitudes
+  std::vector<float> time_;    // stores pulse times
+  auto Qdata = new Double_t[5];  // Linearized charge
+  float tend = 1000 / qie_sf_;   // 1 time sample (in ns)
+  float k_ = pulse_params_[0];
+  float tmax_ = pulse_params_[1];
+  float par0 = (exp(k_ * tmax_) - 1) / (k_ * tmax_) * exp(-k_ * tend);
+  float alpha[2] = {par0,par0*(1-exp(-k_ * tend))};
+  SimQIE qie;
+  auto pulse = new Expo(k_, tmax_);
+
+  for(int i=0;i<adc.size();i++){
+    // Linearize charge, remove pedestal
+    Qdata[i] = qie.ADC2Q(adc[i])-pedestal_;
+
+    // Remove noise (2sigma noise cut)
+    if(Qdata[i]<2*noise_)
+      Qdata[i]=0;
+  }
+  
+  for(int i=0;i<adc.size();i++){
+
+    if(tdc[i]<40) {
+      if (i == sample) poi = npulses;
+      auto tm = tdc[i] / 2;  // measured time
+      auto Qm = Qdata[i]-pulse->Integrate(tend * i, tend * (i + 1));  // remaining measured charge
+      //float Qreco = (par0 * tdc_thr_ * tmax_ * exp(k_ * tm)-Qdata[i])/(par0 * exp(k_ * tm) - 1);
+      float Qreco = (par0 * tdc_thr_ * tmax_ * exp(k_ * tm)-Qm)/(par0 * exp(k_ * tm) - 1);
+      Charge_.push_back(Qreco);
+      npulses++;
+      pulse->AddPulse(tend * i + tdc[i] / 2, Qreco);
+      time_.push_back(tdc[i] / 2);
+    }
+    else if (tdc[i] < 50) {
+      if (i == sample)
+	      poi = npulses;
+      //Charge_.push_back(qie.ADC2Q(adc[i]));
+      Charge_.push_back(Qdata[i]);
+      //pulse->AddPulse(tend * i + tdc[i] / 2, qie.ADC2Q(adc[i]));
+      pulse->AddPulse(tend * i + tdc[i] / 2, Qdata[i]);
+
+      npulses++;
+      time_.push_back(tdc[i] / 2);
+    }
+    else if(i<tdc.size()-1 && tdc[i+1]==62) {
+      if (i == sample)
+	poi = npulses;
+      float Q0 = Qdata[i]-pulse->Integrate(tend*i, tend*(i+1));
+      float Q1 = Qdata[i+1]-pulse->Integrate(tend*(i+1), tend*(i+2));
+      float t0 = log((Q1)/(alpha[1]*Q0+alpha[0]*Q1))/k_;
+      float Qreco = Q0/(1-alpha[0]*exp(k_*t0));      
+      Charge_.push_back(Qreco);
+      time_.push_back(t0);
+      //std::cout << "tdc 62 " << t0 << " " <<  Qreco << std::endl;
+      npulses++;
+      pulse->AddPulse(tend * i + tdc[i] / 2, Qreco);
+    }
+  }
+
+  /////////////// For Debigging purposes
+  if (verbose_) {
+    std::cout << outputCollection_ << std::endl;
+    std::cout << "TS \t|\t0\t|\t1\t|\t2\t|\t3\t|\t4\t|\n"
+              << "---------------------------------------------"
+              << "--------------------------------------------\n"
+              << "tdc \t|";
+    for (int i = 0; i < 5; i++) std::cout << std::setw(10) << tdc[i] << "\t|";
+
+    std::cout << "\nadc \t|";
+    for (int i = 0; i < 5; i++) std::cout << std::setw(10) << adc[i] << "\t|";
+
+    std::cout << "\nQdata\t|";
+    for (int i = 0; i < 5; i++) std::cout << std::setw(10) << Qdata[i] << "\t|";
+
+    std::cout << "\n---------------------------------------------"
+              << "--------------------------------------------";
+    for (int n = 0; n < npulses; n++) {
+      std::cout << std::setw(10) << "\nPulse" << n << "\t|";
+      for (int i = 0; i < 5; i++)
+        std::cout << std::setw(10) << pulse->Integrate(i * tend, (i + 1) * tend)
+                  << "\t|";
+    }
+
+    std::cout << "\n"
+              << "\nnpulses = " << npulses << std::endl
+              << "poi = " << poi << std::endl;
+
+    std::cout << "Charge_[]: ";
+    for (int i = 0; i < Charge_.size(); i++)
+      std::cout << std::setw(10) << Charge_[i] << "\t|";
+    std::cout << "\n";
+  }
+  std::vector<double> poi_vals;
+  if(poi == -1) {
+    poi_vals = {-1,-1};
+  }
+  else {poi_vals = {Charge_[poi], time_[poi]};}
+  return poi_vals; //Charge_[poi];
+}
+}  // namespace trigscint
+
+DECLARE_PRODUCER_NS(trigscint, AnalyticalRecHitProducer);
\ No newline at end of file
diff --git a/TrigScint/src/TrigScint/NumericalRecHitProducer.cxx b/TrigScint/src/TrigScint/NumericalRecHitProducer.cxx
new file mode 100644
index 000000000..7fcc47458
--- /dev/null
+++ b/TrigScint/src/TrigScint/NumericalRecHitProducer.cxx
@@ -0,0 +1,279 @@
+#include "TrigScint/NumericalRecHitProducer.h"
+#include "Framework/Exception/Exception.h"
+#include "Framework/RandomNumberSeedService.h"
+#include "TMath.h"
+
+#include "Math/Minimizer.h"
+#include "Math/Factory.h"
+#include "Math/Functor.h"
+
+#include <iostream>
+#include "TLinearFitter.h"
+#include<iomanip>
+
+namespace trigscint {
+
+  float tm_new=0;
+  std::vector<float> Qm_new;
+  float TimePeriod=0;		// time periodd of 1 TS
+  float k_,tmax_;
+  int npulses=0;
+  double NewCostFunction(const double* params);
+  
+NumericalRecHitProducer::NumericalRecHitProducer(const std::string &name,
+                                                 framework::Process &process)
+    : Producer(name, process) {}
+
+NumericalRecHitProducer::~NumericalRecHitProducer() {}
+
+void NumericalRecHitProducer::configure(
+					framework::config::Parameters &parameters) {
+  // Configure this instance of the producer
+  pedestal_ = parameters.getParameter<double>("pedestal");
+  noise_ = parameters.getParameter<double>("elec_noise");
+  gain_ = parameters.getParameter<double>("gain");
+  mevPerMip_ = parameters.getParameter<double>("mev_per_mip");
+  pePerMip_ = parameters.getParameter<double>("pe_per_mip");
+  inputCollection_ = parameters.getParameter<std::string>("input_collection");
+  inputPassName_ = parameters.getParameter<std::string>("input_pass_name");
+  outputCollection_ = parameters.getParameter<std::string>("output_collection");
+  verbose_ = parameters.getParameter<bool>("verbose");
+  sample_of_interest_ = parameters.getParameter<int>("sample_of_interest");
+  tdc_thr_ = parameters.getParameter<double>("tdc_thr");
+  qie_sf_ = parameters.getParameter<double>("qie_sf");
+
+  input_pulse_shape_ =
+      parameters.getParameter<std::string>("input_pulse_shape");
+  if (input_pulse_shape_ == "Expo") {
+    pulse_params_.clear();
+    pulse_params_.push_back(parameters.getParameter<double>("expo_k"));
+    pulse_params_.push_back(parameters.getParameter<double>("expo_tmax"));
+
+    k_ = pulse_params_[0];
+    tmax_ = pulse_params_[1];
+
+    ldmx_log(debug) << "expo_k =" << pulse_params_[0];
+    ldmx_log(debug) << "expo_tmax =" << pulse_params_[1];
+  }
+
+}
+
+  void NumericalRecHitProducer::produce(framework::Event &event) {
+  // initialize QIE object for linearizing ADCs
+  SimQIE qie;
+
+  // Ensure the sample of interest <4
+  if(sample_of_interest_>3) {
+    ldmx_log(error)<<"sample_of_interest_ should be one of 0,1,2,3\n"
+		   <<"Currently, sample_of_interest = "<<sample_of_interest_
+		   <<"\n";
+    return;
+  }
+
+  // looper over sim hits and aggregate energy depositions
+  // for each detID
+  const auto digis{event.getCollection<trigscint::TrigScintQIEDigis>(
+      inputCollection_, inputPassName_)};
+
+  std::vector<ldmx::TrigScintHit> trigScintHits;
+  for (const auto &digi : digis) {
+    ldmx::TrigScintHit hit;
+    auto adc{digi.getADC()};
+    auto tdc{digi.getTDC()};
+
+    hit.setModuleID(0);
+    hit.setBarID(digi.getChanID());
+    hit.setBeamEfrac(-1.);
+
+
+    if (tdc[sample_of_interest_] > 49)
+      hit.setTime(-999.);
+    else
+      hit.setTime(tdc[sample_of_interest_] * 0.5);
+
+    auto Charge =
+      ChargeReconstruction(adc,tdc,sample_of_interest_);
+
+    hit.setAmplitude(Charge);
+    hit.setEnergy(Charge * 6250. /
+		  gain_ * mevPerMip_ / pePerMip_);  // MeV
+    hit.setPE(Charge * 6250. /gain_);
+    trigScintHits.push_back(hit);
+  }
+  // Create the container to hold the
+  // digitized trigger scintillator hits.
+
+  event.add(outputCollection_, trigScintHits);
+}
+  Double_t NumericalRecHitProducer::ChargeReconstruction
+  (std::vector<int>adc,std::vector<int>tdc,int sample) {
+    
+    npulses = 0;                  // No. of true pulses
+    int poi=0;                    // The pulse of interest
+    std::vector<float> Charge_;	  // stores pulse amplitudes
+    std::vector<double> InitVal;  // Initial values of fitting parameters
+    auto Qdata = new Double_t[5]; // Linearized charge
+    float tend = 1000/qie_sf_;	  // 1 time sample (in ns)
+    TimePeriod = tend;
+    SimQIE qie;
+    auto pulse = new Expo(k_,tmax_);
+
+    // Initialize the minimizer.
+    ROOT::Math::Minimizer* min = 
+      ROOT::Math::Factory::CreateMinimizer("Minuit","Migrad");
+    
+    
+    // if(tdc[sample]==63) return(0); // no signal pulse
+
+    Qm_new.clear();
+    for(int i=0;i<tdc.size();i++) {
+      Qdata[i] = qie.ADC2Q(adc[i]);
+
+      // Clean up the measured charge
+      if(Qdata[i]<pedestal_+2*noise_)
+    	Qm_new.push_back(0);
+      else
+    	Qm_new.push_back(Qdata[i]-pedestal_);
+      
+      // Find the pulses to fit
+      if(tdc[i]<50 ||
+    	 (i<tdc.size()-1 && tdc[i+1]==62)) {
+      	if(i==sample_of_interest_)
+      	  poi=npulses;
+    	if(tdc[i]<50)
+    	  InitVal.push_back(tend*i+tdc[i]/2);
+    	else
+    	  InitVal.push_back(tend*i);
+    	InitVal.push_back(Qm_new[i]);
+    	npulses++;
+      }
+    }
+
+    if(npulses==0) {
+      if(verbose_) {
+	ldmx_log(error)<<"No pulse to fit. npulses = 0\n";
+	ldmx_log(error)<<"tdc "
+		       <<std::setw(10)<<tdc[0]
+		       <<std::setw(10)<<tdc[1]
+		       <<std::setw(10)<<tdc[2]
+		       <<std::setw(10)<<tdc[3]
+		       <<std::setw(10)<<tdc[4];
+	ldmx_log(error)<<"adc "
+		       <<std::setw(10)<<adc[0]
+		       <<std::setw(10)<<adc[1]
+		       <<std::setw(10)<<adc[2]
+		       <<std::setw(10)<<adc[3]
+		       <<std::setw(10)<<adc[4];
+	ldmx_log(error)<<"Qm_new "
+		       <<std::setw(10)<<Qm_new[0]
+		       <<std::setw(10)<<Qm_new[1]
+		       <<std::setw(10)<<Qm_new[2]
+		       <<std::setw(10)<<Qm_new[3]
+		       <<std::setw(10)<<Qm_new[4];
+	ldmx_log(error)<<"returning linear(adc["<<sample
+		       <<"]+adc["<<sample+1<<"])";
+      }
+      return(Qm_new[sample]+Qm_new[sample+1]);
+    }
+    
+    ROOT::Math::Functor f(&NewCostFunction,2*npulses);
+    min->SetFunction(f);
+
+    for(int n = 0;n<npulses;n++) {
+      char* id = new char[5];
+      sprintf(id,"t%d",n+1);
+      min->SetVariable(2*n,id,InitVal[2*n],0.1);
+      sprintf(id,"Q%d",n+1);
+      min->SetVariable(2*n+1,id,InitVal[2*n+1],1);
+    }
+    
+    min->SetMaxFunctionCalls(1000000); // for Minuit/Minuit2 
+    min->SetMaxIterations(10000);  // for GSL
+    min->SetTolerance(0.001);
+    // if(verbose_)
+    //   min->SetPrintLevel(1);
+
+    min->Minimize();
+    const double *pred = min->X();
+
+    // pulse->AddPulse(tend*i+pred[0],pred[1]);
+    // }
+
+    /////////////// For Debigging purposes
+    // if(verbose_ || pred[2*poi+1]<0) {
+    // if(verbose_ || pred[2*poi+1]*6250./gain_ * mevPerMip_ / pePerMip_> 1.6) {
+    if(verbose_){
+      std::cout<<"TS \t|\t0\t|\t1\t|\t2\t|\t3\t|\t4\t|\n"
+	       <<"---------------------------------------------"
+	       <<"--------------------------------------------\n"
+	       <<"tdc \t|";
+      for(int i=0;i<5;i++)
+	std::cout<<std::setw(10)<<tdc[i]<<"\t|";
+    
+      std::cout<<"\nadc \t|";
+      for(int i=0;i<5;i++)
+	std::cout<<std::setw(10)<<adc[i]<<"\t|";
+    
+      std::cout<<"\nQdata\t|";
+      for(int i=0;i<5;i++)
+	std::cout<<std::setw(10)<<Qdata[i]<<"\t|";
+
+      std::cout<<"\nQm_new\t|";
+      for(int i=0;i<5;i++)
+	std::cout<<std::setw(10)<<Qm_new[i]<<"\t|";
+
+      std::cout<<"\n---------------------------------------------"
+	       <<"--------------------------------------------";
+
+      for(int n = 0;n<npulses;n++){
+      	std::cout<<std::setw(10)<<"\nPulse"<<n<<"\t|";
+      	auto pls = new Expo(k_,tmax_);
+      	pls->AddPulse(pred[2*n],pred[2*n+1]);
+      	for(int i=0;i<5;i++)
+      	  std::cout<<std::setw(10)<<pls->Integrate(i*tend,(i+1)*tend)<<"\t|";
+      	std::cout<<" Q= "<<pred[2*n+1];
+      	std::cout<<" t0= "<<pred[2*n];
+      }
+      std::cout<<"\n"
+      	       <<"\nnpulses = "<<npulses<<std::endl
+      	       <<"poi = "<<poi<<std::endl;
+      std::cout<<"\n\n";
+    }
+
+    return pred[2*poi+1];
+  }
+
+  double NewCostFunction(const double* params) {
+    double cost=0;
+    auto pulse = new Expo(k_,tmax_);
+    for(int n=0;n<npulses;n++)
+      pulse->AddPulse(params[2*n],params[2*n+1]);
+
+    // std::cout<<"\nQpred\t|";
+    for(int i=0;i<5;i++){
+      double Qpred = pulse->Integrate(i*TimePeriod,(i+1)*TimePeriod);
+      cost += pow(Qm_new[i]-Qpred,2);
+      // std::cout<<std::setw(10)<<Qpred<<"\t|";
+    }
+
+    // std::cout<<"params: "<<std::setw(10)<<params[0]<<"\t|"<<std::setw(10)<<params[1]<<"\t|";
+    // std::cout<<"cost: "<<std::setw(10)<<cost<<"\n";
+    return cost;
+  }
+
+  double NumericalRecHitProducer::CostFunction(const double* params) {
+    SimQIE qie;
+    auto pulse = new Expo(pulse_params_[0],pulse_params_[1]);
+    pulse->AddPulse(params[0],params[1]);
+    pulse->AddPulse(params[2],params[3]);
+
+    double Qpred = pulse->Integrate(0,1000/qie_sf_);
+    double Tpred = qie.TDC(pulse,0)/2;
+    
+    double cost1 = pow(tm-Tpred,2);
+    double cost2 = pow(Qm-Qpred,2);
+    return(cost1+cost2);
+  }
+}  // namespace trigscint
+
+DECLARE_PRODUCER_NS(trigscint, NumericalRecHitProducer);
\ No newline at end of file