diff --git a/lalapps/src/pulsar/Fstatistic/synthesizeTransientStats.c b/lalapps/src/pulsar/Fstatistic/synthesizeTransientStats.c
index 7e1a8c647b12f1fa1e1391ed0468e77509ecc702..8cd3420d323ed565a3f210c86946702fa640b012 100644
--- a/lalapps/src/pulsar/Fstatistic/synthesizeTransientStats.c
+++ b/lalapps/src/pulsar/Fstatistic/synthesizeTransientStats.c
@@ -71,6 +71,14 @@
 #define CUBE(x) ((x)*(x)*(x))
 #define QUAD(x) ((x)*(x)*(x)*(x))
 
+/* define min/max macros if not already defined */
+#ifndef min
+#define min(a,b) ((a)<(b)?(a):(b))
+#endif
+#ifndef max
+#define max(a,b) ((a)>(b)?(a):(b))
+#endif
+
 /* ---------- local types ---------- */
 
 /** User-variables: can be set from config-file or command-line */
@@ -106,7 +114,8 @@ typedef struct {
   INT4  searchWindow_dtau;       /**< Step-size for search/marginalization over transient-window timescale, in seconds [Default:Tsft] */
 
   /* other parameters */
-  CHAR *IFO;		/**< IFO name */
+  CHAR *IFO;		/**< IFO name [deprecated] */
+  LALStringVector* IFOs; /**< list of detector-names "H1,H2,L1,.." or single detector*/
   INT4 dataStartGPS;	/**< data start-time in GPS seconds */
   INT4 dataDuration;	/**< data-span to generate */
   LALStringVector *timestampsFiles;	/**< per-detector file(s) with timestamps list */
@@ -508,6 +517,10 @@ XLALInitUserVars ( UserInput_t *uvar )
 #define DEFAULT_IFO "H1"
   uvar->IFO = XLALMalloc ( strlen(DEFAULT_IFO)+1 );
   strcpy ( uvar->IFO, DEFAULT_IFO );
+  if ( (uvar->IFOs = XLALCreateStringVector ( "H1", NULL )) == NULL ) {
+    LogPrintf (LOG_CRITICAL, "Call to XLALCreateStringVector() failed with xlalErrno = %d\n", xlalErrno );
+    XLAL_ERROR ( XLAL_ENOMEM );
+  }
 
   /* ---------- transient window defaults ---------- */
 #define DEFAULT_TRANSIENT "rect"
@@ -550,7 +563,8 @@ XLALInitUserVars ( UserInput_t *uvar )
 
   XLALRegisterUvarMember( AmpPriorType,	 	 INT4, 0,  OPTIONAL, "Enumeration of types of amplitude-priors: 0=physical, 1=canonical");
 
-  XLALRegisterUvarMember( IFO,	        STRING, 'I', OPTIONAL, "Detector: 'G1','L1','H1,'H2', 'V1', ... ");
+  XLALRegisterUvarMember( IFO,	        STRING, 0, DEPRECATED, "Detector: 'G1','L1','H1,'H2', 'V1', ... (use --IFOs instead)");
+  XLALRegisterUvarMember( IFOs,	        STRINGVector, 'I', OPTIONAL, "Comma-separated list of detectors, e.g. \"H1,H2,L1,G1, ...\" ");
   XLALRegisterUvarMember( dataStartGPS,	 	 INT4, 0,  OPTIONAL, "data start-time in GPS seconds");
   XLALRegisterUvarMember( dataDuration,	 	 INT4, 0,  OPTIONAL, "data-span to generate (in seconds)");
   XLALRegisterUvarMember( timestampsFiles,	        STRINGVector, 0,  OPTIONAL, "ALTERNATIVE: file(s) to read timestamps from (file-format: lines with <seconds> [<nanoseconds>]; nanoseconds currently ignored; only 1 detector/file currently supported)");
@@ -664,15 +678,16 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
   }
 
   /* init detector info */
