#include <TObject.h>
#include <TMath.h>

#include <RAT/zdab_convert.hh>
#include <RAT/PZdabFile.hh>
#include <RAT/DB.hh>
#include <RAT/DS/Entry.hh>
#include <RAT/DS/Run.hh>
#include <RAT/DS/TrigHeader.hh>
#include <RAT/DS/ECAHeader.hh>
#include <RAT/DS/Digitiser.hh>
#include <RAT/DS/PMT.hh>
#include <RAT/DS/UniversalTime.hh>
#include <RAT/DU/PMTInfo.hh>
#include <RAT/MetaInformation.hh>
#include <RAT/Log.hh>
using RAT::warn;
using RAT::info;
using RAT::detail;
using RAT::debug;
using RAT::Log;

#include <utility>
#include <string>
#include <vector>
#include <map>
#include <algorithm>

#include <CLHEP/Units/SystemOfUnits.h>
using namespace CLHEP;
#include <stdint.h>

namespace ratzdab {

/* RATDB instance to get at PMTINFO, since PMT type is needed to
   populate the RAT DS. */
class RATDB {
    public:
        RATDB() {
            RAT::DB* db = RAT::DB::Get();
            assert(db);

            pmtinfo = db->GetLink("PMTINFO");
            assert(pmtinfo);
            daq_rl = db->GetLink("DAQ_RUN_LEVEL");
            assert(daq_rl);
            // cache these for performance
            pmttype = pmtinfo->GetIArray("pmt_type");
            caen_channels = daq_rl->GetIArray("trigSums");
        }

        virtual ~RATDB() {}

