/*****************************************************************************
******************************************************************************
**<AUTO>
**
** FILE:	targets.cc
**
**<HTML>
**	This file contains C++ functions to select spectroscopic targets.
**</HTML>
**</AUTO>
**
**
** ENVIRONMENT:
**	C++.
**
** REQUIRED PRODUCTS:
**	photoSch
**
******************************************************************************
******************************************************************************
*/
#include <sdb.h>
#include <genericEnum.h>
#include <arFrameStat.h>
#include <arRun.h>
#include <arTrans.h>
#include <arPTrans.h>
#include <arPSCalib.h>
#include <arAstromCalib.h>
#include <arObject.h>
#include <arObjectData.h>
#include <arObjectCont.h>
#include <arObjectList.h>
#include <arFramesPipeRun.h>
#include <arFinalPhotoCalib.h>
#include <arSegment.h>
#include <arChunk.h>
#include <arTarget.h>
#include <dervish.h>
#include <taCalibObj.h>
#include <astrotools.h>
#include <arObjectStatus.h>
#include <taCosmic.h>
#include <atDustGetval.h>
#include "tsPM.h"

static const int nFilters = 5;  // Number of SDSS photometirc filters

/*****************************************************************************
******************************************************************************
**<AUTO EXTRACT>
**
** ROUTINE:	tsSegmentPMFetch
**
**<HTML>
**</HTML>
**</AUTO>
******************************************************************************
******************************************************************************
*/
RET_CODE
tsSegmentPMFetch(ooHandle(arSegment)& segment, /* IN: segment to export */
		 int field0Proc,               /* IN: first field to process
						*     (-1 defaults to first
						*     field in segment/run) */
		 int nFieldsProc,              /* IN: number of fields to
						*     process, starting with
						*     first field in the
						*     segment --- 0=all;
						*     used for debugging */
#ifdef SKIP
		 int split,		       /* IN: one HDU per field */
#endif
		 int scanline,                 /* IN: 1=OK_SCANLINE objects,
						*     0=PRIMARY objects */
                 CHAIN *objChain               /* OUT: chain to append to */)
#ifdef SKIP
		 const char *outputDir
		 /* IN: output directory if exporting.  If it contains
		  * the string %d, the %d is replaced with the camCol. */)
