source: trunk/MagicSoft/Mars/mraw/MRawRunHeader.cc@ 8776

/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2007
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MRawRunHeader
//
// Root storage container for the RUN HEADER information
//
//
//  RAW DATA FORMAT VERSION
//  =======================
//
//  Format Version 9:
//  -----------------
//   + fNumEventsRead;
//   + fSamplingFrequency
//   - fFreqSampling
//   + fFadcResolution;
//   - fNumSignificantBits
//
//  Format Version 8:
//  -----------------
//   + fNumBytesPerSample;
//   + fFreqSampling;
//   + fNumSignificantBits
//   * changes in MRawCrateHeader
//
//  Format Version 7:
//  -----------------
//   - unused
//
//  Format Version 6:
//  -----------------
//   + added CameraVersion
//   + added TelescopeNumber
//   + added ObservationMode
//   + added dummies for TelescopeRa/Dec
//
//  Format Version 5:
//  -----------------
//   - now the sub millisecond information of the time is valid and decoded
//     which enhances the precision from 51.2us to 200ns
//
//  Format Version 4:
//  -----------------
//   - added support for pixels with negative IDs
//
//  Format Version 3:
//  -----------------
//   - ???
//
//  Format Version 2:
//  -----------------
//   - removed mjd from data
//   - added start time
//   - added stop  time
//
//
//  MRawRunHeader CLASS VERSION
//  ===========================
//
//  Format Version 7:
//  -----------------
//   - added fNumEventsRead;
//
//  Format Version 6:
//  -----------------
//   - added fNumBytesPerSample;
//   - added fSamplingFrequency;
//   - added fFadcResolution;
//
//  Class Version 5:
//  -----------------
//   - for compatibility with newer camera versions
//
//  Class Version 4:
//  -----------------
//   - added fCameraVersion
//   - added fTelescopeNumber
//   - changed length of fProjectName to 101
//   - changed length of fSourceName  to 81
//
//  Class Version 3:
//  ----------------
//   - enhanced SourceName and ProjectName by one character, because
//     without telling us the guranteed trailing \0-character has
//     skipped
//
//  Class Version 2:
//  ----------------
//   - removed fMJD, fYear, fMonth, fDay
//   - added fRunStart
//   - added fRunStop
// 
//  Class Version 1:
//  ----------------
//   - first implementation
//
////////////////////////////////////////////////////////////////////////////
#include "MRawRunHeader.h"

#include <fstream>
#include <iomanip>

#include "MLog.h"
#include "MLogManip.h"

#include "MArrayS.h"

ClassImp(MRawRunHeader);

using namespace std;

const UShort_t MRawRunHeader::kMagicNumber      = 0xc0c0;
const Byte_t   MRawRunHeader::kMaxFormatVersion =      9;

// --------------------------------------------------------------------------
//
// Default constructor. Creates array which stores the pixel assignment.
//
//
MRawRunHeader::MRawRunHeader(const char *name, const char *title) : fPixAssignment(NULL)
{
    fName  = name  ? name  : "MRawRunHeader";
    fTitle = title ? title : "Raw Run Header Information";

    fPixAssignment = new MArrayS(0);

    // Remark: If we read old MC data from a root file which do not
    // yet contain one of these variable, the value given here is
    // the default. Do not mix files with and without a value if the
    // files with the value do not match the default!
    fFormatVersion=0;
    fSoftVersion=0;
    fTelescopeNumber=1;
    fCameraVersion=1;
    fFadcType=1;
    fRunType=kRTNone;  // use 0xffff for invalidation, 0 means: Data run
    fRunNumber=0;
    memset(fProjectName,     0, 101);
    memset(fSourceName,      0,  81);
    memset(fObservationMode, 0,  61);
    fSourceEpochChar[0]=0;
    fSourceEpochDate=0;
    fNumCrates=0;
    fNumPixInCrate=0;
    fNumSamplesLoGain=0;
    fNumSamplesHiGain=0;
    fNumEvents=0;
    fNumBytesPerSample=1;
    fSamplingFrequency=300;
    fFadcResolution=8;
}