        RAT::DBLinkPtr pmtinfo;
        RAT::DBLinkPtr daq_rl;
        std::vector<int> pmttype;
        std::vector<int> caen_channels;
};

/* initialize static maps */

std::map<int, int> get_pdg_to_snoman_map() {
    std::map<int, int> m;
    m[0] = 1;  // Photon, G4 also uses 0 for geantino and unknown
    m[22] = 2;  // Gamma
    m[11] = 20;  // Electron
    m[-11] = 21;  // Positron
    m[13] = 22;  // Mu-
    m[-13] = 23;  // Mu+
    m[15] = 24;  // Tau-
    m[-15] = 25;  // Tau+
    m[12] = 30;  // Nu-e
    m[-12] = 31;  // Nu-e-bar
    m[14] = 32;  // Nu-mu
    m[-14] = 33;  // Nu-mu-bar
    m[16] = 34;  // Nu-tau
    m[-16] = 35;  // Nu-tau-bar
    m[111] = 40;  // Pi0
    m[211] = 41;  // Pi+
    m[-211] = 42;  // Pi-
    m[311] = 50;  // K0
    m[-311] = 51;  // K0-bar
    m[321] = 52;  // K+
    m[-321] = 53;  // K-
    m[2212] = 80;  // Proton
    m[2112] = 81;  // Neutron
    return m;
}
static std::map<int, int> pdg_to_snoman = get_pdg_to_snoman_map();

std::map<int, float> get_pdg_to_mass_map() {
    std::map<int, float> m;
    m[0] = 0;  // Photon, G4 also uses 0 for geantino and unknown
    m[22] = 0;  // Gamma
    m[11] = 0.510998928;  // Electron
    m[-11] = 0.510998928;  // Positron
    m[13] = 105.6583715;  // Mu-
    m[-13] = 105.6583715;  // Mu+
    m[15] = 1776.82;  // Tau-
    m[-15] = 1776.82;  // Tau+
    m[12] = 0;  // Nu-e
    m[-12] = 0;  // Nu-e-bar
    m[14] = 0;  // Nu-mu
    m[-14] = 0;  // Nu-mu-bar
    m[16] = 0;  // Nu-tau
    m[-16] = 0;  // Nu-tau-bar
    m[111] = 134.9766;  // Pi0
    m[211] = 139.57018;  // Pi+
    m[-211] = 139.57018;  // Pi-
    m[311] = 497.648;  // K0
    m[-311] = 497.648;  // K0-bar
    m[321] = 493.667;  // K+
    m[-321] = 493.667;  // K-
    m[2212] = 938.272046;  // Proton
    m[2112] = 939.565378;  // Neutron
    return m;
}
static std::map<int, float> pdg_to_mass = get_pdg_to_mass_map();

// Prints a message about the unexpected subrecord named 'type' if we
// haven't warned about that one before. Meant for when we find an
// unexpected, yet valid, subrecord type. We don't want to abort, and we
// don't want to flood the user with messages, but we do want to remark
// upon it. Not written to be efficient, since it isn't supposed to be
// called under normal circumstances.
static void unexpected_subrecord_warn(const char * const type)
{
  static std::vector<std::string> said;
  if(std::find(said.begin(), said.end(), type) == said.end()){
    warn << "Warning: Got " << type << " subrecord, which is never "
            "expected\nin SNO+ ZDAB data files.  Corrupt input?\n";
    said.push_back(type);
  }
}

////////////////////////////////////////////////////////////////////////
//                                                                    //
//                        Unpacking functions                         //
//                                                                    //
////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////
// unpack::event
RAT::DS::Entry* unpack::event(PmtEventRecord* const pev) {
    RAT::DS::Entry* ds = new RAT::DS::Entry;
    RAT::DS::EV ev;
    static RATDB ratdb;

    const unsigned nhit = pev->NPmtHit;
    const uint32_t run_id = pev->RunNumber;
    const uint16_t subrun_id = pev->DaqStatus; // seriously
    //uint32_t evorder = pev->EvNumber;  // unused in RAT
    //uint16_t datatype = pev->DataType;  // unused here...?
    uint32_t mtc_words[6];
    memcpy(mtc_words, &pev->TriggerCardData, 6*sizeof(uint32_t));
    const uint32_t gtid = static_cast<uint32_t>(UNPK_MTC_GT_ID(mtc_words));
    const uint64_t clockcount10 =
      (static_cast<uint64_t>(UNPK_MTC_BC10_2(mtc_words)) << 32) |
      (static_cast<uint64_t>(UNPK_MTC_BC10_1(mtc_words)));
    const uint64_t clockcount50 =
      (static_cast<uint64_t>(UNPK_MTC_BC50_2(mtc_words)) << 11) |
      (static_cast<uint64_t>(UNPK_MTC_BC50_1(mtc_words)));
    const uint32_t trig_error = UNPK_MTC_ERROR(mtc_words);
    const uint32_t esumpeak = UNPK_MTC_PEAK(mtc_words);
    const uint32_t esumdiff = UNPK_MTC_DIFF(mtc_words);
    const uint32_t esumint = UNPK_MTC_INT(mtc_words);

    // trigger word
    const uint32_t trig_type = ((mtc_words[3] & 0xff000000) >> 24) |
                               ((mtc_words[4] & 0x7ffff) << 8);

    ds->SetRunID(run_id);
    ds->SetSubRunID(subrun_id);

    ev.SetClockStat10(0);  // FIXME: what is this?
    ev.SetTrigError(trig_error);
    ev.SetESumPeak(esumpeak);
    ev.SetESumDerivative(esumdiff);
    ev.SetESumIntegral(esumint);
    ev.SetTrigType(trig_type);
    ev.SetGTID(gtid);
    ev.SetClockCount50(clockcount50);
    ev.SetClockCount10(clockcount10);

    // set ut from 10mhz clock
    const uint64_t total = clockcount10 * 100;
    const uint64_t ns = total % (static_cast<uint64_t>(1e9));
    uint64_t s = total / 1e9;
    const uint64_t d = s / 86400;
    s -= (86400 * d);

    ev.SetUniversalTime(RAT::DS::UniversalTime(d, s, ns));

    // loop over sub fields and extract extra info
    PMTCal* calhits = static_cast<PMTCal*>(NULL);
    MonteCarloHeader* mcdata = static_cast<MonteCarloHeader*>(NULL);
    FittedEvent* fitdata = static_cast<FittedEvent*>(NULL);
    uint32_t* caen_data = static_cast<uint32_t*>(NULL);
    uint32_t* tubii_data = static_cast<uint32_t*>(NULL);
    uint32_t* extrahitinfo = static_cast<uint32_t*>(NULL);
    uint32_t* extraeventinfo = static_cast<uint32_t*>(NULL);
    unsigned nfits = 0;

    uint32_t* sub_header = &pev->CalPckType;
    while (*sub_header & SUB_NOT_LAST) {
        sub_header += (*sub_header & SUB_LENGTH_MASK);
        const uint32_t subtype_id = *sub_header >> SUB_TYPE_BITNUM;

        switch(subtype_id){
            case SUB_TYPE_CALIBRATED:
            unexpected_subrecord_warn("SUB_TYPE_CALIBRATED");
            calhits = reinterpret_cast<PMTCal*>(sub_header + 1);
            break;

            case SUB_TYPE_MONTE_CARLO:
            unexpected_subrecord_warn("SUB_TYPE_MONTE_CARLO");
            mcdata = reinterpret_cast<MonteCarloHeader*>(sub_header + 1);
            break;

            case SUB_TYPE_FIT:
            unexpected_subrecord_warn("SUB_TYPE_FIT");
            fitdata = reinterpret_cast<FittedEvent*>(sub_header + 1);
            nfits = ((*sub_header & SUB_LENGTH_MASK) *
                    sizeof(uint32_t) - sizeof(SubFieldHeader)) /
                    sizeof(FittedEvent);
            break;

            case SUB_TYPE_CAEN:
            // CAEN sub-records *are* expected in SNO+ data
            caen_data = static_cast<uint32_t*>(sub_header + 1);
            break;

            case SUB_TYPE_HIT_DATA:
            unexpected_subrecord_warn("SUB_TYPE_HIT_DATA");
            extrahitinfo = static_cast<uint32_t*>(sub_header + 1);
            break;

            case SUB_TYPE_EVENT_DATA:
            unexpected_subrecord_warn("SUB_TYPE_EVENT_DATA");
            extraeventinfo = static_cast<uint32_t*>(sub_header + 1);
            break;

            case SUB_TYPE_TUBII:
            // TUBii sub-record is also expected in SNO+ data
            tubii_data = static_cast<uint32_t*>(sub_header + 1);
            break;

            default:
            warn << dformat("ERROR: Got subrecord with unknown type "
              "0x%x. Corrupt input?\n", subtype_id);
            break;
        }
    }

    // pmt hit data
    uint32_t* hits = reinterpret_cast<uint32_t*>(pev + 1);
    RAT::DS::UncalPMTs& uncalPMTs = ev.GetUncalPMTs();
    double totalQ = 0.0;
    RAT::DS::CalPMTs calPMTs;
    for (unsigned i=0; i<nhit; i++, hits += 3) {
        const unsigned crate_id = UNPK_CRATE_ID(hits);
        const unsigned board_id = UNPK_BOARD_ID(hits);
        const unsigned channel_id = UNPK_CHANNEL_ID(hits);

        const int lcn = 16 * 32 * crate_id + 32 * board_id + channel_id;
        const int type = ratdb.pmttype[lcn];

        RAT::DS::PMTUncal pmtuncal = unpack::pmt(hits, lcn);
        uncalPMTs.AddPMT( pmtuncal, static_cast<RAT::DU::PMTInfo::EPMTType>( type ) );
        totalQ += pmtuncal.GetQHS();

        if (calhits)
          {
            RAT::DS::PMTCal pmtcal;
            pmtcal.SetID(pmtuncal.GetID());
            pmtcal.SetChanFlags(pmtuncal.GetChanFlags());
            pmtcal.SetCellID(pmtuncal.GetCellID());

            pmtcal.SetQHS(calhits->qhs);
            pmtcal.SetQHL(calhits->qhl);
            pmtcal.SetQLX(calhits->qlx);
            pmtcal.SetTime(calhits->tac);
            calPMTs.AddPMT( pmtcal, static_cast<RAT::DU::PMTInfo::EPMTType>( type ) );
          }
    }
    ev.SetNhits( uncalPMTs.GetCount() );
    ev.SetTotalCharge(totalQ);
    if( calPMTs.GetAllCount() > 0 )
      ev.AddCalPMTs( calPMTs, -1 ); // FIXME choose an appropriate calibration type for this, -1 is unknown?

    // caen data
    if (caen_data) {
      RAT::DS::Digitiser dig = unpack::caen(caen_data);
      ev.SetDigitiser(dig);
    }

    // tubii data
    if (tubii_data) {
      uint32_t tubiiWord = unpack::tubii(tubii_data);
      ev.SetTubiiTrig(tubiiWord);
    }

    // extra hit info -- no place in RAT DS currently
    if (extrahitinfo) {
        warn << "ratzdab::unpack::event: Extra hit data not converted "
          "since\nthere is nowhere in RAT::DS to put it." << newline;
    }

    // extra event data -- no place in RAT DS currently
    if (extraeventinfo) {
        ExtraEventData* eed = reinterpret_cast<ExtraEventData*>(extraeventinfo);
        const float value = eed->value;
        warn << "ratzdab::unpack::event: Extra event data with value "
             << value << " not\nconverted since there is nowhere in "
             "RAT::DS to put it." << newline;
    }

    // monte carlo data
    if (mcdata) {
        MonteCarloHeader* mchdr = mcdata;
        MonteCarloVertex* mcvtx = reinterpret_cast<MonteCarloVertex*>(mchdr+1);
        for (unsigned i=0; i<mchdr->nVertices; i++, ++mcvtx) {
            RAT::DS::MCParticle p;
            const int snoman_code = mcvtx->particle;

            // convert cm -> mm
            p.SetPosition(TVector3(mcvtx->x * 10, mcvtx->y * 10, mcvtx->z * 10));

            // convert SNOMAN code to PDG
            int pdgcode = 0;
            bool found_code = false;
            std::map<int, int>::const_iterator it_code;
            for (it_code  = ratzdab::pdg_to_snoman.begin();
                 it_code != ratzdab::pdg_to_snoman.end();
                 it_code++) {
                if (it_code->second == snoman_code) {
                    pdgcode = it_code->first;
                    found_code = true;
                    break;
                }
            }
            if (!found_code) {
                warn << "ratzdab::unpack::event: No PDG code available for "
                  "SNOMAN code " << snoman_code << ", using zero." << newline;
            }
            p.SetPDGCode(pdgcode);

            // convert total energy to kinetic
            float mass = 0;
            std::map<int, float>::const_iterator it_mass =
              ratzdab::pdg_to_mass.find(pdgcode);
            if (it_mass != ratzdab::pdg_to_mass.end()) {
                mass = it_mass->second;
            }
            else {
                warn << "ratzdab::unpack::event: No mass code available for "
                  "PDG code " << pdgcode << ", using zero." << newline;
            }
            float momentum = mcvtx->energy * mcvtx->energy - mass * mass;
            float ke = momentum * momentum / (2 * mass);
            p.SetKineticEnergy(ke);

            // convert direction cosines to cartesian momentum vector
            double px = momentum * TMath::Cos(mcvtx->u);
            double py = momentum * TMath::Cos(mcvtx->v);
            double pz = momentum * TMath::Cos(mcvtx->w);
            p.SetMomentum(TVector3(px, py, pz));

            p.SetTime(mcvtx->time);
            ds->GetMC().AddMCParticle(p);
        }
    }

    // fit results
    if (fitdata) {
      if(nfits > 1000){
        warn << dformat("SUB_TYPE_FIT subrecord claims it has %u "
          "fits.\nDoubtful. Skipping this subrecord.\n", nfits);
      }
      else{
        const size_t pass = RAT::MetaInformation::Get()->GetCurrentPass();
        for (unsigned i=0; i<nfits; i++, ++fitdata) {
          RAT::DS::FitVertex fv;
          fv.SetPosition(TVector3(fitdata->x, fitdata->y, fitdata->z));
          fv.SetDirection(TVector3(fitdata->u, fitdata->v, fitdata->w));
          fv.SetTime(fitdata->time);

          RAT::DS::FitResult fr;
          fr.SetFOM("quality", fitdata->quality);
          fr.SetVertex(0, fv);

          ev.SetFitResult(pass, fitdata->name, fr);
        }
      }
    }

    ds->AddEV( ev );
    return ds;
}

////////////////////////////////////////////////////////////////////////
// unpack::hrdr
RAT::DS::Run* unpack::rhdr(RunRecord* const r) {
    RAT::DS::Run* run = new RAT::DS::Run;

    run->SetDate(r->Date);
    run->SetTime(r->Time);
    run->SetRunID(r->RunNumber);
    run->SetSubRunID(r->SubRunNumber);
    run->SetRunType(r->RunMask);
    run->SetFirstGTID(r->FirstGTID); // First GTID
    run->SetValidGTID(r->ValidGTID); // First valid GTID (HW settings known)
    uint32_t version = r->Spares[0]; // version stored in format 0x04010300 = 4.1.3
    int v1 = (version >> 24);
    int v2 = (version >> 16)%16;
    int v3 = (version >> 8)%16;
    char builderversion[32];
    sprintf(builderversion, "%i.%i.%i",v1,v2,v3);
    run->SetBuilderVersion((std::string)builderversion);
    return run;

}

////////////////////////////////////////////////////////////////////////
// unpack::trig
RAT::DS::TrigHeader* unpack::trig(TriggerInfo* const t) {
    RAT::DS::TrigHeader* triginfo = new RAT::DS::TrigHeader;

    triginfo->SetTrigMask(t->TriggerMask);
    triginfo->SetPulserRate(t->PulserRate);
    triginfo->SetMTC_CSR(t->ControlRegister);
    triginfo->SetLockoutWidth(t->reg_LockoutWidth);
    triginfo->SetPrescaleFrequency(t->reg_Prescale);

    // triggers are separate members in TriggerInfo, array in TrigHeader
    triginfo->SetTriggerThreshold(0, t->n100lo);
    triginfo->SetTriggerThreshold(1, t->n100med);
    triginfo->SetTriggerThreshold(2, t->n100hi);
    triginfo->SetTriggerThreshold(3, t->n20);
    triginfo->SetTriggerThreshold(4, t->n20lb);
    triginfo->SetTriggerThreshold(5, t->esumlo);
    triginfo->SetTriggerThreshold(6, t->esumhi);
    triginfo->SetTriggerThreshold(7, t->owln);
    triginfo->SetTriggerThreshold(8, t->owlelo);
    triginfo->SetTriggerThreshold(9, t->owlehi);
    triginfo->SetTrigZeroOffset(0, t->n100lo_zero);
    triginfo->SetTrigZeroOffset(1, t->n100med_zero);
    triginfo->SetTrigZeroOffset(2, t->n100hi_zero);
    triginfo->SetTrigZeroOffset(3, t->n20_zero);
    triginfo->SetTrigZeroOffset(4, t->n20lb_zero);
    triginfo->SetTrigZeroOffset(5, t->esumlo_zero);
    triginfo->SetTrigZeroOffset(6, t->esumhi_zero);
    triginfo->SetTrigZeroOffset(7, t->owln_zero);
    triginfo->SetTrigZeroOffset(8, t->owlelo_zero);
    triginfo->SetTrigZeroOffset(9, t->owlehi_zero);

    return triginfo;
}

////////////////////////////////////////////////////////////////////////
// unpack::eped
RAT::DS::ECAHeader* unpack::eped(EpedRecord* const e) {

    if (e->Flag == 0x01000000) { //orca subrun record
        TObject* feped = new TObject();
        return reinterpret_cast<RAT::DS::ECAHeader*>(feped);
    }

    uint32_t used_pattern;
    uint32_t used_caltype;
    //If it is SNO+ data
    used_pattern = e->halfCrateID;
    used_caltype = e->CalibrationType;
    //else
    if(e->Flag > 0 && e->Flag < 5 ){ //means that it is SNO data
      //Convert SNO pattern in SNO+ pattern
      if(used_pattern & 128){ //check bit 7
        used_pattern &= ~(1 << 7); //unset bit 7
        used_pattern *=2;
        used_pattern++;
      }
      else{
        used_pattern *=2;
      }
      //Convert SNO calibration type in SNO+ type
      if(used_caltype == 3) used_caltype = 16; //pedestal half-crates calibration
      else if(used_caltype == 2) used_caltype = 26; //tslope half-crates calibration

    }

    RAT::DS::ECAHeader* eped = new RAT::DS::ECAHeader;
    eped->SetGTDelayCoarse(e->ped_delay_coarse);
    eped->SetGTDelayFine(e->ped_delay_fine);
    eped->SetChargePedestalAmplitude(e->qinj_dacsetting);
    eped->SetChargePedestalWidth(e->ped_width);
    eped->SetCalType(used_caltype);
    eped->SetPatternID(used_pattern);
    eped->SetFlag(e->Flag);

    return eped;
}

////////////////////////////////////////////////////////////////////////
// unpack::caen
RAT::DS::Digitiser unpack::caen(uint32_t* const c) {
    RAT::DS::Digitiser d;

    static RATDB ratdb;
    uint32_t channel_mask = UNPK_CAEN_CHANNEL_MASK(c);
    uint32_t magic = UNPK_CAEN_MAGIC(c);
    uint32_t pattern = UNPK_CAEN_PATTERN(c);
    uint32_t event_count = UNPK_CAEN_EVENT_COUNT(c);
    uint32_t clock = UNPK_CAEN_TRIGGER_TIME(c);
    uint32_t pack_flag = UNPK_CAEN_PACK_FLAG(c);

    d.SetBit24(pack_flag);
    d.SetDataFormat(magic);  // FIXME?
    d.SetEventID(event_count);
    d.SetTrigTime(clock);
    d.SetIOPins(pattern);

    unsigned nchannels = 0;
    for (int i=0; i<8; i++) {
        if (channel_mask & (1 << i)) {
            nchannels++;
        }
    }

    if (nchannels > 0) {
        int trace_length = (UNPK_CAEN_WORD_COUNT(c) - 4) / nchannels; // trace_length is the length of the trace in 32bit integers
        int trace_samples = trace_length * 2;  // trace_samples is the number of samples in each trace
        d.SetNWords(trace_length);

        unsigned n = 0;
        for (unsigned i=0; i<8; i++) {
            if (!(channel_mask & (1 << i))) {
                continue;
            }
            std::vector<uint16_t> trace(trace_samples);

            uint32_t* pword = c + 4 + n * trace_length ;
            for (int j=0; j<trace_samples; pword++) {
                //for even j's
                trace[j++] = *(pword) & 0xffff;
                //for odd j's
                trace[j++] = (*(pword) >> 16) & 0xffff;

            }
            n++;

            int id = ratdb.caen_channels[i];
            // If the ratdb table doesn't know what this is just index it by
            // caen channel number
            if(!id){
                id = i;
            }
            d.SetWaveform(id, trace);
        }
    }

    return d;
}

////////////////////////////////////////////////////////////////////////
// unpack::tubii
uint32_t unpack::tubii(uint32_t* const t) {
    uint32_t tubii_word = *t & 0xffffff;  // TUBii's trig word is 24 bit
    return tubii_word;
}

////////////////////////////////////////////////////////////////////////
// unpack::pmt
RAT::DS::PMTUncal unpack::pmt(uint32_t* const hits, const int lcn) {
    const int cell = UNPK_CELL_ID(hits);

    // This can be also be written as
    // (*hits >> 30 & 0x3) | ((*(hits+1)) >> 26 & 0x3c)
    // but the performance gain is unmeasurably small.
    const char flags = (UNPK_CGT_ES_16(hits)    << 0) |
                       (UNPK_CGT_ES_24(hits)    << 1) |
                       (UNPK_MISSED_COUNT(hits) << 2) |
                       (UNPK_NC_CC(hits)        << 3) |
                       (UNPK_LGI_SELECT(hits)   << 4) |
                       (UNPK_CMOS_ES_16(hits)   << 5);

    const uint32_t qhs = UNPK_QHS(hits);
    const uint32_t qhl = UNPK_QHL(hits);
    const uint32_t qlx = UNPK_QLX(hits);
    const uint32_t tac = UNPK_TAC(hits);

    RAT::DS::PMTUncal p;

    p.SetID(lcn);
    p.SetCellID(cell);
    p.SetChanFlags(flags);
    p.SetQHS(qhs);
    p.SetQHL(qhl);
    p.SetQLX(qlx);
    p.SetTime(tac);

    return p;
}

////////////////////////////////////////////////////////////////////////
//                                                                    //
//                         Packing functions                          //
//                                                                    //
////////////////////////////////////////////////////////////////////////

////////////////////////////////////////////////////////////////////////
// pack::event
PmtEventRecord* pack::event(RAT::DS::Entry *ds, int ev_id) {
    if (!ds || ev_id > (int)ds->GetEVCount() - 1) {
        return static_cast<PmtEventRecord*>(NULL);
    }
    RAT::DS::EV& ev = ds->GetEV(ev_id);

    int npmts = ev.GetUncalPMTs().GetAllCount();

    // calculate buffer size (length)
    // there is no GenericRecordHeader here, as they are only used in
    // dispatched data, not in files
    int length = sizeof(PmtEventRecord) + npmts * sizeof(FECReadoutData);

    int nfits = 0;
    int caen_length = 0;

    // calibrated pmt data
    if (ev.GetCalPMTs().GetAllCount() > 0) {
        assert(ev.GetCalPMTs().GetAllCount() == ev.GetUncalPMTs().GetAllCount());
        length += npmts * sizeof(PMTCal) + sizeof(SubFieldHeader);
    }

    // monte carlo data
    if (ds->MCExists()) {
        RAT::DS::MC& mc = ds->GetMC();
        std::vector<unsigned int> mcTrackIDs = mc.GetMCTrackIDs();
        int nvertices = 0;
        if (mc.GetMCTrackCount() > 0) {
            for (size_t itrack=0; itrack < mcTrackIDs.size(); itrack++) {
                RAT::DS::MCTrack& track = mc.GetMCTrack(mcTrackIDs.at(itrack));
                nvertices += track.GetMCTrackStepCount();
            }
        }
        else {
            nvertices = mc.GetMCParticleCount();
        }
        length += nvertices * sizeof(MonteCarloVertex)
                  + sizeof(MonteCarloHeader) + sizeof(SubFieldHeader);
    }

    // fits
    nfits = ev.GetFitNames().size();
    if (nfits > 0) {
        length += nfits * (sizeof(FittedEvent) + sizeof(SubFieldHeader));
    }

    // caen data
    if(ev.DigitiserExists()){
        int caen_total_sample_count = 0;
        RAT::DS::Digitiser & digitizer = ev.GetDigitiser();
        for (int i=0; i < 8; i++)
          if (digitizer.ExistsWaveform(i))
            caen_total_sample_count+=digitizer.GetWaveform(i).size();

        if(caen_total_sample_count){
          const int caen_header_length = 4 * sizeof(uint32_t);
          caen_length = caen_header_length
                      + caen_total_sample_count * sizeof(uint16_t);
          length += sizeof(SubFieldHeader) + caen_length;
        }
    }

    // allocate PmtEventRecord atop char buffer
    char* buffer = new char[length];
    PmtEventRecord* r = reinterpret_cast<PmtEventRecord*>(buffer);

    // get pointer to first MTC and FEC words
    uint32_t* mtc_word = reinterpret_cast<uint32_t*>(&(r->TriggerCardData));
    uint32_t* fec_word = reinterpret_cast<uint32_t*>(r + 1);

    r->RunNumber = ds->GetRunID();
    //r->EvNumber = ev.GetGTID();  // event order not implemented in RAT
    r->NPmtHit = npmts;

    // from xsnoed... ???
    r->PmtEventRecordInfo = PMT_EVR_RECTYPE | PMT_EVR_NOT_MC | PMT_EVR_ZDAB_VER;

    r->DataType = PMT_EVR_DATA_TYPE;
    r->DaqStatus = ds->GetSubRunID();
    r->CalPckType = PMT_EVR_PCK_TYPE | PMT_EVR_CAL_TYPE;

    PMTCal* cal_pmt = static_cast<PMTCal*>(NULL);

    uint32_t* sub_header = &(r->CalPckType);

    // must set the size of this sub-field before calling AddSubField()
    // (from the sub-field header to the end)
    *sub_header |= (static_cast<uint32_t*>(fec_word + npmts * 3) - sub_header);

    // calibrated hit information
    if (ev.GetCalPMTs().GetAllCount() > 0) {
        PZdabFile::AddSubField(&sub_header, SUB_TYPE_CALIBRATED,
                               npmts * sizeof(PMTCal));
        cal_pmt = reinterpret_cast<PMTCal*>(sub_header + 1);
    }

    // add monte carlo data
    if (ds->MCExists()) {
        RAT::DS::MC& mc = ds->GetMC();
        // if tracking is enabled, add track steps. otherwise, just add
        // the primary vertex MCParticles
        if (mc.GetMCTrackCount() > 0) {
            // build table to look up sequential vertex id by track & step index
            std::vector<unsigned int> mcTrackIDs = mc.GetMCTrackIDs();
            std::map<std::pair<int, int>, int> step_to_vertex_id;
            std::map<int, int> track_id_to_index;
            int nvertices = 0;
            for (size_t itrack=0; itrack < mcTrackIDs.size(); itrack++) {
                RAT::DS::MCTrack& track = mc.GetMCTrack(mcTrackIDs.at(itrack));
                track_id_to_index[track.GetTrackID()] = itrack;
                for (size_t s=0; s<track.GetMCTrackStepCount(); s++) {
                    step_to_vertex_id[std::make_pair(itrack, s)] = nvertices++;
                }
            }
            //assert(nvertices == step_to_vertex_id.size());

            if (nvertices > 0) {
                PZdabFile::AddSubField(&sub_header,
                  SUB_TYPE_MONTE_CARLO,
                  sizeof(MonteCarloHeader)+nvertices*sizeof(MonteCarloVertex));

                // get pointer to start of monte-carlo data and vertices
                MonteCarloHeader* mc_hdr =
                  reinterpret_cast<MonteCarloHeader*>(sub_header + 1);
                MonteCarloVertex* mc_vtx =
                  reinterpret_cast<MonteCarloVertex*>(mc_hdr + 1);

                // fill in the monte carlo data
                mc_hdr->nVertices = nvertices;

                int vertex_id = 0;
                for (size_t itrack=1; itrack <= mc.GetMCTrackCount(); itrack++) {
                    RAT::DS::MCTrack& track = mc.GetMCTrack(itrack);
                    for (size_t istep = 0; istep<track.GetMCTrackStepCount();istep++) {
                        RAT::DS::MCTrackStep& step=track.GetMCTrackStep(istep);

                        int pdgcode = track.GetPDGCode();
                        float momentum = step.GetMomentum().Mag();

                        // convert kinetic to total energy
                        double mass = 0;
                        std::map<int, float>::const_iterator it_mass =
                          ratzdab::pdg_to_mass.find(pdgcode);
                        if (it_mass != ratzdab::pdg_to_mass.end()) {
                            mass = it_mass->second;
                        }
                        else {
                            warn << "ratzdab::pack::event: No mass available "
                                 "for PDG code " << pdgcode << ", using zero."
                                 << newline;
                        }

                        // all MeV, c=1
                        mc_vtx->energy=TMath::Sqrt(mass*mass+momentum*momentum);

                        // convert mm -> cm
                        mc_vtx->x = step.GetPosition()[0] / 10;
                        mc_vtx->y = step.GetPosition()[1] / 10;
                        mc_vtx->z = step.GetPosition()[2] / 10;

                        // u, v, w are direction cosines
                        mc_vtx->u = step.GetMomentum()[0] / momentum;
                        mc_vtx->v = step.GetMomentum()[1] / momentum;
                        mc_vtx->w = step.GetMomentum()[2] / momentum;

                        // convert from PDG to SNOMAN code
                        int code = 0;  // 0 indicates unknown particle
                        std::map<int, int>::const_iterator it_code =
                          ratzdab::pdg_to_snoman.find(pdgcode);
                        if (it_code != ratzdab::pdg_to_snoman.end()) {
                            code = it_code->second;
                        }
                        mc_vtx->particle = code;

                        // interaction code; FIXME: 100 = start
                        mc_vtx->int_code = 100;

                        // parent is previous step, or if this is the first
                        // step in the track, look up parent track and
                        // fuzzy match to determine which step
                        int parent = vertex_id - 1;
                        if (istep == 0) {
                            int parent_track_id = track.GetParentID();
                            std::map<int, int>::iterator it_ptid =
                              track_id_to_index.find(parent_track_id);
                            if (it_ptid == track_id_to_index.end()) {
                                parent = -1;
                            }
                            else {
                                RAT::DS::MCTrack& parent_track =
                                  mc.GetMCTrack(it_ptid->second);

                                assert(parent_track.GetMCTrackStepCount() > 0);
                                double closest_distance =
                                  (parent_track.GetMCTrackStep(0).GetPosition()
                                  - step.GetPosition()).Mag();
                                int closest_step_id = 0;
                                for (size_t jstep = 0;
                                     jstep<parent_track.GetMCTrackStepCount();
                                     jstep++) {
                                    RAT::DS::MCTrackStep& parent_step =
                                      parent_track.GetMCTrackStep(jstep);
                                    double dist = (parent_step.GetPosition() -
                                                   step.GetPosition()).Mag();
                                    if (dist < closest_distance) {
                                        closest_distance = dist;
                                        closest_step_id = jstep;
                                    }
                                }
                                parent = step_to_vertex_id[std::make_pair(
                                         it_ptid->second, closest_step_id)];
                            }
                        }

                        mc_vtx->parent = parent;

                        mc_vtx->time = step.GetGlobalTime();
                        mc_vtx->flags = 0;

                        vertex_id++;
                        mc_vtx++;
                    }
                }
            }
        }
        else {
            int nvertices = mc.GetMCParticleCount();
            if (nvertices > 0) {
                PZdabFile::AddSubField(&sub_header,
                  SUB_TYPE_MONTE_CARLO,
                  sizeof(MonteCarloHeader)+nvertices*sizeof(MonteCarloVertex));

                // get pointer to start of monte-carlo data and vertices
                MonteCarloHeader* mc_hdr =
                  reinterpret_cast<MonteCarloHeader*>(sub_header + 1);
                MonteCarloVertex* mc_vtx =
                  reinterpret_cast<MonteCarloVertex*>(mc_hdr + 1);

                // fill in the monte carlo data values
                mc_hdr->nVertices = nvertices;
                for (int i=0; i<nvertices; i++, mc_vtx++) {
                    RAT::DS::MCParticle& p = mc.GetMCParticle(i);
                    int pdgcode = p.GetPDGCode();
                    float momentum = p.GetMomentum().Mag();

                    // convert kinetic to total energy
                    double mass = 0;
                    std::map<int, float>::const_iterator it_mass =
                      ratzdab::pdg_to_mass.find(pdgcode);
                    if (it_mass != ratzdab::pdg_to_mass.end()) {
                        mass = it_mass->second;
                    }
                    else {
                        warn << "ratzdab::pack::event: No mass available for "
                          "PDG code " << pdgcode << ", using zero." << newline;
                    }

                    // all MeV, c=1
                    mc_vtx->energy = TMath::Sqrt(mass*mass + momentum*momentum);

                    // convert mm -> cm
                    mc_vtx->x = p.GetPosition()[0] / 10;
                    mc_vtx->y = p.GetPosition()[1] / 10;
                    mc_vtx->z = p.GetPosition()[2] / 10;

                    // u, v, w are direction cosines
                    mc_vtx->u = p.GetMomentum()[0] / momentum;
                    mc_vtx->v = p.GetMomentum()[1] / momentum;
                    mc_vtx->w = p.GetMomentum()[2] / momentum;

                    // convert from PDG to SNOMAN code
                    int code = 0;  // 0 indicates unknown particle
                    std::map<int, int>::const_iterator it_code =
                      ratzdab::pdg_to_snoman.find(pdgcode);
                    if (it_code != ratzdab::pdg_to_snoman.end()) {
                        code = it_code->second;
                    }
                    mc_vtx->particle = code;

                    // interaction code; FIXME: 100 = start
                    mc_vtx->int_code = 100;

                    mc_vtx->parent = -1;  // p->GetIndex() in QSNO
                    mc_vtx->time = p.GetTime();
                    mc_vtx->flags = 0;
                }
            }
        }
    }

    // add all available fits
    const std::vector<std::string>& fitNames = ev.GetFitNames();
    for( size_t iFit = 0; iFit < fitNames.size(); iFit++ ) {
        PZdabFile::AddSubField(&sub_header, SUB_TYPE_FIT, sizeof(FittedEvent));

        // get pointer to start of fit data
        FittedEvent* fit = reinterpret_cast<FittedEvent*>(sub_header + 1);
        const RAT::DS::FitResult& rfit = ev.GetFitResult( fitNames[iFit] );

        if (rfit.GetVertexCount() == 0) {
            continue;
        }

        // FIXME can zdab handle >1 vertices?
        RAT::DS::FitVertex rv = rfit.GetVertex(0);

        if (rv.ContainsPosition()) {
            fit->x = rv.GetPosition()[0]/10;
            fit->y = rv.GetPosition()[1]/10;
            fit->z = rv.GetPosition()[2]/10;
        }
        else {
            fit->x = -9999;
            fit->y = -9999;
            fit->z = -9999;
        }

        if (rv.ContainsDirection()) {
            fit->u = rv.GetDirection()[0];
            fit->v = rv.GetDirection()[1];
            fit->w = rv.GetDirection()[2];
        }
        else {
            fit->u = 0;
            fit->v = 0;
            fit->w = 0;
        }

        if (rv.ContainsTime()) {
            fit->time = rv.GetTime();
        }
        else {
            fit-> time = -9999;
        }

        try{
          fit->quality = rfit.GetFOM(fitNames[iFit]);
        }
        catch(std::out_of_range &o){
          fit->quality = 0;
        }

        fit->npmts = npmts;
        fit->spare = 0;
        char buff[256];
        snprintf(buff, 256, "%s", fitNames[iFit].c_str());
        memset(fit->name, 0, 32);
        strncpy(fit->name, buff, 31);  // copy fit name (31 chars max)
    }

    // caen digitizer waveforms
    if (caen_length) {
        PZdabFile::AddSubField(&sub_header, SUB_TYPE_CAEN, caen_length);
        uint32_t* caen = sub_header + 1;

        RAT::DS::Digitiser& d = ev.GetDigitiser();
        uint16_t magic = d.GetDataFormat();  // FIXME ???

        uint16_t channel_mask = 0;
        for (int i=0; i<8; i++)
            if(d.ExistsWaveform(i))
              channel_mask |= (1 << i);

        uint16_t pattern = d.GetIOPins();
        uint16_t pack_flag = d.GetBit24();
        uint32_t word_count = d.GetNWords();
        uint32_t board_id = d.GetBoardID();
        uint32_t event_id = d.GetEventID();
        uint32_t trigger_time = d.GetTrigTime();

        // pack header
        *(caen++) = ((magic & 0xf) << 28) | word_count;
        *(caen++) = (board_id << 25) | ((pack_flag & 1) << 24) |
                    ((pattern & 0xffff) << 8) | (channel_mask & 0xff);
        *(caen++) = event_id & 0xffffff;
        *(caen++) = trigger_time;

        // copy samples
        uint16_t* psample = reinterpret_cast<uint16_t*>(caen);
        for (int i = 0; i < 8; i++) {
            if (d.ExistsWaveform(i)) {
                const std::vector<uint16_t> s = d.GetWaveform(i);
                for (unsigned int j = 0; j < s.size(); j++)
                    *(psample++) = s[j];
            }
        }
    }

    // trigger
    uint32_t gtid = ev.GetGTID();
    uint32_t trigger = ev.GetTrigType();
    uint32_t error = ev.GetTrigError();

    uint32_t peak = ev.GetESumPeak() & 0x3ff;
    uint32_t inte = ev.GetESumIntegral() & 0x3ff;
    uint32_t diff = ev.GetESumDerivative() & 0x3ff;

    *(mtc_word++) = static_cast<uint32_t>(ev.GetClockCount10());
    *(mtc_word++) = ((ev.GetClockCount10() >> 32) & 0x7fffff) |
                    ((ev.GetClockCount50() & 0x7ff) << 21);
    *(mtc_word++) = static_cast<uint32_t>(ev.GetClockCount50() >> 11);
    *(mtc_word++) = ((trigger & 0x000000ff) << 24) | (gtid & 0xffffff);
    *(mtc_word++) = ((trigger & 0x07ffff00) >> 8) | (peak << 19) | (diff << 29);
    *(mtc_word++) = ((error & 0x3fff) << 17) | (inte << 7) | (diff >> 3);

    // pmts
    for (int i=0; i<npmts; i++) {
        const RAT::DS::PMTUncal& pmt = ev.GetUncalPMTs().GetAllPMT(i);

        uint32_t crate_id = static_cast<uint32_t>((pmt.GetID() >> 9) & 0x1f);
        uint32_t card_id = static_cast<uint32_t>((pmt.GetID() >> 5) & 0xf);
        uint32_t channel_id = static_cast<uint32_t>(pmt.GetID() & 0x1f);
        uint32_t cell_id = static_cast<uint32_t>(pmt.GetCellID());

        *(fec_word++) = (card_id << 26) | (crate_id << 21) |
                        (channel_id << 16) | (gtid & 0xffffUL);

        *(fec_word++) = ((static_cast<uint32_t>(pmt.GetQHS()) ^ 0x800) << 16) |
             (cell_id << 12) | ((static_cast<uint32_t>(pmt.GetQLX()) ^ 0x800));

        *(fec_word++) = ((gtid & 0x00f00000UL) << 8) |
            ((static_cast<uint32_t>(pmt.GetTime()) ^ 0x800) << 16) |
            ((gtid & 0x000f0000UL) >> 4) |
            ((static_cast<uint32_t>(pmt.GetQHL()) ^ 0x800));

        // fill in calibrated PmtEventRecord entries if available
        if (cal_pmt) {
            const RAT::DS::PMTCal& rpmtcal = ev.GetCalPMTs().GetAllPMT(i);
            cal_pmt->tac = rpmtcal.GetTime();
            cal_pmt->qhs = rpmtcal.GetQHS();
            cal_pmt->qhl = rpmtcal.GetQHL();
            cal_pmt->qlx = rpmtcal.GetQLX();
            cal_pmt++;
        }
    }

    return r;
}

////////////////////////////////////////////////////////////////////////
// pack::rhdr
RunRecord* pack::rhdr(RAT::DS::Run* const run) {
    int i;

    RunRecord* rr = new RunRecord;

    rr->Date = run->GetDate();
    rr->Time = run->GetTime();
    rr->SubRunNumber = run->GetSubRunID();
    rr->RunNumber = run->GetRunID();
    rr->CalibrationTrialNumber = 0; /* RAT doesn't currently store this. */
    rr->SourceMask = 0x0;//FIXME PHIL from db run->GetSrcMask();
    rr->RunMask = run->GetRunType();
    rr->GTCrateMsk = 0x0;//FIXME PHIL from db run->GetCrateMask();
    rr->FirstGTID = run->GetFirstGTID();
    rr->ValidGTID = run->GetValidGTID();

    for (i = 0; i < 8; i++) {
        rr->Spares[i] = 0;
    }

    return rr;
}

////////////////////////////////////////////////////////////////////////
// pack::trig
TriggerInfo* pack::trig(RAT::DS::TrigHeader* const triginfo,
                        const uint32_t gtid) {
    TriggerInfo* ti = new TriggerInfo;

    ti->TriggerMask = triginfo->GetTrigMask();

    ti->n100lo  = triginfo->GetTriggerThreshold(0);
    ti->n100med = triginfo->GetTriggerThreshold(1);
    ti->n100hi  = triginfo->GetTriggerThreshold(2);
    ti->n20     = triginfo->GetTriggerThreshold(3);
    ti->n20lb   = triginfo->GetTriggerThreshold(4);
    ti->esumlo  = triginfo->GetTriggerThreshold(5);
    ti->esumhi  = triginfo->GetTriggerThreshold(6);
    ti->owln    = triginfo->GetTriggerThreshold(7);
    ti->owlelo  = triginfo->GetTriggerThreshold(8);
    ti->owlehi  = triginfo->GetTriggerThreshold(9);

    ti->n100lo_zero  = triginfo->GetTriggerZeroOffset(0);
    ti->n100med_zero = triginfo->GetTriggerZeroOffset(1);
    ti->n100hi_zero  = triginfo->GetTriggerZeroOffset(2);
    ti->n20_zero     = triginfo->GetTriggerZeroOffset(3);
    ti->n20lb_zero   = triginfo->GetTriggerZeroOffset(4);
    ti->esumlo_zero  = triginfo->GetTriggerZeroOffset(5);
    ti->esumhi_zero  = triginfo->GetTriggerZeroOffset(6);
    ti->owln_zero    = triginfo->GetTriggerZeroOffset(7);
    ti->owlelo_zero  = triginfo->GetTriggerZeroOffset(8);
    ti->owlehi_zero  = triginfo->GetTriggerZeroOffset(9);

    ti->PulserRate = triginfo->GetPulserRate();
    ti->ControlRegister = triginfo->GetMTC_CSR();
    ti->reg_LockoutWidth = triginfo->GetLockoutWidth();
    ti->reg_Prescale = triginfo->GetPrescaleFrequency();
    ti->GTID = gtid;

    return ti;
}

////////////////////////////////////////////////////////////////////////
// pack::eped
EpedRecord* pack::eped(RAT::DS::ECAHeader* const epedinfo,
                       const uint32_t gtid) {
    EpedRecord* er = new EpedRecord;

    er->ped_width = epedinfo->GetChargePedestalWidth();
    er->ped_delay_coarse = epedinfo->GetGTDelayCoarse();
    er->ped_delay_fine = epedinfo->GetGTDelayFine();
    er->qinj_dacsetting = epedinfo->GetChargePedestalAmplitude();
    er->halfCrateID = epedinfo->GetPatternID();
    er->CalibrationType = epedinfo->GetCalType();
    er->GTID = gtid;
    er->Flag = 0;

    return er;
}

}  // namespace ratzdab