#endif
{
  // Get handle to Frames Pipeline Run for this segment
  ooHandle(arFramesPipeRun) run = segment->framesPipelineRun(oocRead);
  uint16 field0 = segment->field0();
  uint16 nFields = segment->nFields();
  uint16 segField0 = field0;
  int skyVersion = segment->chunk()->skyVersion();
  uint32 runId = run->run();
  uint8 camCol = run->camCol();
  uint16 rerun= run->rerun();
  if (field0Proc > field0)
    {
      nFields -= (field0Proc - field0);
      field0 = field0Proc;
    }
  if (nFieldsProc > 0 && nFieldsProc < nFields) nFields = nFieldsProc;

#ifdef SKIP
  // If outputDir contains the string %d, replace it with the camera column
  char outputDirString[160];
  char *newOutputDir = (char *) outputDir;
  if (outputDir != NULL)
    if (strstr(outputDir, "%d") != NULL)
      {
	sprintf(outputDirString, outputDir, camCol);
	newOutputDir = outputDirString;
      }
#endif

  // Get current best astrometric calibration
  ooHandle(arAstromCalib) astrom =
    run->currentAstrometricCalibration(oocRead);
  if (! astrom.isValid())
    {
      shErrStackPush("tsSegmentFetch: no current astrometric calibration");
      return SH_GENERIC_ERROR;
    }
  if (fabs(astrom->equinox() - 2000.) > 0.000001)
    {
      shErrStackPush("tsSegmentCalibrate: equinox not 2000");
      return SH_GENERIC_ERROR;
    }
  int astromField0 = astrom->field0();
  if (astromField0 > field0 || astromField0+astrom->nFields() < field0+nFields)
    {
      shErrStackPush("tsSegmentCalibrate: astrometric calibration doesn't cover field range of the frames pipeline run");
      return SH_GENERIC_ERROR;
    }
  double node = astrom->node();
  double incl = astrom->incl();

  // Set up array of camRow versus filter, in same order as filter order in
  // array fields in TA_CALIB_OBJ (ugriz).  This is used when getting TRANS
  // structures for specific filters.
  const int camRow[5] = {3, 5, 1, 2, 4};

  // All Frames Pipeline Runs were required when stuffed to have r' as their
  // reference filter for the primary positions.
  char refFilter = 'r';
  int refIndex = 2;

  // Get final photometric calibration
  ooHandle(arFinalPhotoCalib) photom =
    run->currentPhotometricCalibration(oocRead);
  if (! photom.isValid())
    {
      shErrStackPush("tsSegmentCalibrate: no current photometric calibration");
      return SH_GENERIC_ERROR;
    }
  int photomField0 = photom->field0();
  if (photomField0 > field0 || photomField0+photom->nFields() < field0+nFields)
    {
      shErrStackPush("tsSegmentCalibrate: photometric calibration doesn't cover field range of the frames pipeline run");
      return TCL_ERROR;
    }

  // Get postage stamp calibration
  ooHandle(arPSCalib) ps = run->postageStampCalibration(oocRead);

  // Get the drift scan
  ooHandle(arRun) scan = run->imagingRun();

  // Get exposure time (seconds)
  // HARWIRE: all photometric CCDs are 2048 rows long
  int nRowsPhys = 2048;
  double expTime = scan->c_obs() * nRowsPhys * 1.e-6;

#ifdef SKIP
  // Create output PDUs
  char fitsFile[100];
  HDR *hdr;
      // Stuff header
      hdr = shHdrNew();
      if (hdr == NULL)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't allocate header");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertLogical(hdr, "SIMPLE", 1, "") != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword SIMPLE");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertInt(hdr, "BITPIX", 8, "") != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword BITPIX");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertInt(hdr, "NAXIS", 0, "") != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword NAXIS");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertLogical(hdr, "EXTEND", 1, "") != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword EXTEND");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertInt(hdr, "RUN", runId, "Imaging run number.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword RUN");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertInt(hdr, "CAMCOL", camCol, "Column in imaging camera.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword CAMCOL");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertInt(hdr, "RERUN", rerun, "Rerun number.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword RERUN");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertInt(hdr, "FIELD0", field0, "First field reduced.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FIELD0");
	  return SH_GENERIC_ERROR;
	}
      if (split)
	{
	  if (shHdrInsertInt(hdr, "NFIELDS", 1,
			     "Number of fields reduced.")!= SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword NFIELDS");
	      return SH_GENERIC_ERROR;
	    }
	}
      else
	{
	  if (shHdrInsertInt(hdr, "NFIELDS", nFields,
			     "Number of fields reduced.")!= SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword NFIELDS");
	      return SH_GENERIC_ERROR;
	    }
	}
      if (shHdrInsertInt(hdr, "STRIPE", scan->stripe(), "Stripe number.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword STRIPE");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "STRIP",
			   scan->strip() == AR_STRIP_NORTH ? "N" :
			   (scan->strip() == AR_STRIP_SOUTH ? "S" : "O"),
			   "N - north, S - south, O - other.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword STRIP");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertDbl(hdr, "EQUINOX", astrom->equinox(),
			 "Equinox of great circle scanned (Julian years).")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword EQUINOX");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertDbl(hdr, "NODE", astrom->node(),
			 "Ascending node of great circle scanned (deg).")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword NODE");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertDbl(hdr, "INCL", astrom->incl(),
			   "Inclination of great circle scanned (deg).")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword INCL");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertDbl(hdr, "C_OBS", scan->c_obs(),
			 "CCD clock rate (usec/unbinned-row).")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword C_OBS");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertDbl(hdr, "TRACKING", scan->tracking(),
			 "Tracking rate (arcsec/TAI-sec).")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword TRACKING");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "PHOTO_ID", run->id(),
			   "Photometric pipeline run id.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword PHOTO_ID");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "PHOT_VER", run->version(),
			   "Version of frames pipeline used.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword PHO_VER");
	  return SH_GENERIC_ERROR;
	}
      ooHandle(arAstromCalib) targetAstrom =
	segment->targetAstrometricCalibration(oocRead);
      if (shHdrInsertAscii(hdr, "TASTR_ID", targetAstrom->id(),
			   "Astrometric pipeline run id for target.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword ASTRO_ID");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "TAST_VER", targetAstrom->version(),
			   "Version of astrometric pipeline for target.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword ASTR_VER");
	  return SH_GENERIC_ERROR;
	}
      ooHandle(arFinalPhotoCalib) targetPhotom =
	segment->targetPhotometricCalibration(oocRead);
      if (shHdrInsertAscii(hdr, "TFCAL_ID", targetPhotom->id(),
			   "Photometric calibration run id for target.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FCALI_ID");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "TFCA_VER", targetPhotom->version(),
			   "Version of fcalib for target.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FCAL_VER");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "EASTR_ID", astrom->id(),
			   "Astrometric pipeline run id for export.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword ASTRO_ID");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "EAST_VER", astrom->version(),
			   "Version of astrometric pipeline for export.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword ASTR_VER");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "EFCAL_ID", photom->id(),
			   "Photometric calibration run id for export.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FCALI_ID");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "EFCA_VER", photom->version(),
			   "Version of fcalib for export.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FCAL_VER");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertAscii(hdr, "FILTERS", "u g r i z",
			   "Filter order in arrays.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FILTERS");
	  return SH_GENERIC_ERROR;
	}

      // Finish object header
      if (shHdrInsertInt(hdr, "REF_BAND", 2, "Band used for canonical centre.")
	  != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword REF_BAND");
	  return SH_GENERIC_ERROR;
	}
      int lineNo;
      if (shHdrGetLineno(hdr, "REF_BAND", &lineNo) != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't find header keyword GAIN");
	  return SH_GENERIC_ERROR;
	}
      if (shHdrInsertLine(hdr, ++lineNo, "END") != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword END");
	  return SH_GENERIC_ERROR;
	}

      // Write out PDUs
      if (!split)
	{
	  if (newOutputDir == NULL)
	    sprintf(fitsFile, "tsObj-%06d-%d-%d-%04d.fit", runId,
		    camCol, rerun, field0);
	  else
	    sprintf(fitsFile, "%s/tsObj-%06d-%d-%d-%04d.fit", newOutputDir,
		    runId, camCol, rerun, field0);
	  if (shHdrWriteAsFits(hdr, fitsFile, STANDARD, 0, DEF_NONE,
			       NULL, NULL, NULL) != SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentCalibrate: couldn't write object file PDU");
	      return SH_GENERIC_ERROR;
	    }
	}
      shHdrDel(hdr);
