mcsmear: FDC cathode strip readout thresholds on raw_amp = q/a_scale, per strip from CCDB

src/programs/Simulation/mcsmear/FDCSmearer.cc

@@ -8,6 +8,8 @@
 #include <sstream>
 
 bool fdc_config_t::FDC_EFFVSDOCA(1);
+double fdc_config_t::FDC_RAW_AMP_THR(-1.);
+const double fdc_config_t::FDC_RAW_AMP_THR_DEFAULT(50.);
 
 //-----------
 // fdc_config_t  (constructor)
@@ -37,6 +39,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
@@ -69,6 +80,30 @@ fdc_config_t::fdc_config_t(const std::shared_ptr<const JEvent>& event)
 		}
 	}
 
+	// readout thresholds of the cathode strips on raw_amp = q / FDC_ADC_ASCALE, in ADC counts, per strip from
+	// /FDC/packageN/strip_thresholds (12 rows, as strip_mc_efficiency); FDC_RAW_AMP_THR_DEFAULT where a table is absent
+	for(int package=1; package<=4; package++) {
+		vector< vector<double> > new_thresholds;
+		char ccdb_str[100];
+		sprintf(ccdb_str, "/FDC/package%d/strip_thresholds", package);
+		bool have_table = false;
+		if (FDC_RAW_AMP_THR < 0) {
+			have_table = !DEvent::GetCalib(event, ccdb_str, new_thresholds) && new_thresholds.size() == 12;
+			if (have_table)
+				jout << "FDC strip raw-amplitude thresholds of package " << package << " from " << ccdb_str << endl;
+			else
+				jout << "FDC strip raw-amplitude thresholds of package " << package << ": no " << ccdb_str
+				     << ", default " << FDC_RAW_AMP_THR_DEFAULT << " ADC counts" << endl;
+		}
+		for(int chamber=0; chamber<6; chamber++) {
+			channel_raw_amp_thresholds.push_back( have_table? new_thresholds[2*chamber+1] : vector<double>() );
+			channel_raw_amp_thresholds.push_back( vector<double>() );
+			channel_raw_amp_thresholds.push_back( have_table? new_thresholds[2*chamber] : vector<double>() );
+		}
+	}
+	if (FDC_RAW_AMP_THR >= 0)
+		jout << "FDC strip raw-amplitude threshold FDC:RAW_AMP_THRESHOLD = " << FDC_RAW_AMP_THR << " ADC counts for every strip" << endl;
+
     FDC_EFFVSDOCA_PAR[0] = 0.999;
     FDC_EFFVSDOCA_PAR[1] = 3.75e-4;
     FDC_EFFVSDOCA_PAR[2] = 0.506;
@@ -109,10 +144,19 @@ void FDCSmearer::SmearEvent(hddm_s::HDDM *record)
              	q += gDRandom.SampleGaussian(fdc_config->FDC_PED_NOISE);
              	t += gDRandom.SampleGaussian(fdc_config->FDC_TDRIFT_SIGMA)*1.0e9;
 			}
+            // readout threshold on the raw amplitude raw_amp = q / FDC_ADC_ASCALE, in ADC counts
+            if (fdc_config->FDC_ADC_ASCALE > 0
+                  && q / fdc_config->FDC_ADC_ASCALE < fdc_config->GetRawAmpThreshold(siter))
+               continue;
             if (q > threshold && t > config->TRIGGER_LOOKBACK_TIME && t < t_max) {
                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,11 +18,16 @@ 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;
+    static double FDC_RAW_AMP_THR;        // FDC:RAW_AMP_THRESHOLD in ADC counts; < 0 (default): per strip from CCDB /FDC/packageN/strip_thresholds,
+                                          // FDC_RAW_AMP_THR_DEFAULT where a table is absent; >= 0: this value for every strip (0: no cut)
+    static const double FDC_RAW_AMP_THR_DEFAULT;
     double FDC_EFFVSDOCA_PAR[4];
 
 	vector< vector<double > > channel_efficiencies;
+	vector< vector<double > > channel_raw_amp_thresholds;   // same row order as channel_efficiencies; empty rows for the anode wires
 
 	double GetEfficiencyCorrectionFactor(hddm_s::FdcCathodeStripList::iterator &siter) {
 		// cathode strips
@@ -37,6 +42,17 @@ class fdc_config_t
 		return channel_efficiencies.at(gPlane).at(element-1);
 	}
 
+	double GetRawAmpThreshold(hddm_s::FdcCathodeStripList::iterator &siter) {
+		// readout threshold of a cathode strip on raw_amp = q / FDC_ADC_ASCALE, in ADC counts
+		if (FDC_RAW_AMP_THR >= 0)
+			return FDC_RAW_AMP_THR;
+		int gPlane = 9*(siter->getModule()-1) + 3*(siter->getLayer()-1) 
+             		+ (siter->getPlane()-1);               // starts counting at 0
+		int element = siter->getStrip();
+		const vector<double> &row = channel_raw_amp_thresholds.at(gPlane);
+		return (element >= 1 && element <= (int)row.size())? row[element-1] : FDC_RAW_AMP_THR_DEFAULT;
+	}
+
 	double GetEfficiencyCorrectionFactor(hddm_s::FdcAnodeWireList::iterator &witer) {
 		// anode wires
 		int gPlane = 9*(witer->getModule()-1) + 3*(witer->getLayer()-1) + 1;  // starts counting at 0

src/programs/Simulation/mcsmear/MyProcessor.cc

@@ -162,6 +162,12 @@ void MyProcessor::Init(void)
    // By default, the empirical fdc DOCA-dependent efficiency
    // is applied to fdc wire hits in FDCSmearer. This can be
    // disabled by this commandline option.
+   // FDC cathode strip readout threshold on raw_amp = q / FDC_ADC_ASCALE (ADC counts). By default the
+   // per-strip values of /FDC/packageN/strip_thresholds are used, or 50 where a table is absent.
+   app->SetDefaultParameter("FDC:RAW_AMP_THRESHOLD", fdc_config_t::FDC_RAW_AMP_THR,
+                            "FDC strip dropped if raw_amp = q / FDC_ADC_ASCALE is below this, in ADC counts;"
+                            " < 0 = per strip from CCDB /FDC/packageN/strip_thresholds, 50 where absent (default);"
+                            " >= 0 = this value for every strip (0 = no cut)");
    app->SetDefaultParameter("FDC:USE_EFFVSDOCA", fdc_config_t::FDC_EFFVSDOCA,
                           "Turn on/off fdc doca-dependent wire hit efficiency,"
                           " \"0\"=no efficiency correction,"

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.