Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9616) @@ -34,4 +34,5 @@ ///////////////////////////////////////////////////////////////////////////// #include "MCorsikaEvtHeader.h" +#include "MCorsikaFormat.h" #include @@ -121,16 +122,13 @@ // return FALSE if there is no header anymore, else TRUE // -Int_t MCorsikaEvtHeader::ReadEvt(std::istream &fin) +Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat * fInFormat) { - char evth[4]; - fin.read(evth, 4); - if (memcmp(evth, "EVTH", 4)) - { - fin.seekg(-4, ios::cur); + + if (!fInFormat->SeekNextBlock("EVTH", 1202)) return kFALSE; - } Float_t f[273]; - fin.read((char*)&f, 273*4); + if (!fInFormat->ReadData(272, f)) + return kFALSE; fEvtNumber = TMath::Nint(f[0]); @@ -160,6 +158,6 @@ if (n!=1) { - *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; - return kFALSE; + *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; + *fLog << warn << "WARNING - This error is replaced by a warning." << endl; } @@ -167,7 +165,5 @@ fY = -f[97]; //fY = f[117]; - fin.seekg(1088-273*4, ios::cur); - - return !fin.eof(); + return !fInFormat->Eof(); } @@ -175,10 +171,16 @@ // this member function is for reading the event end block -Bool_t MCorsikaEvtHeader::ReadEvtEnd(std::istream &fin) +Bool_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) { + + if (!fInFormat->SeekNextBlock("EVTE", 1209)) + return kFALSE; + //fin.seekg(-1088,ios::cur); Float_t f[2]; - fin.read((char*)&f, 2*4); + + if (!fInFormat->ReadData(2, f)) + return kFALSE; const UInt_t evtnum = TMath::Nint(f[0]); @@ -192,7 +194,5 @@ fWeightedNumPhotons = f[1]; - fin.seekg(1080,ios::cur); - - return !fin.eof(); + return !fInFormat->Eof(); } Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9616) @@ -13,4 +13,6 @@ //class ifstream; #include + +class MCorsikaFormat; class MCorsikaEvtHeader : public MParContainer @@ -60,6 +62,6 @@ Double_t GetImpact() const; - Int_t ReadEvt(istream& fin); // read in event header block - Bool_t ReadEvtEnd(istream& fin); // read in event end block + Int_t ReadEvt(MCorsikaFormat * fInFormat); // read in event header block + Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); // read in event end block ClassDef(MCorsikaEvtHeader, 1) // Parameter Conatiner for raw EVENT HEADER Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaFormat.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaFormat.cc (revision 9616) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaFormat.cc (revision 9616) @@ -0,0 +1,420 @@ +#include +#include + + +#include "MLogManip.h" + +#include "MCorsikaFormat.h" + + +using namespace std; + +MCorsikaFormat * CorsikaFormatFactory(MLog * log, const char * fileName) +{ + ifstream * fileIn = new ifstream(fileName); + + const Bool_t noexist = !(*fileIn); + if (noexist) + { + *log << err << "Cannot open file " << fileName << ": "; + *log << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl; + delete fileIn; + fileIn = NULL; + return NULL; + } + + + if (!fileIn->is_open()) + { + *log << err << "Failed to open file " << fileName << endl; + delete fileIn; + fileIn = NULL; + return NULL; + } + + char buffer[4]; + fileIn->read(buffer, 4); + fileIn->seekg(-4, ios::cur); + + int * syncmarker = reinterpret_cast(buffer); + if (memcmp(buffer, "RUNH", 4) == 0) + return new MCorsikaFormatRaw(log, fileIn); + + else if (*syncmarker == -736130505) + return new MCorsikaFormatEventIO(log, fileIn); + + *log << err << "File " << fileName << + " is neither a corsica file, nor an eventio file" << endl; + delete fileIn; + fileIn = NULL; + + return NULL; +} + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + +Bool_t MCorsikaFormat::ReadData(Int_t numValues, Float_t * buffer, + Int_t minSeekValues) +{ + fPrevPos = fIn->tellg(); + + fIn->read((char*)buffer, numValues * sizeof(Float_t)); + + if (numValues < minSeekValues) + // skip the remaining data of this block + fIn->seekg( (minSeekValues - numValues) * 4, ios::cur); + + return !fIn->fail(); + +} + +void MCorsikaFormat::UnreadLastData() +{ + fIn->seekg(fPrevPos, ios::beg); +} + +void MCorsikaFormat::StorePos() +{ + fPos = fIn->tellg(); +//*fLog << all << "storePos: " << fPos << endl; +} + +void MCorsikaFormat::ResetPos() +{ + fIn->seekg(fPos, ios::beg); +} + + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + + +MCorsikaFormat::~MCorsikaFormat() +{ + delete fIn; +} + + +// -------------------------------------------------------------------------- +// +// Reads the next 4 bytes which should be the id. (for example "RUNH", RUNE". +// "EVTH") +// Returns kTRUE if the functions finds the id +// kFALSE if the functions does not find the "id" as the next 4 +// bytes in the file. +// After this call the file position pointer points just after the 4 bytes +// of the id. +// +Bool_t MCorsikaFormatRaw::SeekNextBlock(const char * id, unsigned short type) +{ + char blockHeader[5]; + fIn->read(blockHeader, 4); + + if (memcmp(blockHeader, id, 4)) + { + blockHeader[4] = 0; + if (strcmp(id, "EVTH") != 0 || strcmp(blockHeader, "RUNE") != 0 ) + { + // at the end of a file we are looking for the next Event header, + // but find the end of a run. This is expected, therefor no error + // message. + *fLog << err << "ERROR - Wrong identifier: " << id << " expected." + << " But read " << blockHeader << " from file." << endl; + } + return kFALSE; + } + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +void MCorsikaFormatRaw::UnreadLastHeader() +{ + fIn->seekg(-4, ios::cur); +} + +Bool_t MCorsikaFormatRaw::SeekEvtEnd() +{ + + // 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 kTRUE; +} + + +// -------------------------------------------------------------------------- +// +// Returns one event (7 Float values) after the other until the EventEnd +// block is found. +// If a read error occurred, the readError is set to kTRUE and the function +// returns kFALSE; +Bool_t MCorsikaFormatRaw::GetNextEvent(Float_t ** buffer, Bool_t & readError) +{ + static Float_t data[273]; + static Float_t * next = data + 273; + + if (next == data + 273) + { + // read next block of events + if (!ReadData(273, data)) + { + readError = kTRUE; + return kFALSE; + } + + if (!memcmp(data, "EVTE", 4)) + { + // we found the end of list of events + UnreadLastData(); + return kFALSE; + } + + next = data; + } + + *buffer = next; + next += 7; + + return kTRUE; + +} + + +/////////////////////////////////////////////////////////////////////////////// +/////////////////////////////////////////////////////////////////////////////// + +// -------------------------------------------------------------------------- +// +// Jumps from one top level object to the next until if finds the object with +// correct type and correct id. The id is identifier of the raw corsika block, +// for example "RUNH", RUNE", "EVTH" +// +// Returns kTRUE if the functions finds the type / id +// kFALSE if the functions does not find the type / id. +// +// After this call the file position pointer points just after the 4 bytes +// of the id. +// +Bool_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) +{ + int blockHeader[7]; + + while (1) + { + // we read - synchronisation marker + // - type / version field + // - identification field + // - length + // - unknown field + // - id (first 4 bytes of data field) + fIn->read((char*)blockHeader, 6 * sizeof(int)); + + if (fIn->eof()) + { + *fLog << err << "ERROR - Missing identifier: " << id << + " type: " << type << endl; + return kFALSE; + } + + unsigned short fileType = blockHeader[1] & 0xFFFF; + if (/*memcmp(blockHeader+5, id, 4) == 0 && */ + type == fileType ) + { +//*fLog << all << "found " << id << " type: " << type << endl; + break; + } + + if (type == 1202 && fileType == 1210) + // we are looking for a event header, but found a runEnd. + // This will happen at the end of a run. We stop here looking for + // the next event header + return kFALSE; + + // a unknown block, we jump to the next one +//*fLog << "unknown: " << id << " type: " << fileType << " sub-blocks: " << (blockHeader[3]>>29); +//*fLog << " length: " << (blockHeader[3] & 0x3fffffff) << endl; + int length = blockHeader[3] & 0x3fffffff; + fIn->seekg(length - 2 * sizeof(int), ios::cur); + } + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +void MCorsikaFormatEventIO::UnreadLastHeader() +{ + fIn->seekg( (int)(-6 * sizeof(int)), ios::cur); +} + +// -------------------------------------------------------------------------- +// +Bool_t MCorsikaFormatEventIO::SeekEvtEnd() +{ + if (SeekNextBlock("RUNE", 1210)) + { + UnreadLastHeader(); + return kTRUE; + } + + return kFALSE; +} + +// -------------------------------------------------------------------------- +// +// Returns one event (7 Float values) after the other until the EventEnd +// block is found. +// If a read error occurred, the readError is set to kTRUE and the function +// returns kFALSE; +Bool_t MCorsikaFormatEventIO::GetNextEvent(Float_t ** buffer, + Bool_t & readError) +{ + + static Float_t * data = NULL; + static Float_t * next; + static Int_t topLevelLength = 0; + static Int_t eventLength = 0; + + + while (eventLength == 0) + { + delete [] data; + data = NULL; + + if (topLevelLength == 0) + { + if (!NextTopLevelBlock(topLevelLength, readError)) + return kFALSE; + } + + if (!NextEventBlock(eventLength, readError)) + return kFALSE; + topLevelLength -= eventLength + 3 * sizeof(int); + + // read next block of events + data = new Float_t [eventLength / sizeof(Float_t)]; + if (!ReadData(eventLength / sizeof(Float_t), data, 0)) + { + delete [] data; + data = NULL; + readError = kTRUE; + return kFALSE; + } + next = data + 3; + eventLength -= 3 * sizeof(Float_t); + } + + eventLength -= 8 * sizeof(Float_t); + *buffer = next; + next += 8; + + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Looks for the next Block with type 1204 and return kTRUE. +// The function also stops moving forward in the file, if it finds a +// EventEnd block (1209). In this case kFALSE is returned +Bool_t MCorsikaFormatEventIO::NextTopLevelBlock(Int_t & length, + Bool_t & readError) +{ + Int_t blockHeader[4]; + + while (1) + { + // we read - synchronisation marker + // - type / version field + // - identification field + // - length + fIn->read((char*)blockHeader, 4 * sizeof(Int_t)); + + if (fIn->eof()) + { + *fLog << err << "ERROR - Missing identifier: 1204 or 1209" << endl; + readError = kTRUE; + return kFALSE; + } + + length = blockHeader[3] & 0x3fffffff; + unsigned short fileType = blockHeader[1] & 0xFFFF; + if (fileType == 1204) +{ +*fLog << all << "found new top level block 1204" << endl; + return kTRUE; +} + if (fileType == 1209) + { + // we found an eventEnd block, reset file pointer + fIn->seekg( (Int_t)(-4 * sizeof(Int_t)), ios::cur); + length = 0; + return kFALSE; + } + + + // a unknown block, we jump to the next one +*fLog << all << "found block " << fileType << endl; + fIn->seekg(length, ios::cur); + } + + return kTRUE; +} + +Bool_t MCorsikaFormatEventIO::NextEventBlock(Int_t & length, + Bool_t & readError) +{ + Int_t blockHeader[3]; + + // we read - synchronisation marker + // - type / version field + // - identification field + // - length + fIn->read((char*)blockHeader, 3 * sizeof(Int_t)); + + if (fIn->eof()) + { + *fLog << err << "ERROR - Missing identifier: 1205" << endl; + readError = kTRUE; + return kFALSE; + } + + unsigned short fileType = blockHeader[0] & 0xFFFF; + if (fileType != 1205) + { + *fLog << err << "ERROR - Unexpected type: " << fileType << "expected 1205" << endl; + readError = kTRUE; + return kFALSE; + } + + unsigned short version = (blockHeader[0] >> 20) & 0x0FFF; + if (version != 0) + { + *fLog << err << "ERROR - Unexpected version: " << version << "expected: 0" << endl; + readError = kTRUE; + return kFALSE; + } + + length = blockHeader[2] & 0x3fffffff; + + return kTRUE; + +} + + Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaFormat.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaFormat.h (revision 9616) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaFormat.h (revision 9616) @@ -0,0 +1,95 @@ +#ifndef MARS_MDataFormat +#define MARS_MDataFormat + + +#ifndef ROOT_Rtypes +#include +#endif + + +#include + +#include "MLog.h" + + +class MCorsikaFormat +{ +protected: + MLog * fLog; + std::istream * fIn; + + std::streampos fPrevPos; // file position before previous read + std::streampos fPos; + +public: + MCorsikaFormat(MLog * log, std::istream * in) + : fIn(in) {fLog = log;} + virtual ~MCorsikaFormat(); + + virtual Bool_t SeekNextBlock(const char * id, unsigned short type) = 0; + virtual void UnreadLastHeader() = 0; + + virtual Bool_t ReadData(Int_t numValues, Float_t * buffer, + Int_t minSeekValues = 272); + virtual void UnreadLastData(); + + virtual Bool_t SeekEvtEnd() = 0; + virtual void StorePos(); + virtual void ResetPos(); + + virtual Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError) = 0; + virtual Bool_t IsEventioFormat() = 0; + + virtual Bool_t Eof() {return fIn->eof();} + + + std::streampos GetCurrPos() {return fIn->tellg();} +}; + + +class MCorsikaFormatRaw : public MCorsikaFormat +{ +private: + +public: + MCorsikaFormatRaw(MLog * log, std::istream * in) + : MCorsikaFormat(log, in) {} + + Bool_t SeekNextBlock(const char * id, unsigned short type); + void UnreadLastHeader(); + + Bool_t SeekEvtEnd(); + + Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError); + Bool_t IsEventioFormat() {return kFALSE;} +}; + + +class MCorsikaFormatEventIO : public MCorsikaFormat +{ +private: + +public: + MCorsikaFormatEventIO(MLog * log, std::istream * in) + : MCorsikaFormat(log, in) {} + + Bool_t SeekNextBlock(const char * id, unsigned short type); + void UnreadLastHeader(); + + Bool_t SeekEvtEnd(); + + Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError); + Bool_t IsEventioFormat() {return kTRUE;} + +private: + Bool_t NextTopLevelBlock(Int_t & length, Bool_t & readError); + Bool_t NextEventBlock(Int_t & length, Bool_t & readError); + +}; + + +MCorsikaFormat * CorsikaFormatFactory(MLog * log, const char * fileName); + + +#endif + Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc (revision 9616) @@ -49,4 +49,5 @@ #include "MStatusDisplay.h" +#include "MCorsikaFormat.h" #include "MCorsikaRunHeader.h" #include "MCorsikaEvtHeader.h" @@ -97,5 +98,5 @@ : fRunHeader(0), fEvtHeader(0), fEvent(0), /*fEvtData(0),*/ fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0), - fIn(0), fParList(0) + fIn(0), fInFormat(0), fParList(0) { fName = name ? name : "MRead"; @@ -118,4 +119,6 @@ if (fIn) delete fIn; + if (fInFormat) + delete fInFormat; } @@ -161,5 +164,5 @@ Bool_t MCorsikaRead::ReadEvtEnd() { - if (!fRunHeader->SeekEvtEnd(*fIn)) + if (!fInFormat->SeekEvtEnd()) { *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl; @@ -168,5 +171,5 @@ } - if (!fRunHeader->ReadEvtEnd(*fIn)) + if (!fRunHeader->ReadEvtEnd(fInFormat)) { *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl; @@ -184,4 +187,5 @@ Int_t MCorsikaRead::OpenNextFile(Bool_t print) { + // // open the input stream and check if it is really open (file exists?) @@ -191,4 +195,8 @@ fIn = NULL; + if (fInFormat) + delete fInFormat; + fInFormat = NULL; + // // Check for the existance of a next file to read @@ -204,32 +212,22 @@ const char *expname = gSystem->ExpandPathName(name); - fIn = new ifstream(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) - { - // Show the number of the last event after - // which we now open a new file - TString txt = GetFileName(); - txt += " @ "; - txt += GetNumExecutions()-1; - fDisplay->SetStatusLine2(txt); - } - } - + fInFormat = CorsikaFormatFactory(fLog, expname); delete [] expname; - if (noexist) + if (fInFormat == NULL) return kERROR; + *fLog << inf << "Open file: '" << name << "'" << endl; + + if (fDisplay) + { + // Show the number of the last event after + // which we now open a new file + TString txt = GetFileName(); + txt += " @ "; + txt += GetNumExecutions()-1; + fDisplay->SetStatusLine2(txt); + } + fNumFile++; @@ -237,13 +235,13 @@ // Read RUN HEADER (see specification) from input stream // - if (!fRunHeader->ReadEvt(*fIn)) + if (!fRunHeader->ReadEvt(fInFormat)) return kERROR; // if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader)) // return kERROR; - const streampos pos = fIn->tellg(); + fInFormat->StorePos(); if (!ReadEvtEnd()) return kERROR; - fIn->seekg(pos, ios::beg); + fInFormat->ResetPos(); @@ -395,5 +393,5 @@ // Read a single event from the stream // -Bool_t MCorsikaRead::ReadEvent(istream &fin) +Bool_t MCorsikaRead::ReadEvent() { // @@ -401,9 +399,9 @@ // if there is no next event anymore stop eventloop // - Int_t rc = fEvtHeader->ReadEvt(fin); //read event header block + Int_t rc = fEvtHeader->ReadEvt(fInFormat); //read event header block if (!rc) return kFALSE; - rc = fEvent->ReadCorsikaEvt(fin); + rc = fEvent->ReadCorsikaEvt(fInFormat); /* @@ -418,5 +416,5 @@ */ - return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fin) : rc; + return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fInFormat) : rc; } @@ -433,8 +431,8 @@ while (1) { - if (fIn) + if (fInFormat) { // Read a single event from file - const Bool_t rc = ReadEvent(*fIn); + const Bool_t rc = ReadEvent(); // kFALSE means: end of file (try next one) @@ -442,5 +440,7 @@ return rc; - if (!fRunHeader->ReadEvtEnd(*fIn)) + + fInFormat->UnreadLastHeader(); + if (!fRunHeader->ReadEvtEnd(fInFormat)) if (!fForceMode) return kERROR; Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h (revision 9616) @@ -13,4 +13,5 @@ class MCorsikaEvtHeader; class MPhotonEvent; +class MCorsikaFormat; class MCorsikaRead : public MRead @@ -29,4 +30,5 @@ ifstream *fIn; //! input stream (file to read from) + MCorsikaFormat * fInFormat; //! access to input corsika data MParList *fParList; //! tasklist to call ReInit from @@ -40,5 +42,5 @@ Int_t OpenNextFile(Bool_t print=kTRUE); Bool_t CalcNumTotalEvents(); - Bool_t ReadEvent(istream &fin); + Bool_t ReadEvent(); Int_t PreProcess(MParList *pList); Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc (revision 9616) @@ -46,4 +46,5 @@ #include "MCorsikaRunHeader.h" +#include "MCorsikaFormat.h" #include @@ -78,16 +79,12 @@ // 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[272*4]; - fin.read((char*)f, 272*4); +Bool_t MCorsikaRunHeader::ReadEvt(MCorsikaFormat * fInFormat) +{ + if (!fInFormat->SeekNextBlock("RUNH", 1200)) + return kFALSE; + + Float_t f[272]; + if (!fInFormat->ReadData(272, f)) + return kFALSE; fRunNumber = TMath::Nint(f[0]); @@ -157,24 +154,19 @@ // f[145] Muon multiple scattering flag - char evth[4]; - fin.read(evth, 4); - if (memcmp(evth, "EVTH", 4)) - { - *fLog << err << "ERROR - Wrong identifier: EVTH expected." << endl; - return kFALSE; - } + if (!fInFormat->SeekNextBlock("EVTH", 1202)) + return kFALSE; Float_t g[273]; - fin.read((char*)&g, 273*4); - if (fin.eof()) - return kFALSE; - - fin.seekg(-274*4, ios::cur); + if (!fInFormat->ReadData(272, g)) + return kFALSE; + + fInFormat->UnreadLastData(); + fInFormat->UnreadLastHeader(); const Int_t n = TMath::Nint(g[96]); // Number i of uses of each cherenkov event if (n!=1) { - *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; - return kFALSE; + *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; + *fLog << warn << "WARNING - This error is replaced by a warning." << endl; } @@ -223,16 +215,13 @@ } -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; - } +Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) +{ + + if (!fInFormat->SeekNextBlock("RUNE", 1210)) + return kFALSE; Float_t f[2]; - fin.read((char*)f, 2*4); + if (!fInFormat->ReadData(2, f)) + return kFALSE; const UInt_t runnum = TMath::Nint(f[0]); @@ -245,6 +234,4 @@ fNumEvents = TMath::Nint(f[1]); - - fin.seekg(270*4, ios::cur); // skip the remaining data of this block return kTRUE; Index: trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h =================================================================== --- trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h (revision 9616) @@ -10,4 +10,6 @@ #include "MTime.h" #endif + +class MCorsikaFormat; class MCorsikaRunHeader : public MParContainer @@ -119,6 +121,6 @@ // I/O - Bool_t ReadEvt(istream& fin); - Bool_t ReadEvtEnd(istream& fin); + Bool_t ReadEvt(MCorsikaFormat * fInFormat); + Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); Bool_t SeekEvtEnd(istream &fin); Index: trunk/MagicSoft/Mars/mcorsika/Makefile =================================================================== --- trunk/MagicSoft/Mars/mcorsika/Makefile (revision 9615) +++ trunk/MagicSoft/Mars/mcorsika/Makefile (revision 9616) @@ -22,4 +22,5 @@ SRCFILES = MCorsikaRunHeader.cc \ + MCorsikaFormat.cc \ MCorsikaEvtHeader.cc \ MCorsikaRead.cc Index: trunk/MagicSoft/Mars/msim/MPhotonData.cc =================================================================== --- trunk/MagicSoft/Mars/msim/MPhotonData.cc (revision 9615) +++ trunk/MagicSoft/Mars/msim/MPhotonData.cc (revision 9616) @@ -245,4 +245,33 @@ // -------------------------------------------------------------------------- // +// Set the data member according to the 8 floats read from a eventio-file. +// This function MUST reset all data-members, no matter whether these are +// contained in the input stream. +// +// Currently we exchange x and y and set y=-y to convert Corsikas coordinate +// system into our own. +// +Int_t MPhotonData::FillEventIO(Float_t f[8]) +{ + fPosX = f[1]; // xpos relative to telescope [cm] + fPosY = -f[0]; // ypos relative to telescope [cm] + fCosU = f[3]; // cos to x + fCosV = -f[2]; // cos to y + fTime = f[4]; // a relative arival time [ns] + fProductionHeight = f[5]; // altitude of emission [cm] + fNumPhotons = f[6]; // photons in this bunch + fWavelength = f[7]; // so far always zeor = unspec. [nm] + + + // Now reset all data members which are not in the stream + fPrimary = MMcEvtBasic::kUNDEFINED; + fTag = -1; + fWeight = 1; + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// // Read seven floats from the stream and call FillCorsika for them. // Index: trunk/MagicSoft/Mars/msim/MPhotonData.h =================================================================== --- trunk/MagicSoft/Mars/msim/MPhotonData.h (revision 9615) +++ trunk/MagicSoft/Mars/msim/MPhotonData.h (revision 9616) @@ -133,4 +133,5 @@ Int_t FillCorsika(Float_t f[7]); + Int_t FillEventIO(Float_t f[7]); Int_t FillRfl(Float_t f[8]); Index: trunk/MagicSoft/Mars/msim/MPhotonEvent.cc =================================================================== --- trunk/MagicSoft/Mars/msim/MPhotonEvent.cc (revision 9615) +++ trunk/MagicSoft/Mars/msim/MPhotonEvent.cc (revision 9616) @@ -121,4 +121,5 @@ #include "MPhotonData.h" +#include "MCorsikaFormat.h" ClassImp(MPhotonEvent); @@ -452,5 +453,5 @@ // Read the Event section from the file // -Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin) +Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat * fInFormat) { Int_t n = 0; @@ -482,24 +483,63 @@ // 1.06GB/ 3s CPU // 1.06GB/22s REAL + Bool_t readError = kFALSE; + Float_t * buffer; + + if ( fInFormat->IsEventioFormat() ) + { + while (fInFormat->GetNextEvent(&buffer, readError)) + { + + const Int_t rc = Add(n).FillEventIO(buffer); + switch (rc) + { + case kCONTINUE: continue; // No data in this bunch... skip it. + case kERROR: return kERROR; // Error occured + //case kFALSE: return kFALSE; // End of stream + } + + // This is a photon we would like to keep later. + // Increase the counter by one + n++; + } + } + else + { + while (fInFormat->GetNextEvent(&buffer, readError)) + { + + const Int_t rc = Add(n).FillCorsika(buffer); + switch (rc) + { + case kCONTINUE: continue; // No data in this bunch... skip it. + case kERROR: return kERROR; // Error occured + //case kFALSE: return kFALSE; // End of stream + } + + // This is a photon we would like to keep later. + // Increase the counter by one + n++; + } + } + if (readError) return kFALSE; + +/* while (1) { - // Stprage for one block - char c[273*4]; - - // Read the first four byte to check whether the next block - // doen't belong to the event anymore - fin.read(c, 4); - if (!fin) + Float_t buffer[273]; + Float_t * ptr = buffer; + + + if (!fInFormat->ReadData(273, buffer)) return kFALSE; - // Check if the event is finished - if (!memcmp(c, "EVTE", 4)) + if (!memcmp(ptr, "EVTE", 4)) + { + + fInFormat->UnreadLastData(); break; - - // Now read the rest of the data - fin.read(c+4, 272*4); - - Float_t *ptr = reinterpret_cast(c); + } + Float_t *end = ptr + 273; @@ -525,4 +565,89 @@ } } + +*/ + Resize(n); + fData.UnSort(); + + SetReadyToSave(); + + //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; + return kTRUE; + +} + +Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin) +{ + Int_t n = 0; + + // --- old I/O --- + // Read only + Reflector (no absorption) + // Muons: 1.06GB/115s = 9.2MB/s (100kEvs) + // Gammas: 1.57GB/275s = 5.7MB/s ( 1kEvs) + + // Read only: + // Gammas: 1.57GB/158s = 9.9MB/s ( 1kEvs) + // Muons: 1.06GB/ 77s = 13.8MB/s (100kEvs) + + // --- new I/O --- + // Read only (don't allocate storage space): + // Gammas: 1.57GB/143s = 11.0MB/s ( 1kEvs) + // Muons: 1.06GB/ 77s = 13.8MB/s (100kEvs) + + // Read only in blocks (with storage allocation): + // Gammas: 1.57GB/28s = 56MB/s ( 1kEvs) + // Muons: 1.06GB/5.2s = 204MB/s (100kEvs) + + // Read only in blocks (without storage allocation): + // similar to just copy + + // Copy with cp + // 1.57GB/ 5s CPU + // 1.57GB/28s REAL + // 1.06GB/ 3s CPU + // 1.06GB/22s REAL + + while (1) + { + // Stprage for one block + char c[273*4]; + + // Read the first four byte to check whether the next block + // doen't belong to the event anymore + fin.read(c, 4); + if (!fin) + return kFALSE; + + // Check if the event is finished + if (!memcmp(c, "EVTE", 4)) + break; + + // Now read the rest of the data + fin.read(c+4, 272*4); + + Float_t *ptr = reinterpret_cast(c); + Float_t *end = ptr + 273; + + // Loop over all sub-blocks (photons) in the block and if + // they contain valid data add them to the array + while (ptr