#endif
  
  // Initialize arrays of TRANS structures.
  TRANS astroTrans[nFilters];

#ifdef SKIP
  // Initialize output chain
  CHAIN *objChain;
#endif

  // Mode to open objects in
  ooMode mode = oocRead;;

  /*adjacent filter to use in colorterms*/
  const int adj[nFilters]  = {1, 2, 3, 4, 3};
  /* sign to use for color terms */
  const int sign[nFilters] = {1, 1, 1, 1, -1}; 

  // For each field ...
  int lastField = field0 + nFields - 1;
  printf("run %d, camCol %d, rerun %d, fields %d - %d\n", runId, camCol,
	 rerun, field0, lastField);
  for (int field = field0; field <= lastField; field++)
    {
      printf("%4d", field);
      fflush(stdout);

      // Get object list for this field
      ooHandle(arObjectList) objList = run->objectList(field);

      // If this is a MISSING or HOLE field skip it.
      if (objList->quality(skyVersion) == AR_DQ_MISSING ||
	  objList->quality(skyVersion) == AR_DQ_HOLE)
	continue;

      // Get the astrometric and photometric transformations for this field
      arFieldPTrans fieldPTrans = photom->ptrans(field);
      int k;
      for (k = 0; k < nFilters; k++)
	{
	  int ccd = camRow[k] * 10 + camCol;
	  astroTrans[k] = astrom->astrotoolsTrans(ccd, field);
	}
      arFieldPhotoTrans psTrans = ps->trans(field);

      // Convert photometric transformations to form callable by
      // "target" routines.
      TA_PTRANS ptrans;
      ptrans.field = field;
      for (k = 0; k < nFilters; k++)
	{
	  arPTrans dbPTrans = fieldPTrans.filter((arFilter)k);
	  ptrans.a[k].val = dbPTrans.a().value();
	  ptrans.a[k].err = dbPTrans.a().err();
	  ptrans.b[k].val = dbPTrans.b().value();
	  ptrans.b[k].err = dbPTrans.b().err();
	  ptrans.c[k].val = dbPTrans.c().value();
	  ptrans.c[k].err = dbPTrans.c().err();
	  ptrans.k[k].val = dbPTrans.k().value();
	  ptrans.k[k].err = dbPTrans.k().err();
	  double calMag, calMagErr;
	  int caliStatus;
	  taLCaliApply(k,0,psTrans.filter((arFilter)k).flux20().value(),0.,
		       -1000., -1000., expTime, expTime,
		       astroTrans[k].airmass, ptrans.a[k].val,
		       ptrans.k[k].val, 0., 0., sign[k], 0., 0.,
		       &calMag, &calMagErr, &caliStatus);
	  ptrans.dZero[k] = calMag - 20.;
	}

#ifdef SKIP
      // Create chain for calibrated objects for this field
      objChain = shChainNew("TS_PM_OBJ");
#endif
	  
      // Initialize iterator over objects
      ooItr(arObject) objItr;
      if (objList->objects(objItr, mode) != oocSuccess)
	{
	  shErrStackPush("tsSegmentFetch: couldn't initialize iterator over objects");
	  return SH_GENERIC_ERROR;
	}

      // For each object ...
      while (objItr.next())
	{
	  // Output only PRIMARY, or if requested, OK_SCANLINE objects
	  if (! (objItr->status(skyVersion) & AR_OBJECT_STATUS_PRIMARY) &&
	      ! (scanline &&
		 (objItr->status(skyVersion) & AR_OBJECT_STATUS_OK_SCANLINE)))
	    continue;
	  if (field == segField0)
	    {
	      if (!(objItr->status(skyVersion) & AR_OBJECT_STATUS_FIRST_FIELD))
		continue;
	    }
	  else if (objItr->status(skyVersion) & AR_OBJECT_STATUS_FIRST_FIELD)
	    {
	      continue;
	    }
	  
	  // Calibrate the object, transforming instrument quantities to
	  // calibrated quantities in the TA_CALIB_OBJ C structure.
	  TS_PM_OBJ *objPtr =
	    (TS_PM_OBJ*) shMalloc(sizeof(TS_PM_OBJ));
	  objPtr->run = runId;
	  objPtr->camCol = camCol;
	  objPtr->rerun = rerun;
	  objPtr->field = field;
	  objPtr->id = objItr->id();
	  objPtr->objc_type = objItr->objc_type();
	  objPtr->objc_flags = objItr->objc_flags();
	  objPtr->objc_flags2 = objItr->objc_flags2();
	  objPtr->rowv = objItr->rowv().value();
	  objPtr->rowvErr = objItr->rowv().err();
	  objPtr->colv = objItr->colv().value();
	  objPtr->colvErr = objItr->colv().err();
	  objPtr->status = objItr->status(skyVersion);
	  objPtr->primTarget = objItr->primaryType(skyVersion);
	  if (objItr->properMotionMatch(skyVersion))
	    {
	      objPtr->properMotionMatch = 1;
              objPtr->properMotionDelta = objItr->pm(skyVersion);
              objPtr->properMotion = objItr->properMotion(skyVersion);
	      objPtr->properMotionAngle = objItr->pmAngle(skyVersion);
	      objPtr->usnoBlue = objItr->blue(skyVersion);
	      objPtr->usnoRed = objItr->red(skyVersion);
	    }
	  else
	    {
	      objPtr->properMotionMatch = 0;
	      objPtr->properMotionDelta = 0.0;
	      objPtr->properMotion = 0.0;
	      objPtr->properMotionAngle = 0.0;
	      objPtr->usnoBlue = 0.0;
	      objPtr->usnoRed = 0.0;
	    }
          objPtr->mjd = astroTrans[refIndex].mjd;
	  for (int m = 0; m < 5; m++)
	    {
	      int c = cosmicIdx[m];
	      double mag, magErr;
	      int status;
	      if (objPtr->objc_type == AR_OBJECT_TYPE_STAR)
		taLCaliApply(m,0,objItr->psfCounts((arFilter)m).value(),
			     objItr->psfCounts((arFilter)m).err(),
			     objItr->psfCounts((arFilter)adj[m]).value(),
			     objItr->psfCounts((arFilter)adj[m]).err(),
			     expTime, expTime, astroTrans[m].airmass,
			     ptrans.a[m].val, ptrans.k[m].val,
			     ptrans.b[m].val,ptrans.c[m].val,
			     sign[m], cosmicColor[c], cosmicScatter[c],
			     &mag, &magErr, &status);
	      else
		taLCaliApply(m,0,objItr->counts_model((arFilter)m).value(),
			     objItr->counts_model((arFilter)m).err(),
			     objItr->counts_model((arFilter)adj[m]).value(),
			     objItr->counts_model((arFilter)adj[m]).err(),
			     expTime, expTime, astroTrans[m].airmass,
			     ptrans.a[m].val, ptrans.k[m].val,
			     ptrans.b[m].val,ptrans.c[m].val,
			     sign[m], cosmicColor[c], cosmicScatter[c],
			     &mag, &magErr, &status);
	      objPtr->mag[m] = mag;
	      objPtr->magErr[m] = magErr;
	    }
	  atTransApply(&astroTrans[refIndex], refFilter,
		       objItr->objc_rowc().value(),
		       0., objItr->objc_colc().value(), 0.,
		       objPtr->mag, objPtr->magErr, &objPtr->mu, NULL,
		       &objPtr->nu, NULL);
	  atGCToEq(objPtr->mu, objPtr->nu, &(objPtr->ra), &(objPtr->dec),
		   node, incl);
	  atEqToGal(objPtr->ra, objPtr->dec, &(objPtr->l), &(objPtr->b));
          if (atExtinction(objPtr->l, objPtr->b, &(objPtr->reddening[0]),
	      &(objPtr->reddening[1]), &(objPtr->reddening[2]),
	      &(objPtr->reddening[3]), &(objPtr->reddening[4])) != SH_SUCCESS)
             return SH_GENERIC_ERROR;

          // Skip objects fainter than 21.5 in r' and 21.5 in g'
          if (objPtr->mag[2] > 21.5 && objPtr->mag[1] > 21.5)
            {
              shFree(objPtr);
              continue;
            }
#ifdef SKIP
          // Skip objects fainter than 21 in all bands
          if (objPtr->mag[0] > 21 &&
	      objPtr->mag[1] > 21 &&
	      objPtr->mag[2] > 21 &&
	      objPtr->mag[3] > 21 &&
	      objPtr->mag[4] > 21)
            {
              shFree(objPtr);
              continue;
            }
#endif

	  // Append to exported chain of calibrated objects.
	  if (shChainElementAddByPos(objChain, objPtr, "TS_PM_OBJ",
				     TAIL, AFTER)
	      != SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentFetch: couldn't append object chain");
	      return SH_GENERIC_ERROR;
	    }
	} // end loop over objects

#ifdef SKIP
      // Write the objects file
      // Set file name
      if (newOutputDir == NULL)
	sprintf(fitsFile, "tsObj-%06d-%d-%d-%04d.fit", runId,
		camCol, rerun, (split ? field : field0));
      else
	sprintf(fitsFile, "%s/tsObj-%06d-%d-%d-%04d.fit", newOutputDir,
		runId, camCol, rerun, (split ? field : field0));

      if (split)
	{
	  // One file per field.  Write the PDU.
	  if (shHdrDelByKeyword(hdr, "FIELD0") != SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentCalibrate: couldn't delete header keyword");
	      return SH_GENERIC_ERROR;
	    }
	  if (shHdrInsertInt(hdr, "FIELD0", field, "First field reduced.")
	      != SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentCalibrate: couldn't insert header keyword FIELD0");
	      return SH_GENERIC_ERROR;
	    }
	  if (shHdrWriteAsFits(hdr, fitsFile, STANDARD, 0, DEF_NONE,
			       NULL, NULL, NULL) != SH_SUCCESS)
	    {
	      shErrStackPush("tsSegmentCalibrate: couldn't write object file PDU");
	      return SH_GENERIC_ERROR;
	    }
	}

      // Write the objects HDU
      TBLCOL *tblCol;
      if ((tblCol = taCalibObjsWriteToTbl(objChain)) == NULL)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't append object HDU");
	  return SH_GENERIC_ERROR;
	}
      taGenericChainDestroy(objChain);
      shHdrInsertInt(&(tblCol->hdr), "FIELD", field,
		     "Field sequence number within the run");
      shHdrInsertInt(&(tblCol->hdr), "QUALITY", objList->quality(skyVersion),
		       "Field quality");
      shHdrInsertInt(&(tblCol->hdr), "CULLED",
		     objList->culled(skyVersion),
		     "Culling bit mask");
      HANDLE handle;
      handle.ptr = (void *) tblCol;
      handle.type = shTypeGetFromName("TBLCOL");
      if (shFitsWrite(&handle, (char *)fitsFile, DEF_NONE, f_hduBINTABLE,
		      SH_FITS_PDU_NONE, 1, STANDARD) != SH_SUCCESS)
	{
	  shErrStackPush("tsSegmentCalibrate: couldn't append object HDU");
	  return SH_GENERIC_ERROR;
	}
      shTblcolDel(tblCol);
#endif
    } // end loop over fields
  printf("\n");

  // Clean up
#ifdef SKIP
if (split) shHdrDel(hdr);
#endif
  
  return SH_SUCCESS;
}