// --------------------------------------------------------------------------
//
// Destructor. Deletes the 'pixel-assignment-array'
//
MRawRunHeader::~MRawRunHeader()
{
    delete fPixAssignment;
}

// --------------------------------------------------------------------------
//
// Swap the assignment of the pixels with hardware id id0 and id1
//
Bool_t MRawRunHeader::SwapAssignment(Short_t id0, Short_t id1)
{
    const Int_t n = fPixAssignment->GetSize();

    // Look-up-table
    UShort_t *lut = fPixAssignment->GetArray();

    // Search for one of the hardware indices to get exchanged
    int i;
    for (i=0; i<n; i++)
        if (lut[i]==id0 || lut[i]==id1)
            break;

    // Check if one of the two pixels were found
    if (i==n)
    {
        *fLog << warn << "WARNING - Assignment of pixels with hardware ID " << id0 << " and " << id1;
        *fLog << " should be exchanged, but none were found." << endl;
        return kTRUE;
    }

    // Store first index
    const Int_t idx0 = i;

    // Search for the other hardware indices to get exchanged
    for (i++; i<n; i++)
        if (lut[i]==id0 || lut[i]==id1)
            break;

    // Check if the second pixel was found
    if (i==n)
    {
        *fLog << err << "ERROR - Assignment of pixels with hardware ID " << id0 << " and " << id1;
        *fLog << " should be exchanged, but only one was found." << endl;
        return kFALSE;
    }

    const Int_t idx1 = i;

    const Short_t p0 = lut[idx0];
    const Short_t p1 = lut[idx1];

    lut[idx0] = p1;
    lut[idx1] = p0;

    *fLog << inf << "Assignment of pixels with hardware ID " << id0 << " and " << id1 << " exchanged." << endl;

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// This implements a fix of the pixel assignment before 25.9.2005
//
// From magic_online (the pixel numbers are hardware indices):
// --------------------------------------------------------------------------
//   From: Florian Goebel <fgoebel@mppmu.mpg.de>
//   Date: Sun Sep 25 2005 - 05:13:19 CEST
//   
//   [...]
//   - problem 2: pixels were swaped
//   - cause: opical fibers were wrongly connected to receiver board
//      554 <-> 559
//      555 <-> 558
//      556 <-> 557
//   - action: reconnect correctly
//   - result: ok now
//   - comment: I don't know when this error was introduced, so backward
//      correction of the data is difficult.
//      Fortunately the effect is not too large since the affected 6 pixels are
//      on the outermost edge of the camera
//      Since this board has special pixels the optical fibers have been
//      unplugged several times.
//   !!!!! Whenever you unplug and reconnect optical fibers make really !!!!!
//   !!!!! sure you connect them back correctly. !!!!!
//   !!!!! Always contact me before you do so and if you have any doubts !!!!!
//   !!!!! how to reconnect the fibers ask me. !!!!!
//      These swapped pixels have only been found by chance when doing the
//      flatfielding.
//   [...]
//
// --------------------------------------------------------------------------
//
// MAGIC runbook CC_2006_04_22_22_28_52.rbk
//
// [2006-04-22 23:14:13]
//
// [...]
// Found 2 pairs of swapped pixels.
// We corrected swapped pixels 54<->55 in the receiver boards. We probably
// swapped today in the camera.
// We did not correct 92<-<93 which are likely swapped. Will check and correct
// tomorrow.
//
// ---
//
// comments:
// - 54<->55 were corrected but have not been swapped, hence they are swapped
//   since then (run 88560 ok, run 88669 swapped; between them mostly dummy and
//   test runs)
// - 92<->93 are never swapped, always ok.
//
// --------------------------------------------------------------------------
//
// MAGIC runbook CC_2006_08_28_19_40_18.rbk
//
// [2006-08-28 23:09:07]
// While doing a flatfielding we have found out that the signals for pixels
// 119 and 120 were swapped. We have fixed it by exchanging the corresponding
// fibers at the input of the receivers (not before the splitters!).
//
// ---
//
// MAGIC runbook CC_2006_08_29_15_19_14.rbk
//
// [2006-08-29 16:43:09]
// In the last hours we have found out and fixed a good number of pixels which
// were swapped: 119-120, 160-161-162 and 210-263. According to Florian,
// 160-161-162 and 210-263 were swapped since November 2005.
//
// ---
//
// mail Florian Goebel (08/30/2006 03:13 PM):
//
// As far as I can tell pixels 161 and 162 as well as 263 and 210 were
// swapped in the trigger. This leads to some inefficiency of the trigger.
// However, they were not swapped in the readout. So, you don't have to
// correct anything in the data for these pixels.
//
// ---
//
// comments:
// - 119-120 swapped between run 93251 (not swapped) and 93283 (swapped)
//   (only testruns between these runs)
//   corrected since run 99354 (== runbook [2006-08-28 23:09:07])
// - 160 never swapped
// - 161-162 were only swapped in the trigger, but nevertheless were
//   "corrected" also in the signal. Hence signal swapped since 99354
//
// --------------------------------------------------------------------------

Bool_t MRawRunHeader::FixAssignment()
{
    if (fRunNumber>=53300 && fRunNumber<=68754)
    {
        if (!SwapAssignment(554, 559))
            return kFALSE;
        if (!SwapAssignment(555, 558))
            return kFALSE;
        if (!SwapAssignment(556, 557))
            return kFALSE;
    }

    if (fRunNumber>=93283 && fRunNumber<99354)
    {
        if (!SwapAssignment(119, 120))
            return kFALSE;
    }

    if (fRunNumber>=99354 && fRunNumber<=101789)
    {
        if (!SwapAssignment(161, 162))
            return kFALSE;
        if (!SwapAssignment(210, 263))
            return kFALSE;
    }

    if (fRunNumber>=88669)
    {
        if (!SwapAssignment(54, 55))
            return kFALSE;
    }

    if (fRunNumber>=200000)
    {
        if (!SwapAssignment(428, 429))
            return kFALSE;
    }

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Read in one run header from the binary file
//
Bool_t MRawRunHeader::ReadEvt(istream& fin)
{
    //
    // read one RUN HEADER from the input stream
    //
    fMagicNumber = 0;

    fin.read((char*)&fMagicNumber, 2);          // Total=2

    //
    // check whether the the file has the right file type or not
    //
    if (fMagicNumber != kMagicNumber && fMagicNumber != kMagicNumber+1)
    {
        *fLog << err << "ERROR - Wrong Magic Number (0x" << hex << fMagicNumber << "): Not a Magic File!" << endl;
        return kFALSE;
    }

    if (fMagicNumber == kMagicNumber+1)
        *fLog << warn << "WARNING - This file maybe broken (0xc0c1) - DAQ didn't close it correctly!" << endl;

    Byte_t dummy[16];

    // ----- File format version -----
    fin.read((char*)&fFormatVersion, 2);     // Total=4
    if (fFormatVersion>kMaxFormatVersion)
    {
        *fLog << err << "ERROR - File format V" << fFormatVersion << " not implemented!" << endl;
        return kFALSE;
    }

    if (fFormatVersion==7)
    {
        *fLog << err << "ERROR - File format V7 was for testing only and is not correctly implemented!" << endl;
        return kFALSE;
    }

    // ----- DAQ software format version -----
    fin.read((char*)&fSoftVersion, 2);     // Total=6


    fFadcType = 1;
    if (fFormatVersion>7)
        fin.read((char*)&fFadcType, 2);

    // ----- Camera geometry and telescope number -----
    fCameraVersion   = 1;
    fTelescopeNumber = 1;
    if (fFormatVersion>5)
    {
        fin.read((char*)&fCameraVersion,   2); // (+2)
        fin.read((char*)&fTelescopeNumber, 2); // (+2)
    }

    // ----- Run information -----
    fin.read((char*)&fRunType,     2); // Total=8

    fin.read((char*)&fRunNumber,   4); // Total=12
    fin.read((char*)&fProjectName, fFormatVersion>5?100:22); // Total=34  (+78)
    fin.read((char*)&fSourceName,  fFormatVersion>5? 80:12); // Total=46  (+58)

    if (fFormatVersion>5)
        fin.read((char*)fObservationMode, 60); // (+60)
    // Maybe we should set fObservationMode to something
    // in case of fFormatVersion<6

    // ----- Source position -----
    fin.seekg(fFormatVersion>5 ? 16 : 8, ios::cur);
    /*
    if (fFormatVersion>5)
    {
        fin.read((char*)&fSourceRa,  4);  // F32 SourceRA;  Total=48
        fin.read((char*)&fSourceDec, 4);  // F32 SourceDEC; Total=52
    }
    fin.read((char*)&fTelescopeRa,  4);   // F32 TelescopeRA;  (+4)
    fin.read((char*)&fTelescopeDec, 4);   // F32 TelescopeDEC; (+4)
    */

    // Maybe we should set these to something
    // in case of fFormatVersion<6
    fin.read((char*)&fSourceEpochChar,  2);     // Total=56
    fin.read((char*)&fSourceEpochDate,  2);     // Total=58

    // ----- Old Start time -----
    if (fFormatVersion<2)                       // Total += 10
    {
        UShort_t y, m, d;
        fin.seekg(4, ios::cur);    // Former fMJD[4],
        fin.read((char*)&y,    2); // Former fDateYear[2]
        fin.read((char*)&m,    2); // Former fDateMonth[2]
        fin.read((char*)&d,    2); // Former fDateDay[2]
        fRunStart.Set(y, m, d, 0, 0, 0, 0);
    }

    // ----- Data Geometry -----

    fin.read((char*)&fNumCrates,        2); // MUX: number of channels
    fin.read((char*)&fNumPixInCrate,    2); // MUX: number of pix in channel
    fin.read((char*)&fNumSamplesLoGain, 2); // MUX: dummy (must be 0 for MUX data)
    fin.read((char*)&fNumSamplesHiGain, 2); // MUX: Number of samples per pixel

    if (fFormatVersion>8)
        fin.read((char*)dummy, 4); // 2xU16 (NumSamplesRemovedHead and NumSamplesRemovedTail)

    // ----- Number of events -----
    fin.read((char*)&fNumEvents, 4);     // Total=70

    if (fFormatVersion>8)
        fin.read((char*)&fNumEventsRead, 4);     // Total=70

    // New in general features: (should they be included in new MAGIC1 formats, too?)
    fNumBytesPerSample  = 1;      // 2 for MUX DATA
    fSamplingFrequency  = 300;
    fFadcResolution = 8;
    if (fFormatVersion>7)
    {
        fin.read((char*)&fNumBytesPerSample,  2);
        fin.read((char*)&fSamplingFrequency,  2); // [MHz], 2000 for MuxFadc
        fin.read((char*)&fFadcResolution, 1); // nominal resolution [# Bits], 10 for MuxFadc

        if (fNumBytesPerSample!=2)
        {
            *fLog << err << "ERROR - " << fNumBytesPerSample << " bytes per sample are not supported!" << endl;
            return kFALSE;
        }
    }

    // ----- Start/Stop time -----
    if (fFormatVersion>1)
    {
        fRunStart.ReadBinary(fin);              // Total += 7
        fRunStop.ReadBinary(fin);               // Total += 7
    }

    //
    // calculate size of array, create it and fill it
    //
    const Int_t nPixel = fNumCrates*fNumPixInCrate;
    fPixAssignment->Set(nPixel);

    // ----- Pixel Assignement -----
    fin.read((char*)fPixAssignment->GetArray(), nPixel*2);

    if (fFormatVersion<7)
        fin.read((char*)&dummy, 16);

    return FixAssignment();
}

// --------------------------------------------------------------------------
//
// Return the run type as string ("Data", "Pedestal", ...), for example
// to print it as readable text.
//
const char *MRawRunHeader::GetRunTypeStr() const
{
    switch (fRunType)
    {
    case kRTData:
        return "Data";
    case kRTPedestal:
        return "Pedestal";
    case kRTCalibration:
        return "Calibration";
    case kRTPointRun:
        return "Point-Run";
    case kRTMonteCarlo:
        return "Monte Carlo";
    case kRTNone:
        return "<none>";
    default:
        return "<unknown>";
    }
}

// --------------------------------------------------------------------------
//
// print run header information on *fLog. The option 'header' supresses
// the pixel index translation table.
//
void MRawRunHeader::Print(Option_t *t) const
{
    *fLog << all << endl;
    *fLog << "MagicNumber:  0x" << hex << fMagicNumber << " - ";
    switch (fMagicNumber)
    {
    case kMagicNumber:   *fLog << "OK";               break;
    case kMagicNumber+1: *fLog << "File not closed!"; break;
    default:             *fLog << "Wrong!";           break;
    }
    *fLog << endl;
    *fLog << "Versions:     " << dec << "Format=" << fFormatVersion << "  ";
    *fLog << "Software=" << fSoftVersion << "  ";
    if (fFormatVersion>5)
        *fLog << "Camera=" << fCameraVersion;
    *fLog << endl;
    if (fFormatVersion>5)
        *fLog << "Telescope:    " << fTelescopeNumber << endl;
    if (fFormatVersion>7)
    {
        *fLog << "FadcType:     " << fFadcType << " (";
        switch (fFadcType)
        {
        case 1: *fLog << "Siegen"; break;
        case 2: *fLog << "MUX"; break;
        default: *fLog << "unknown";
        }
        *fLog << ")" << endl;
    }
    *fLog << "RunNumber:    " << fRunNumber << " (Type=" << GetRunTypeStr() << ")" << endl;
    *fLog << "ProjectName: '" << fProjectName << "'" << endl;
    if (fFormatVersion>5)
        *fLog << "Observation: '" << fObservationMode << "'" << endl;
    *fLog << "Source:      '" << fSourceName << "' " << "  ";
    *fLog << fSourceEpochChar << dec << fSourceEpochDate << endl;
    *fLog << "Run Start:    " << fRunStart << endl;
    *fLog << "Run Stop:     " << fRunStop << endl;
    *fLog << "Crates:       " << fNumCrates << " x " << fNumPixInCrate << " Pixel/Crate = " << fNumCrates*fNumPixInCrate << " Pixel/Evt" << endl;
    *fLog << "Num Pixels:   " << GetNumNormalPixels() << " (normal) + " << GetNumSpecialPixels() << " (special) = " << GetNumConnectedPixels() << " (total)" << endl;
    if (fFormatVersion>6)
        *fLog << "Sampling:     " << fSamplingFrequency << "MHz with " << (int)fFadcResolution << " significant bits" << endl;
    *fLog << "Samples:      " << fNumSamplesHiGain << "/" << fNumSamplesLoGain << " (hi/lo) * " << fNumBytesPerSample << "B/sample = " << (fNumSamplesLoGain+fNumSamplesHiGain) * fNumCrates * fNumPixInCrate * fNumBytesPerSample/1000 << "kB/Evt" << endl;
    *fLog << "Evt Counter:  " << fNumEvents;
    if (fFormatVersion>8)
        *fLog << " (read=" << fNumEventsRead << ")";
    *fLog << endl;

    if (TString(t).Contains("header", TString::kIgnoreCase))
        return;

    *fLog << inf3 << "Assignment:" << hex << endl;
    for (int i=0; i<GetNumPixel(); i++)
        *fLog << setfill('0') << setw(3) << (*fPixAssignment)[i] << " ";

    *fLog << dec << setfill(' ') << endl;
}

// --------------------------------------------------------------------------
//
// Return the assigned pixel number for the given FADC channel
//
Short_t MRawRunHeader::GetPixAssignment(UShort_t i) const
{
    // FIXME: Do we need a range check here?
    return (Short_t)(*fPixAssignment)[i];
}

// --------------------------------------------------------------------------
//
// Return the number of pixel which are markes as connected in the
// pix assignment (!=0)
//
UShort_t MRawRunHeader::GetNumConnectedPixels() const
{
    const Int_t num = fPixAssignment->GetSize();

    UShort_t rc = 0;
    for (int i=0; i<num; i++)
    {
        if (GetPixAssignment(i)!=0)
            rc++;
    }
    return rc;
}

// --------------------------------------------------------------------------
//
// Return the number of pixel which are markes as connected and so-called
// 'normal' pixels in the pix assignment (>0)
//
UShort_t MRawRunHeader::GetNumNormalPixels() const
{
    const Int_t num = fPixAssignment->GetSize();

    UShort_t rc = 0;
    for (int i=0; i<num; i++)
    {
        if (GetPixAssignment(i)>0)
            rc++;
    }
    return rc;
}

// --------------------------------------------------------------------------
//
// Return the number of pixel which are markes as connected and so-called
// 'special' pixels in the pix assignment (<0)
//
UShort_t MRawRunHeader::GetNumSpecialPixels() const
{
    const Int_t num = fPixAssignment->GetSize();

    UShort_t rc = 0;
    for (int i=0; i<num; i++)
    {
        if (GetPixAssignment(i)<0)
            rc++;
    }
    return rc;
}

// --------------------------------------------------------------------------
//
// Return the maximum id which exists in the pix assignment
//
UShort_t MRawRunHeader::GetMaxPixId() const
{
    const Int_t num = fPixAssignment->GetSize();

    Short_t rc = 0;
    for (int i=0; i<num; i++)
        rc = TMath::Max(GetPixAssignment(i), rc);

    return rc;
}

// --------------------------------------------------------------------------
//
// Return minus th minimum id which exists in the pix assignment
//
UShort_t MRawRunHeader::GetMinPixId() const
{
    const Int_t num = fPixAssignment->GetSize();

    Short_t rc = 0;
    for (int i=0; i<num; i++)
        rc = TMath::Min(GetPixAssignment(i), rc);

    return (UShort_t)-rc;
}

// --------------------------------------------------------------------------
//
// Return the number of pixel in this event.
//
// WARNING: This is the number of pixels stored in this file which is
//          a multiple of the number of pixels per crate and in general
//          a number which is larger than the camera size!
//
//          To know the range of the pixel indices please use the geometry
//          container!
//
UShort_t MRawRunHeader::GetNumPixel() const
{
    return fPixAssignment->GetSize();
}

// --------------------------------------------------------------------------
//
// Returns absolute size in bytes of the run header as read from a raw file.
// This must be done _after_ the header is read, because the header doesn't
// have a fixed size (used in MRawSocketRead)
//
Int_t MRawRunHeader::GetNumTotalBytes() const
{
    switch (fFormatVersion)
    {
    case 1:
        return 80+fNumCrates*fNumPixInCrate*2+16;
    case 2:
    case 3:
    case 4:
    case 5:
        return 84+fNumCrates*fNumPixInCrate*2+16;
    case 6:
        return 84+fNumCrates*fNumPixInCrate*2+16 +4+78+58+60+8;
    case 7:
        return 84+fNumCrates*fNumPixInCrate*2+16 +4+78+58+60+8 +3-16;
    }
    return 0;
}