-  LALDetector *site;
-  if ( (site = XLALGetSiteInfo ( uvar->IFO )) == NULL ) {
-    XLALPrintError ("%s: Failed to get site-info for detector '%s'\n", __func__, uvar->IFO );
-    XLAL_ERROR ( XLAL_EFUNC );
+  if ( XLALUserVarWasSet ( &uvar->IFO ) && XLALUserVarWasSet ( &uvar->IFOs ) ) {
+    XLALPrintError ("%s: Please do not use both the deprecated --IFO and the newer --IFOs option.\n", __func__ );
+    XLAL_ERROR ( XLAL_EINVAL );
+  }
+  else if ( XLALUserVarWasSet ( &uvar->IFO ) ) {
+    strcpy(uvar->IFOs->data[0], uvar->IFO);
   }
+  UINT4 numDetectors = uvar->IFOs->length;
   MultiLALDetector multiDet;
-  multiDet.length = 1;
-  multiDet.sites[0] = (*site); 	/* copy! */
-  XLALFree ( site );
+  XLAL_CHECK ( XLALParseMultiLALDetector ( &multiDet, uvar->IFOs ) == XLAL_SUCCESS, XLAL_EFUNC );
 
   /* TIMESTAMPS: either from --startTime+duration OR --timestampsFiles */
   BOOLEAN have_startTime = XLALUserVarWasSet ( &uvar->dataStartGPS );
@@ -685,41 +700,43 @@ XLALInitCode ( ConfigVariables *cfg, const UserInput_t *uvar )
   XLAL_CHECK ( !have_timestampsFiles || !have_startTime, XLAL_EINVAL, "--timestampsFiles incompatible with --dataStartGPS and --dataDuration}\n");
   MultiLIGOTimeGPSVector * multiTS;
   UINT4 numSteps = 0;
-  INT4 dataStartGPS = 0;
+  INT4 dataStartGPS = LAL_INT4_MAX;
+  INT4 dataEndGPS = 0;
   INT4 dataDuration = 0;
   /* now handle the two mutually-exclusive cases */
   if ( have_timestampsFiles )
     {
       /* NOTE: nanoseconds will be truncated! */
       XLAL_CHECK ( (multiTS = XLALReadMultiTimestampsFiles ( uvar->timestampsFiles )) != NULL, XLAL_EFUNC );
-      XLAL_CHECK ( (multiTS->length > 0) && (multiTS->data != NULL), XLAL_EINVAL, "Got empty timestamps-list from XLALReadMultiTimestampsFiles()\n" );
-      XLAL_CHECK ( (multiTS->length == 1) , XLAL_EINVAL, "Can currently deal only with one --timestampsFiles entry for a single detector!\n" );
-      XLAL_CHECK ( (multiTS->data[0]->length > 0) && (multiTS->data[0]->data != NULL), XLAL_EINVAL, "Got empty timestamps-list for detector %s", uvar->IFO );
-      numSteps = multiTS->data[0]->length;
-      multiTS->data[0]->deltaT = uvar->TAtom; // Tsft information not given by timestamps-file
-      dataStartGPS = multiTS->data[0]->data[0].gpsSeconds;
-      dataDuration = uvar->TAtom + multiTS->data[0]->data[numSteps-1].gpsSeconds - dataStartGPS; // difference of last and first timestamp, plus one extra atom length
+      XLAL_CHECK ( (multiTS->length > 0) && (multiTS->data != NULL), XLAL_EINVAL, "Got empty timestamps-list from XLALReadMultiTimestampsFiles()" );
+      XLAL_CHECK ( (multiTS->length == numDetectors) , XLAL_EINVAL, "Number of timestamps files does not match number of IFOs! (%d!=%d)", multiTS->length, numDetectors );
+      for ( UINT4 X=0; X < numDetectors; X++ ) {
+          XLAL_CHECK ( (multiTS->data[X]->length > 0) && (multiTS->data[X]->data != NULL), XLAL_EINVAL, "Got empty timestamps-list for detector %s", uvar->IFOs->data[X] );
+          numSteps = multiTS->data[X]->length;
+          multiTS->data[X]->deltaT = uvar->TAtom; // Tsft information not given by timestamps-file
+          dataStartGPS = min(dataStartGPS, multiTS->data[X]->data[0].gpsSeconds);
+          dataEndGPS = max(dataEndGPS, multiTS->data[X]->data[numSteps-1].gpsSeconds);
+      }
+      dataDuration = uvar->TAtom + dataEndGPS - dataStartGPS; // difference of last and first timestamp, plus one extra atom length
     } // endif have_timestampsFiles
   else
     {
       /* init timestamps vector covering observation time */
-      XLAL_CHECK ( (multiTS = XLALCalloc ( 1, sizeof(*multiTS))) != NULL, XLAL_ENOMEM );
-      multiTS->length = 1;
-      XLAL_CHECK ( (multiTS->data = XLALCalloc (1, sizeof(*multiTS->data))) != NULL, XLAL_ENOMEM );
+      XLAL_CHECK ( (multiTS = XLALCreateMultiLIGOTimeGPSVector (numDetectors)) != NULL, XLAL_EINVAL, "XLALCreateMultiLIGOTimeGPSVector(%d) failed.", numDetectors );
       dataStartGPS = uvar->dataStartGPS;
       dataDuration = uvar->dataDuration;
       numSteps = (UINT4) ceil ( dataDuration / uvar->TAtom );
-      XLAL_CHECK ( (multiTS->data[0] = XLALCreateTimestampVector (numSteps)) != NULL, XLAL_EFUNC, "XLALCreateTimestampVector(%d) failed.", numSteps );
-      multiTS->data[0]->deltaT = uvar->TAtom;
-      UINT4 i;
-      for ( i=0; i < numSteps; i ++ )
-        {
-          UINT4 ti = dataStartGPS + i * uvar->TAtom;
-          multiTS->data[0]->data[i].gpsSeconds = ti;
-          multiTS->data[0]->data[i].gpsNanoSeconds = 0;
-        }
+      XLAL_CHECK ( numSteps>=2, XLAL_EINVAL, "Need timestamps covering at least 2 atoms!" );
+      for ( UINT4 X=0; X < numDetectors; X++ ) {
+          XLAL_CHECK ( (multiTS->data[X] = XLALCreateTimestampVector (numSteps)) != NULL, XLAL_EFUNC, "XLALCreateTimestampVector(%d) failed.", numSteps );
+          multiTS->data[X]->deltaT = uvar->TAtom;
+          for ( UINT4 i=0; i < numSteps; i ++ ) {
+              UINT4 ti = dataStartGPS + i * uvar->TAtom;
+              multiTS->data[X]->data[i].gpsSeconds = ti;
+              multiTS->data[X]->data[i].gpsNanoSeconds = 0;
+          }
+      }
     } // endif !have_noiseSFTs
-  XLAL_CHECK(numSteps>=2,XLAL_EINVAL,"Need timestamps covering at least 2 atoms!");
 
   /* get detector states */
   if ( (cfg->multiDetStates = XLALGetMultiDetectorStates ( multiTS, &multiDet, edat, 0.5 * uvar->TAtom )) == NULL ) {
diff --git a/lalapps/src/pulsar/Fstatistic/test_synthesizeTransientStats.sh b/lalapps/src/pulsar/Fstatistic/test_synthesizeTransientStats.sh
index 6491827da7183db99d274a0dd19386318e41efa8..08e4e305d71774b79bf922dcc0e243e08d8d4f9d 100644
--- a/lalapps/src/pulsar/Fstatistic/test_synthesizeTransientStats.sh
+++ b/lalapps/src/pulsar/Fstatistic/test_synthesizeTransientStats.sh
@@ -1,5 +1,8 @@
 ## functions
 
+awk_reldev='{printf "%.2e", sqrt(($1-$2)*($1-$2))/(0.5*($1+$2)) }'
+awk_isgtr='{if($1>$2) {print "1"}}'
+
 function check_average {
     local file="$1"
     local col="$2"
@@ -22,30 +25,195 @@ function check_average {
     awk "${awk_script}" "${file}"
 }
 
-## common argments
+function twoF_from_atoms(){
+    local file="$1"
+    local awk_script='
+      /^%/ { next }
+      {
+        A += $2;
+        B += $3;
+        C += $4;
+        Fa_re += $5;
+        Fa_im += $6;
+        Fb_re += $7;
+        Fb_im += $8;
+        count += 1
+      }
+      END {
+        Dinv = 1. / ( A * B - C * C );
+        twoF = 2. * Dinv * ( B * ( Fa_re * Fa_re + Fa_im * Fa_im ) + A * ( Fb_re * Fb_re + Fb_im * Fb_im ) - 2.0 * C * ( Fa_re * Fb_re + Fa_im * Fb_im ) );
+        print twoF
+      }
+    '
+    twoF=$(awk "${awk_script}" "${file}")
+    echo $twoF
+}
+
+## common arguments
 
 timestamp1=818845553
 duration=86400
 numDraws=1000
-common_args="--IFO=H1 --dataStartGPS=$timestamp1 --dataDuration=$duration --numDraws=$numDraws --randSeed=1"
+common_args="--IFOs=H1 --dataStartGPS=$timestamp1 --dataDuration=$duration --numDraws=$numDraws --randSeed=1"
 
 ## lalapps_synthesizeTransientStats with a single square window should give the same results as lalapps_synthesizeLVStats
 
-synthLV="lalapps_synthesizeLVStats ${common_args} --computeBSGL=false"
-synthTS="lalapps_synthesizeTransientStats ${common_args} --injectWindow-type=rect --injectWindow-tauDays=1 --injectWindow-tauDaysBand=0 --injectWindow-t0Days=0 --injectWindow-t0DaysBand=0 --searchWindow-type=rect --searchWindow-tauDays=1 --searchWindow-tauDaysBand=0 --searchWindow-t0Days=0 --searchWindow-t0DaysBand=0"
+synthLVcode="lalapps_synthesizeLVStats"
+synthTScode="lalapps_synthesizeTransientStats"
+
+transient_args="--injectWindow-type=rect --injectWindow-tauDays=1 --injectWindow-tauDaysBand=0 --injectWindow-t0Days=0 --injectWindow-t0DaysBand=0 --searchWindow-type=rect --searchWindow-tauDays=1 --searchWindow-tauDaysBand=0 --searchWindow-t0Days=0 --searchWindow-t0DaysBand=0"
+
+CLsynthLV="$synthLVcode ${common_args} --computeBSGL=false"
+CLsynthTS="$synthTScode ${common_args} ${transient_args}"
 
 # --- first try with Gaussian noise
 
-${synthLV} --fixedSNR=0 --outputStats="synthLV_stats_g.dat"
-${synthTS} --fixedSNR=0 --outputStats="synthTS_stats_g.dat"
+cmdline="${CLsynthLV} --fixedSNR=0 --outputStats=synthLV_stats_H1_g.dat"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthLVcode' ..."
+    exit 1
+fi
+cmdline="${CLsynthTS} --fixedSNR=0 --outputStats=synthTS_stats_H1_g.dat"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthTScode' ..."
+    exit 1
+fi
 
-check_average "synthLV_stats_g.dat" 5 "2F" 4.0
-check_average "synthTS_stats_g.dat" 9 "2F" 4.0
+check_average "synthLV_stats_H1_g.dat" 5 "2F" 4.0
+check_average "synthTS_stats_H1_g.dat" 9 "2F" 4.0
 
 # --- now try with signals
 
-${synthLV} --fixedSNR=4 --outputStats="synthLV_stats_s.dat"
-${synthTS} --fixedSNR=4 --outputStats="synthTS_stats_s.dat"
+cmdline="${CLsynthLV} --fixedSNR=4 --outputStats=synthLV_stats_H1_s.dat"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthLVcode' ..."
+    exit 1
+fi
+cmdline="${CLsynthTS} --fixedSNR=4 --outputStats=synthTS_stats_H1_s.dat"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthTScode' ..."
+    exit 1
+fi
+
+check_average "synthLV_stats_H1_s.dat" 5 "2F" 20.0
+check_average "synthTS_stats_H1_s.dat" 9 "2F" 20.0
+
+# --- now try multi-IFO mode
+
+common_args="--IFOs=H1,L1 --dataStartGPS=$timestamp1 --dataDuration=$duration --numDraws=$numDraws --randSeed=1"
+cmdline="$synthTScode ${common_args} ${transient_args} --fixedSNR=4 --outputStats=synthTS_stats_H1L1_s.dat"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthTScode' ..."
+    exit 1
+fi
+check_average "synthTS_stats_H1L1_s.dat" 9 "2F" 20.0
+
+# --- multi-IFO and timestampsfiles
+
+tsH1="ts_H1.txt"
+tsL1="ts_L1.txt"
+TSFT=1800
+Nsfts=48
+Tend=$(($timestamp1 + $duration - $TSFT)) # last SFT should *end* on first+duration
+for i in `seq $timestamp1 $TSFT $Tend`; do
+    echo "$i 0" >> ${tsH1}
+done
+for i in `seq $timestamp1 $TSFT $Tend`; do
+    echo "$(($i + 900)) 0" >> ${tsL1}
+done
+
+common_args="--IFOs=H1,L1 --timestampsFiles=$tsH1,$tsL1 --numDraws=$numDraws --randSeed=1"
+cmdline="$synthTScode ${common_args} ${transient_args} --fixedSNR=4 --outputStats=synthTS_stats_H1L1_s_TS.dat"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthTScode' ..."
+    exit 1
+fi
+check_average "synthTS_stats_H1L1_s_TS.dat" 9 "2F" 20.0
+
+# test single- and multi-IFO atoms on a single draw
+# very strong signal but avoiding SignalOnly with its custom +-4 handling
+statsH1="synthTS_stats_H1_s_N1.dat"
+statsH1L1="synthTS_stats_H1L1_s_N1.dat"
+atomsH1="synthTS_atoms_H1_s"
+atomsH1L1="synthTS_atoms_H1L1_s_TS"
+snr=4000
+common_args="--numDraws=1 --randSeed=1 --fixedSNR=$snr --computeFtotal"
+cmdline="$synthTScode ${common_args} ${transient_args} --IFOs=H1 --timestampsFiles=$tsH1 --outputStats=$statsH1 --outputAtoms=$atomsH1"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthTScode' ..."
+    exit 1
+fi
+cmdline="$synthTScode ${common_args} ${transient_args} --IFOs=H1,L1 --timestampsFiles=$tsH1,$tsL1 --outputStats=$statsH1L1 --outputAtoms=$atomsH1L1"
+echo $cmdline
+if ! eval "$cmdline"; then
+    echo "Error.. something failed when running '$synthTScode' ..."
+    exit 1
+fi
+
+echo "Checking number of timestamps in atoms output..."
+atomsH1="${atomsH1}_0001_of_0001.dat"
+atomsH1L1="${atomsH1L1}_0001_of_0001.dat"
+
+NatomsH1=$(grep -v % $atomsH1 | wc -l)
+NatomsH1L1=$(grep -v % $atomsH1L1 | wc -l)
+if [ $NatomsH1 -ne $Nsfts ]; then
+    echo "Error.. number of atoms $NatomsH1 in H1 output file does not match expectation $Nsfts ..."
+    exit 1
+fi
+Natoms2=$(echo $Nsfts | awk '{printf "%d", 2*$1}')
+if [ $NatomsH1L1 -ne $Natoms2 ]; then
+    echo "Error.. number of atoms $NatomsH1L1 in H1L1 output file does not match expectation 2*$Nsfts=$Natoms2 ..."
+    exit 1
+fi
+
+echo "...got $NatomsH1 atoms in H1 output file and $NatomsH1L1 in H1L1 output file, as expected."
+
+# compare the TOTAL F-stat against summing up the atoms ourselves
+# (in this example, we expect maxTwoF=twoF, but in general deciding how many atoms to sum up for maxTwoF would be more work)
+# and also compare recovered SNRs from sqrt(2F-4) against injected value
+twoF_H1=$(grep -v % $statsH1 | tr -s " " | sed "s/^ *//g" | cut -d " " -f 6)
+snr_H1=$(echo $twoF_H1 | awk '{printf "%.6f", sqrt($1-4)}')
+reldev0=$(echo $snr_H1 $snr | awk "$awk_reldev")
+echo "2F from stats file (H1 only): $twoF_H1 (snr=$snr_H1 off from injected $snr by $reldev0)"
+twoF_H1L1=$(grep -v % $statsH1L1 | tr -s " " | sed "s/^ *//g" | cut -d " " -f 6)
+snr_H1L1=$(echo $twoF_H1L1 | awk '{printf "%.6f", sqrt($1-4)}')
+reldev1=$(echo $snr_H1L1 $snr | awk "$awk_reldev")
+echo "2F from stats file (H1L1): $twoF_H1L1 (snr=$snr_H1L1 off from injected $snr by $reldev1)"
+
+twoF_H1_from_atoms=$(twoF_from_atoms ${atomsH1})
+reldev2=$(echo $twoF_H1_from_atoms $twoF_H1 | awk "$awk_reldev")
+echo "2F from atoms (H1 only): $twoF_H1_from_atoms (off from stats file by $reldev2)"
+twoF_H1L1_from_atoms=$(twoF_from_atoms ${atomsH1L1})
+reldev3=$(echo $twoF_H1L1_from_atoms $twoF_H1L1 | awk "$awk_reldev")
+echo "2F from atoms (H1L1): $twoF_H1L1_from_atoms (off from stats file by $reldev3)"
+
+grep -v % $atomsH1L1 | head -n 48 > "${atomsH1L1}_firsthalf"
+twoF_H1_from_H1L1_atoms=$(twoF_from_atoms "${atomsH1L1}_firsthalf")
+echo "2F from atoms (H1 only from first half of H1L1 file): $twoF_H1_from_H1L1_atoms"
+grep -v % $atomsH1L1 | tail -n 48 > "${atomsH1L1}_secondhalf"
+twoF_L1_from_H1L1_atoms=$(twoF_from_atoms "${atomsH1L1}_secondhalf")
+echo "2F from atoms (L1 only from second half of H1L1 file): $twoF_L1_from_H1L1_atoms"
+snr_summed=$(echo $twoF_H1_from_H1L1_atoms $twoF_L1_from_H1L1_atoms | awk '{printf "%.6f", sqrt($1-4+$2-4) }')
+reldev4=$(echo $snr_summed $snr | awk "$awk_reldev")
+echo "sqrt(sum(snr2)): $snr_summed (off from injected $snr by $reldev4)"
 
-check_average "synthLV_stats_s.dat" 5 "2F" 20.0
-check_average "synthTS_stats_s.dat" 9 "2F" 20.0
+tolerance=1e-3
+echo "Checking at tolerance $tolerance..."
+fail0=$(echo $reldev0 $tolerance | awk "$awk_isgtr")
+fail1=$(echo $reldev1 $tolerance | awk "$awk_isgtr")
+fail2=$(echo $reldev2 $tolerance | awk "$awk_isgtr")
+fail3=$(echo $reldev3 $tolerance | awk "$awk_isgtr")
+fail4=$(echo $reldev4 $tolerance | awk "$awk_isgtr")
+if [ "$fail0" -o "$fail1" -o "$fail2" -o "$fail3" -o "$fail4" ]; then
+    echo " ==> *FAILED*"
+    exit 1
+else
+    echo " ==> OK"
+fi