Index: trunk/MagicSoft/Mars/mcorsika/CorsikaLinkDef.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/CorsikaLinkDef.h (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/CorsikaLinkDef.h (revision 9182) @@ -0,0 +1,11 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ class MCorsikaRunHeader+; +#pragma link C++ class MCorsikaEvtHeader+; +#pragma link C++ class MCorsikaRead+; + +#endif Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9182) @@ -0,0 +1,199 @@ +/* ======================================================================== *\ +! +! * +! * 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 11/2008 +! +! Copyright: Software Development, 2008 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MCorsikaEvtHeader +// +///////////////////////////////////////////////////////////////////////////// +#include "MCorsikaEvtHeader.h" + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +#include "MMcEvt.hxx" +#include "MCorsikaRunHeader.h" + +ClassImp(MCorsikaEvtHeader); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. Create the array to store the data. +// +MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title) +{ + fName = name ? name : "MCorsikaEvtHeader"; + fTitle = title ? title : "Raw Event Header Information"; +} + +// -------------------------------------------------------------------------- +// +// This member function prints all Data of one Event to *fLog. +// +void MCorsikaEvtHeader::Print(Option_t *o) const +{ + *fLog << all; + *fLog << "Event Number: " << dec << fEvtNumber << endl; + *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; + *fLog << "Energy: " << fTotalEnergy << "GeV" << endl; + *fLog << "Starting Altitude: " << fStartAltitude << "g/cm²" << endl; + *fLog << "Number of 1st Target: " << fFirstTargetNum << endl, + *fLog << "Height of 1st Interaction: " << fFirstInteractionHeight/100. << "m" << endl; + *fLog << "Momentum X/Y/Z (GeV/c): " << fMomentumX << "/" << fMomentumY << "/" << fMomentumZ << endl; + *fLog << "Zenith/Azimuth Angle: " << fZd*TMath::RadToDeg() << "°/" << fAz*TMath::RadToDeg() << "°" << endl; + *fLog << "Impact X/Y: " << fX/100. << "m/" << fY/100. << "m (r=" << TMath::Hypot(fX, fY)/100. << "m)" << endl; + *fLog << "Weighted Num Photons: " << fWeightedNumPhotons << endl; + *fLog << endl; +} + +// -------------------------------------------------------------------------- +// +// read the EVTH information from the input stream +// return FALSE if there is no header anymore, else TRUE +// +Int_t MCorsikaEvtHeader::ReadEvt(std::istream &fin, MCorsikaRunHeader &header) +{ + char evth[4]; + fin.read(evth, 4); + if (memcmp(evth, "EVTH", 4)) + { + fin.seekg(-4, ios::cur); + return kFALSE; + } + + Float_t f[273]; + fin.read((char*)&f, 273*4); + + fEvtNumber = TMath::Nint(f[0]); + fParticleID = TMath::Nint(f[1]); + + fTotalEnergy = f[2]; + fStartAltitude = f[3]; + fFirstTargetNum = f[4]; + fFirstInteractionHeight = f[5]; + + // Pointing opposite particle direction + // (x=north, y=west, z=upwards) + //fMomentumX = f[6]; + //fMomentumY = f[7]; + //fMomentumZ = f[8]; + // Pointing along particle direction + // (x=east, y=north, z=upwards) + fMomentumX = -f[7]; + fMomentumY = f[6]; + fMomentumZ = -f[8]; + + fZd = f[9]; + fAz = TMath::Pi()-f[10]; + + const Int_t n = TMath::Nint(f[96]); + if (n!=1) + { + *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; + return kFALSE; + } + + //fX = f[97]; + //fY = f[117]; + + fX = f[117]; + fY = -f[97]; + + if (fEvtNumber==1) + { + header.fZdMin = f[79]; + header.fZdMax = f[80]; + header.fAzMin = 180-f[81]; + header.fAzMax = 180-f[82]; + + header.fViewConeInnerAngle = f[151]; + header.fViewConeOuterAngle = f[152]; + + header.Print(); + } + + fin.seekg(1088-273*4, ios::cur); + + // f[76] Cherenkov flag: + // bit(1) : CERENKOV option compiled in + // bit(2) : IACT option compiled in + // bit(3) : CEFFIC option compiled in + // bit(4) : ATMEXT option compiled in + // bit(5) : ATMEXT option used with refaction enabled + // bit(6) : VOLUMEDET option compiled in + // bit(7) : CURVED option compiled in + // bit(9) : SLATN option compiled in + // 11-21 : table number for externam athmosphere (but<1024) + // f[78] Curved athmosphere? (0=flat, 1=curved) + + // f[80] lower edge of theta in ° + // f[81] upper edge of theta in ° + // f[82] lower edge of phi in ° + // f[83] upper edge of phi in ° + // f[84] cherenkov bunch size + + // f[93] flag for additinal muon information of particle output file + + // f[95] Cherenkov bandwidth lower end in nm + // f[96] Cherenkov bandwidth upper end in nm + + // f[97] Numbr i of uses of each cherenkov event + // f[145] Muon multiple scattering flag + // f[152] !! inner angle of view cone (°) + // f[153] !! outer angle of view cone (°) + // f[156] altitude of horizontal shower axis + + return !fin.eof(); +} + + +// this member function is for reading the event end block + +Bool_t MCorsikaEvtHeader::ReadEvtEnd(std::istream &fin) +{ + //fin.seekg(-1088,ios::cur); + + Float_t f[2]; + fin.read((char*)&f, 2*4); + + const UInt_t evtnum = TMath::Nint(f[0]); + if (evtnum!=fEvtNumber) + { + *fLog << err << "ERROR - Mismatch in stream: Event number in EVTE ("; + *fLog << evtnum << ") doesn't match EVTH (" << fEvtNumber << ")." << endl; + return kERROR; + } + + fWeightedNumPhotons = f[1]; + + fin.seekg(1080,ios::cur); + + return !fin.eof(); +} + Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9182) @@ -0,0 +1,65 @@ +#ifndef MARS_MCorsikaEvtHeader +#define MARS_MCorsikaEvtHeader + +#ifndef MARS_MParContainer +#include "MParContainer.h" +#endif + +#ifndef ROOT_TVector3 +#include +#endif + +// gcc 3.2 +//class ifstream; +#include + +class MCorsikaRunHeader; + +class MCorsikaEvtHeader : public MParContainer +{ +private: + UInt_t fEvtNumber; // Event number + UInt_t fParticleID; // Particle ID (see MMcEvtBasic or CORSIKA manual) + Float_t fTotalEnergy; // [GeV] + + Float_t fStartAltitude; // [g/cm²] + Float_t fFirstTargetNum; // Number of first target if fixed + Float_t fFirstInteractionHeight; // [cm] z coordinate, first intercation height + + Float_t fMomentumX; // [GeV/c] + Float_t fMomentumY; // [GeV/c] + Float_t fMomentumZ; // [GeV/c] + + Float_t fZd; // [rad] Zenith distance + Float_t fAz; // [rad] Azimuth (north=0; west=90) + + Float_t fX; // [cm] Position of telescope on ground x / - impact parameter x + Float_t fY; // [cm] Position of telescope on gorund y / - impact parameter y + + Float_t fWeightedNumPhotons; // weighted number of photons arriving at observation level + +public: + MCorsikaEvtHeader(const char *name=NULL, const char *title=NULL); + + //void Clear(Option_t * = NULL); + void Print(Option_t * = NULL) const; + + TVector3 GetMomentum() const { return TVector3(fMomentumX, fMomentumY, fMomentumZ); } + TVector2 GetImpactPos() const { return TVector2(fX, fY); } + + Float_t GetTotalEnergy() const { return fTotalEnergy; } + Float_t GetFirstInteractionHeight() const { return fFirstInteractionHeight; } + + Float_t GetZd() const { return fZd; } + Float_t GetAz() const { return fAz; } + + Float_t GetX() const { return fX; } + Float_t GetY() const { return fY; } + + Int_t ReadEvt(istream& fin, MCorsikaRunHeader &header); // read in event header block + Bool_t ReadEvtEnd(istream& fin); // read in event end block + + ClassDef(MCorsikaEvtHeader, 1) // Parameter Conatiner for raw EVENT HEADER +}; + +#endif Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc (revision 9182) @@ -0,0 +1,497 @@ +/* ======================================================================== *\ +! +! * +! * 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 +! +! Copyright: MAGIC Software Development, 2000-2007 +! +! +\* ======================================================================== */ + +////////////////////////////////////////////////////////////////////////////// +// +// MCorsikaRead +// +// This tasks reads the raw binary file like specified in the TDAS??? +// and writes the data in the corresponding containers which are +// either retrieved from the parameter list or created and added. +// +// Use SetInterleave() if you don't want to read all events, eg +// SetInterleave(5) reads only each 5th event. +// +// Input Containers: +// -/- +// +// Output Containers: +// MCorsikaRunHeader, MCorsikaEvtHeader, MCorsikaEvtData, MCorsikaCrateArray, MCorsikaEvtTime +// +////////////////////////////////////////////////////////////////////////////// +#include "MCorsikaRead.h" + +#include +#include + +#include + +#include "MLog.h" +#include "MLogManip.h" + +//#include "MZlib.h" +//#include "MTime.h" +#include "MParList.h" +#include "MStatusDisplay.h" + +#include "MCorsikaRunHeader.h" +#include "MCorsikaEvtHeader.h" + +//#include "MPhotonData.h" +#include "MPhotonEvent.h" + +ClassImp(MCorsikaRead); + +using namespace std; + +/* ----------- please don't delete and don't care about (Thomas) ------------ +#define kBUFSZ 64 //1024*1024*64 +#include +class bifstream : public istream, public streambuf +{ +private: + char fBuffer[kBUFSZ]; //! + FILE *fd; + + int sync() + { + memset(fBuffer, 0, kBUFSZ); + return 0; + } + int underflow() + { + int sz=fread(fBuffer, kBUFSZ, 1, fd); + //int sz=fread(fBuffer, 1, kBUFSZ, fd); + setg(fBuffer, fBuffer, fBuffer+kBUFSZ); + + return sz==1 ? *(unsigned char*)fBuffer : EOF;//EOF; + //return sz==kBUFSZ ? *(unsigned char*)fBuffer : EOF;//EOF; + } +public: + bifstream(const char *name) : istream(this) + { + fd = fopen(name, "rb"); + setbuf(fBuffer, kBUFSZ); + } +}; +*/ + +// -------------------------------------------------------------------------- +// +// Default constructor. It tries to open the given file. +// +MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title) + : fRunHeader(0), fEvtHeader(0), fEvent(0), fEvtData(0), fForceMode(kFALSE), + fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0), + fIn(0), fParList(0) +{ + fName = name ? name : "MRead"; + fTitle = title ? title : "Read task to read DAQ binary files"; + + fFileNames = new TList; + fFileNames->SetOwner(); + + if (fname!=NULL) + AddFile(fname); +} + +// -------------------------------------------------------------------------- +// +// Destructor. Delete input stream. +// +MCorsikaRead::~MCorsikaRead() +{ + delete fFileNames; + if (fIn) + delete fIn; +} + +/* +Byte_t MCorsikaRead::IsFileValid(const char *name) +{ + MZlib fin(name); + if (!fin) + return 0; + + Byte_t c[4]; + fin.read((char*)c, 4); + if (!fin) + return 0; + + if (c[0]!=0xc0) + return 0; + + if (c[1]==0xc0) + return 1; + + if (c[1]==0xc1) + return 2; + + return 0; +} +*/ + +// -------------------------------------------------------------------------- +// +// Add a new file to a list of files to be processed, Returns the number +// of files added. (We can enhance this with a existance chack and +// wildcard support) +// +Int_t MCorsikaRead::AddFile(const char *fname, Int_t entries) +{ + TNamed *name = new TNamed(fname, ""); + fFileNames->AddLast(name); + return 1; + +} + +Bool_t MCorsikaRead::ReadEvtEnd() +{ + if (!fRunHeader->SeekEvtEnd(*fIn)) + { + *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl; + if (!fForceMode) + return kFALSE; + } + + if (!fRunHeader->ReadEvtEnd(*fIn)) + { + *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl; + if (!fForceMode) + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// This opens the next file in the list and deletes its name from the list. +// +Int_t MCorsikaRead::OpenNextFile(Bool_t print) +{ + // + // open the input stream and check if it is really open (file exists?) + // + if (fIn) + delete fIn; + fIn = NULL; + + // + // Check for the existance of a next file to read + // + TObject *file = fFileNames->At(fNumFile); + if (!file) + return kFALSE; + + // + // open the file which is the first one in the chain + // + const char *name = file->GetName(); + + const char *expname = gSystem->ExpandPathName(name); + fIn = new fstream(expname); + + const Bool_t noexist = !(*fIn); + if (noexist) + { + *fLog << err << "Cannot open file " << expname << ": "; + *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl; + } + else + { + *fLog << inf << "Open file: '" << name << "'" << endl; + + if (fDisplay) + { + TString txt = GetFileName(); + txt += " @ "; + txt += GetNumExecutions()-1; + fDisplay->SetStatusLine2(txt); + } + } + + delete [] expname; + + if (noexist) + return kERROR; + + fNumFile++; + + // + // Read RUN HEADER (see specification) from input stream + // + if (!fRunHeader->ReadEvt(*fIn)) + return kERROR; + + if (!(*fIn)) + { + *fLog << err << "Error: Accessing file '" << name << "'" << endl; + return kERROR; + } + + const streampos pos = fIn->tellg(); + + if (!ReadEvtEnd()) + return kERROR; + + fIn->seekg(pos, ios::beg); + + + fNumEvents += fRunHeader->GetNumEvents(); + fRunHeader->SetReadyToSave(); + + // + // Print Run Header + // We print it after the first event was read because + // we still miss information which is stored in the event header?!? + //if (print) + // fRunHeader->Print(); + + if (!fParList) + return kTRUE; + + // + // Search for MTaskList + // + MTask *tlist = (MTask*)fParList->FindObject("MTaskList"); + if (!tlist) + { + *fLog << err << dbginf << "MTaskList not found... abort." << endl; + return kERROR; + } + + // + // A new file has been opened and new headers have been read. + // --> ReInit tasklist + // + return tlist->ReInit(fParList) ? kTRUE : kERROR; +} + +// -------------------------------------------------------------------------- +// +// Return file name of current file. +// +TString MCorsikaRead::GetFullFileName() const +{ + const TObject *file = fFileNames->At(fNumFile-1); + return file ? file->GetName() : ""; +} + +// -------------------------------------------------------------------------- +// +// Restart with the first file +// +Bool_t MCorsikaRead::Rewind() +{ + fNumFile=0; + fNumEvents=0; + return OpenNextFile()==kTRUE; +} + +Bool_t MCorsikaRead::CalcNumTotalEvents() +{ + fNumTotalEvents = 0; + + Bool_t rc = kTRUE; + + while (1) + { + switch (OpenNextFile(kFALSE)) + { + case kFALSE: + break; + case kERROR: + rc = kFALSE; + break; + case kTRUE: + if (!ReadEvtEnd()) + { + rc = kFALSE; + break; + } + + fNumTotalEvents += fRunHeader->GetNumEvents(); + continue; + } + break; + } + + if (fIn) + delete fIn; + fIn = NULL; + + return rc; +} + +// -------------------------------------------------------------------------- +// +// The PreProcess of this task checks for the following containers in the +// list: +// MCorsikaRunHeader if not found it is created +// MCorsikaEvtHeader if not found it is created +// MCorsikaEvtData if not found it is created +// MCorsikaCrateArray if not found it is created +// MCorsikaEvtTime if not found it is created (MTime) +// +// If all containers are found or created the run header is read from the +// binary file and printed. If the Magic-Number (file identification) +// doesn't match we stop the eventloop. +// +// Now the EvtHeader and EvtData containers are initialized. +// +Int_t MCorsikaRead::PreProcess(MParList *pList) +{ + // + // open the input stream + // first of all check if opening the file in the constructor was + // successfull + // + fParList = 0; + + if (!OpenStream()) + return kFALSE; + + // + // check if all necessary containers exist in the Parameter list. + // if not create one and add them to the list + // + fRunHeader = (MCorsikaRunHeader*)pList->FindCreateObj("MCorsikaRunHeader"); + if (!fRunHeader) + return kFALSE; + + fEvtHeader = (MCorsikaEvtHeader*)pList->FindCreateObj("MCorsikaEvtHeader"); + if (!fEvtHeader) + return kFALSE; + + fEvtData = (MPhotonData*)pList->FindCreateObj("MPhotonData"); + if (!fEvtData) + return kFALSE; + + fEvent = (MPhotonEvent*)pList->FindCreateObj("MPhotonEvent"); + if (!fEvent) + return kFALSE; + + //*fLog << inf << "Maintaining " << fEvent->GetClassName() << " found in MCorsikaEvent." << endl; + + *fLog << inf << "Calculating number of total events..." << flush; + if (!CalcNumTotalEvents()) + return kFALSE; + *fLog << inf << " " << fNumTotalEvents << " found." << endl; + + fNumFile=0; + fNumEvents=0; + + fParList = pList; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read a single event from the stream +// +Bool_t MCorsikaRead::ReadEvent(istream &fin) +{ + // + // Read in the next EVENT HEADER (see specification), + // if there is no next event anymore stop eventloop + // + Int_t rc = fEvtHeader->ReadEvt(fin, *fRunHeader); //read event header block + if (!rc) + return kFALSE; + + rc = fEvent->ReadCorsikaEvt(fin); + + /* + // Check if we are allowed to stop reading in the middle of the data + if (rc==kFALSE && !fForce) + { + *fLog << err; + *fLog << "ERROR - End of file in the middle of the data stream!" << endl; + *fLog << " Set the force-option to force processing." << endl; + return kERROR; + } + */ + + return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fin) : rc; +} + +// -------------------------------------------------------------------------- +// +// The Process reads one event from the binary file: +// - The event header is read +// - the run header is read +// - all crate information is read +// - the raw data information of one event is read +// +Int_t MCorsikaRead::Process() +{ + while (1) + { + if (fIn) + { + // Read a single event from file + const Bool_t rc = ReadEvent(*fIn); + + // kFALSE means: end of file (try next one) + if (rc!=kFALSE) + return rc; + + if (!fRunHeader->ReadEvtEnd(*fIn)) + if (!fForceMode) + return kERROR; + } + + // + // If an event could not be read from file try to open new file + // + const Int_t rc = OpenNextFile(); + if (rc!=kTRUE) + return rc; + } + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Close the file. Check whether the number of read events differs from +// the number the file should containe (MCorsikaRunHeader). Prints a warning +// if it doesn't match. +// +Int_t MCorsikaRead::PostProcess() +{ + // + // Sanity check for the number of events + // + if (fNumEvents==GetNumExecutions()-1 || GetNumExecutions()==0) + return kTRUE; + + *fLog << warn << dec; + *fLog << "Warning - number of read events (" << GetNumExecutions()-1; + *fLog << ") doesn't match number in run header(s) ("; + *fLog << fNumEvents << ")." << endl; + + return kTRUE; +} Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h (revision 9182) @@ -0,0 +1,68 @@ +#ifndef MARS_MCorsikaRead +#define MARS_MCorsikaRead + +#ifndef MARS_MRead +#include "MRead.h" +#endif + +class TList; +//class MZlib; +class MTaskList; + +class MCorsikaRunHeader; +class MCorsikaEvtHeader; +class MPhotonData; +class MPhotonEvent; + +class MCorsikaRead : public MRead +{ +private: + MCorsikaRunHeader *fRunHeader; // run header information container to fill from file + MCorsikaEvtHeader *fEvtHeader; // event header information container to fill from file + MPhotonEvent *fEvent; // event information + MPhotonData *fEvtData; // raw evt data information container to fill from file + + Bool_t fForceMode; // Force mode skipping defect events + + TList *fFileNames; // list of file names + UInt_t fNumFile; //! number of next file + UInt_t fNumEvents; //! input stream (file to read from) + UInt_t fNumTotalEvents; //! total number of events in all files + + fstream *fIn; //! input stream (file to read from) + + MParList *fParList; //! tasklist to call ReInit from + + //UInt_t fInterleave; + //Bool_t fForce; + + virtual Bool_t OpenStream() { return kTRUE; } + + Bool_t ReadEvtEnd(); + Int_t OpenNextFile(Bool_t print=kTRUE); + Bool_t CalcNumTotalEvents(); + Bool_t ReadEvent(istream &fin); + + Int_t PreProcess(MParList *pList); + Int_t Process(); + Int_t PostProcess(); + +public: + MCorsikaRead(const char *filename=NULL, const char *name=NULL, const char *title=NULL); + ~MCorsikaRead(); + + //static Byte_t IsFileValid(const char *name); + + //void SetInterleave(UInt_t i) { fInterleave = i; } + void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; } + + TString GetFullFileName() const; + + Int_t AddFile(const char *fname, Int_t entries=-1); + Bool_t Rewind(); + UInt_t GetEntries() { return fNumTotalEvents/*/fInterleave*/; } + + ClassDef(MCorsikaRead, 0) // Task to read the raw data binary file +}; + +#endif Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc (revision 9182) @@ -0,0 +1,258 @@ +/* ======================================================================== *\ +! +! * +! * 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 + Qi Zhe, 06/2007 + +! Copyright: MAGIC Software Development, 2000-2004 +! +! +\* ======================================================================== */ + +///////////////////////////////////////////////////////////////////////////// +// +// MCorsikaRunHeader +// +// Root storage container for the RUN HEADER information +// +// +// RAW DATA FORMAT VERSION +// ======================= +// +// 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 +// +// +// MCorsikaRunHeader CLASS VERSION +// =========================== +// +// Format Version 6: +// ----------------- +// - added fNumBytesPerSample; +// - added fFreqSampling; +// - added fNumSignificantBits; +// +// 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 "MCorsikaRunHeader.h" + +#include +#include + +#include "MLog.h" +#include "MLogManip.h" + +ClassImp(MCorsikaRunHeader); + +using namespace std; + +// -------------------------------------------------------------------------- +// +// Default constructor. Creates array which stores the pixel assignment. +// +// +MCorsikaRunHeader::MCorsikaRunHeader(const char *name, const char *title) + : fNumObsLevel(0), + fZdMin(0), fZdMax(-1), fAzMin(0), fAzMax(0), +fViewConeInnerAngle(0), fViewConeOuterAngle(-1) +{ + fName = name ? name : "MCorsikaRunHeader"; + fTitle = title ? title : "Raw Run Header Information"; +} + +// -------------------------------------------------------------------------- +// +// Read in one run header from the binary file +// +Bool_t MCorsikaRunHeader::ReadEvt(istream& fin) +{ + char runh[4]; + fin.read(runh, 4); + if (memcmp(runh, "RUNH", 4)) + { + *fLog << err << "ERROR - Wrong identifier: RUNH expected." << endl; + return kFALSE; + } + + Float_t f[68]; + fin.read((char*)f, 68); + + fRunNumber = TMath::Nint(f[0]); + fNumEvents = 0; + + fRunStart.SetCorsikaTime(f[1]); + + fProgramVersion = f[2]; //FIXME: INT??? + fNumObsLevel = TMath::Nint(f[3]); + + memset(fObsLevel, 0, 10*4); + memcpy(fObsLevel, f+4, fNumObsLevel*4); + + fSlopeSpectrum = f[14]; + fEnergyMin = f[15]; + fEnergyMax = f[16]; + + fin.seekg(1020, ios::cur); // skip the remaining data of this block + + return kTRUE; +} + +Bool_t MCorsikaRunHeader::ReadEvtEnd(istream& fin) +{ + char runh[4]; + fin.read(runh, 4); + if (memcmp(runh, "RUNE", 4)) + { + *fLog << err << "ERROR - Wrong identifier: RUNE expected." << endl; + return kFALSE; + } + + Float_t f[2]; + fin.read((char*)f, 2*4); + + const UInt_t runnum = TMath::Nint(f[0]); + if (runnum!=fRunNumber) + { + *fLog << err << "ERROR - Mismatch in stream: Run number in RUNE ("; + *fLog << runnum << ") doesn't match RUNH (" << fRunNumber << ")." << endl; + return kFALSE; + } + + fNumEvents = TMath::Nint(f[1]); + + fin.seekg(270*4, ios::cur); // skip the remaining data of this block + + return kTRUE; +} + +Bool_t MCorsikaRunHeader::SeekEvtEnd(istream &fin) +{ + // Search subblockwise backward (Block: 5733*4 = 21*273*4) + for (int i=1; i<22; i++) + { + fin.seekg(-i*273*4, ios::end); + + char runh[4]; + fin.read(runh, 4); + + if (!memcmp(runh, "RUNE", 4)) + { + fin.seekg(-4, ios::cur); + return kTRUE; + } + } + + return kFALSE; +} + +// -------------------------------------------------------------------------- +// +// print run header information on *fLog. The option 'header' supresses +// the pixel index translation table. +// +void MCorsikaRunHeader::Print(Option_t *t) const +{ + *fLog << all << endl; + *fLog << "Run Number: " << fRunNumber << " (" << fRunStart.GetStringFmt("%d.%m.%Y") << ", V" << fProgramVersion << ")" << endl; + if (fNumEvents>0) + *fLog << "Num Events: " << fNumEvents << endl; + *fLog << "Obs Level: "; + for (Byte_t i=0; i0) + *fLog << "ViewCone: " << fViewConeInnerAngle << "deg-" << fViewConeOuterAngle << "deg" << endl; + + if (fZdMax>=0) + { + *fLog << "Zd/Az: " << fZdMin << "deg"; + if (fZdMin==fZdMax) + *fLog << " (fixed)"; + else + *fLog << "-" << fZdMax << "deg"; + *fLog << " / " << fAzMin << "deg"; + if (fAzMin==fAzMax) + *fLog << " (fixed)"; + else + *fLog << "-" << fAzMax << "deg"; + *fLog << endl; + } +} + Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h (revision 9182) @@ -0,0 +1,60 @@ +#ifndef MARS_MCorsikaRunHeader +#define MARS_MCorsikaRunHeader +/////////////////////////////////////////////////////////////////////// +// // +// MRunHeader // +// // +/////////////////////////////////////////////////////////////////////// + +#ifndef MARS_MTime +#include "MTime.h" +#endif + +class MCorsikaRunHeader : public MParContainer +{ + friend class MCorsikaEvtHeader; + +private: + UInt_t fRunNumber; // Run number + UInt_t fNumEvents; // Number of events + MTime fRunStart; // Date of begin (yymmdd) + Float_t fProgramVersion; // Version of program + + Byte_t fNumObsLevel; // Number of observation levels + Float_t fObsLevel[10]; //[fNumObsLevel] Observation levels [cm] + + Float_t fSlopeSpectrum; // Slope of energy spectrum + Float_t fEnergyMin; // Lower limit of energy range + Float_t fEnergyMax; // Upper limit of energy range + + Float_t fZdMin; // [rad] Zenith distance + Float_t fZdMax; // [rad] Zenith distance + Float_t fAzMin; // [rad] Azimuth (north=0; west=90) + Float_t fAzMax; // [rad] Azimuth (north=0; west=90) + + Float_t fViewConeInnerAngle; // [deg] + Float_t fViewConeOuterAngle; // [deg] + +public: + MCorsikaRunHeader(const char *name=NULL, const char *title=NULL); + + UInt_t GetFormatVersion() const { return (UInt_t)fProgramVersion; } + Bool_t HasViewCone() const { return fViewConeOuterAngle>0; } + + Float_t GetZdMin() const { return fZdMin; } + Float_t GetZdMax() const { return fZdMax; } + + Float_t GetAzMin() const { return fAzMin; } + Float_t GetAzMax() const { return fAzMax; } + + UInt_t GetNumEvents() const { return fNumEvents; } + + void Print(Option_t *t=NULL) const; + + Bool_t ReadEvt(istream& fin); + Bool_t ReadEvtEnd(istream& fin); + Bool_t SeekEvtEnd(istream &fin); + + ClassDef(MCorsikaRunHeader, 1) // storage container for general info +}; +#endif Index: trunk/MagicSoft/Mars/mcorsika/Makefile =================================================================== --- trunk/MagicSoft/Mars/mcorsika/Makefile (revision 9182) +++ trunk/MagicSoft/Mars/mcorsika/Makefile (revision 9182) @@ -0,0 +1,33 @@ +################################################################## +# +# makefile +# +# for the MARS software +# +################################################################## +include ../Makefile.conf.$(OSTYPE) +include ../Makefile.conf.general + +#------------------------------------------------------------------------------ + +# +# Handling name of the Root Dictionary Files +# +CINT = Corsika + +# +# connect the include files defined in the config.mk file +# +INCLUDES = -I. -I../mbase -I../mfileio -I../mmc -I../msim -I../mgui + +SRCFILES = MCorsikaRunHeader.cc \ + MCorsikaEvtHeader.cc \ + MCorsikaRead.cc + +############################################################ + +all: $(OBJS) + +include ../Makefile.rules + +mrproper: clean rmbak