mcsmear: write the FDC strip peakAmp and carry it through the random-trigger merge

src/programs/Simulation/mcsmear/FDCSmearer.cc

@@ -37,6 +37,15 @@ fdc_config_t::fdc_config_t(const std::shared_ptr<const JEvent>& event)
        	FDC_THRESH_KEV 		  = fdcparms["FDC_THRESH_KEV"]; 
 	}
 
+	// scale for the strip peak amplitude written in fdcDigihit (0 = do not write fdcDigihit)
+	FDC_ADC_ASCALE = 0.;
+	map<string, double> digi_scales;
+	if(DEvent::GetCalib(event, "FDC/digi_scales", digi_scales) || digi_scales.find("FDC_ADC_ASCALE") == digi_scales.end())
+		jerr << "Problem loading FDC/digi_scales FDC_ADC_ASCALE from CCDB! fdcDigihit peakAmp will not be written." << endl;
+	else
+		FDC_ADC_ASCALE = digi_scales["FDC_ADC_ASCALE"];
+	jout << "FDC strip peakAmp = q / FDC_ADC_ASCALE, FDC_ADC_ASCALE = " << FDC_ADC_ASCALE << endl;
+
    	// Calculate ped noise level based on position resolution
    	//   FDC_PED_NOISE = -0.004594 + 0.008711*FDC_CATHODE_SIGMA +
    	//                    0.000010*FDC_CATHODE_SIGMA*FDC_CATHODE_SIGMA; //pC
@@ -113,6 +122,11 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record)
                hddm_s::FdcCathodeHitList hits = siter->addFdcCathodeHits();
                hits().setQ(q);
                hits().setT(t);
+               // strip peak amplitude on the data scale, read by DEventSourceHDDM as DFDCHit::pulse_height
+               if (fdc_config->FDC_ADC_ASCALE > 0) {
+                  hddm_s::FdcDigihitList digihit = hits().addFdcDigihits();
+                  digihit().setPeakAmp(q / fdc_config->FDC_ADC_ASCALE);
+               }
             }
          }
 

src/programs/Simulation/mcsmear/FDCSmearer.h

@@ -18,6 +18,7 @@ class fdc_config_t
 	//double FDC_HIT_DROP_FRACTION = 0.0; // 1000.0E-9
 	double FDC_TIME_WINDOW;
 	double FDC_THRESH_KEV;        // fdc anode discriminator threshold
+	double FDC_ADC_ASCALE;        // /FDC/digi_scales FDC_ADC_ASCALE; strip peakAmp = q / FDC_ADC_ASCALE, as in data (DFDCHit_factory: q = a_scale * pulse_height)
 
     static bool FDC_EFFVSDOCA;
     double FDC_EFFVSDOCA_PAR[4];

src/programs/Simulation/mcsmear/hddm_s_merger.cc

@@ -960,16 +960,33 @@ hddm_s::FdcCathodeStripList &operator+=(hddm_s::FdcCathodeStripList &dst,
    return dst;
 }
 
+// Peak amplitude of an FDC strip hit, as the reconstruction reads it: the
+// fdcDigihit peakAmp if there is one, otherwise q. If the hit has no
+// fdcDigihit, one holding q is added so that the merged amplitude can be
+// stored in it.
+static double fdc_strip_peakamp(hddm_s::FdcCathodeHit &hit)
+{
+   if (hit.getFdcDigihits().size() > 0)
+      return hit.getFdcDigihit().getPeakAmp();
+   hddm_s::FdcDigihitList digihit = hit.addFdcDigihits();
+   digihit().setPeakAmp(hit.getQ());
+   return hit.getQ();
+}
+
 hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst,
                                       hddm_s::FdcCathodeHitList &src)
 {
+   if (src.size() == 0)
+      return dst;
 
    // There should be no more than one hit in dst - assume this is the case.  (otherwise, would have to pick the earliest)
 
    // Find the earliest (random) hit in src, ignore the later ones, if any. 
    // If it is before the original hit, replace the original hit's time and peak height. 
    // If it is at the same time sample as the original hit, sum the pulse height.
    // If it is after the original hit, leave the original hit's time and peak height alone. 
+   // The strip peak amplitude (fdcDigihit peakAmp, read by the reconstruction as the
+   // pulse height) follows the same rule as q.
 
 
    hddm_s::FdcCathodeHitList::iterator iter;
@@ -990,11 +1007,17 @@ hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst,
    if (t > fdc_strips_readout_window_ns) return dst; // random is too late
 
    double newQ = iter_firsthit->getQ();
+   bool newA_exists = (iter_firsthit->getFdcDigihits().size() > 0);
+   double newA = (newA_exists)? iter_firsthit->getFdcDigihit().getPeakAmp() : newQ;
 
    if (dst.size() == 0) { // no hits exist, add the random
      dst.add(1, -1);
      dst(0).setT(t);
      dst(0).setQ(newQ);
+     if (newA_exists) {
+        hddm_s::FdcDigihitList digihit = dst(0).addFdcDigihits();
+        digihit().setPeakAmp(newA);
+     }
 
    } else {
 
@@ -1009,13 +1032,17 @@ hddm_s::FdcCathodeHitList &operator+=(hddm_s::FdcCathodeHitList &dst,
 
       if (src_hitsample == dst_hitsample) { // same sample, add pulse heights
 
+         double origA = fdc_strip_peakamp(dst(0));
  	 if (t < origT) dst(0).setT(t);
 	 dst(0).setQ(origQ + newQ);
+         dst(0).getFdcDigihit().setPeakAmp(origA + newA);
 
       } else if (src_hitsample < dst_hitsample) { // random arrives earlier, replace hit time and pulse height 
 
+         fdc_strip_peakamp(dst(0));
          dst(0).setT(t);
 	 dst(0).setQ(newQ);           
+         dst(0).getFdcDigihit().setPeakAmp(newA);
 
       } // if the random hit is after the original one, it is too late, ignore it.