#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; }