added analytical and numerical rec hit producers

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

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

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"""