#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 (type == fileType ) 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) << " pos: " << fIn->tellg() << endl; int length = blockHeader[3] & 0x3fffffff; fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur); } return kTRUE; } // -------------------------------------------------------------------------- // void MCorsikaFormatEventIO::UnreadLastHeader() { fIn->seekg( -6 * (Int_t)(sizeof(int)), ios::cur); } // -------------------------------------------------------------------------- // Bool_t MCorsikaFormatEventIO::SeekEvtEnd() { if (fRunePos != streampos(0)) { fIn->seekg(fRunePos, ios::beg); return kTRUE; } else { // it is the first time we are looking for the RUNE block // is the RUNE block at the very end of the file? std::streampos currentPos = fIn->tellg(); fIn->seekg(-32, ios::end); unsigned int blockHeader[4]; fIn->read((char*)blockHeader, 4 * sizeof(int)); if ( blockHeader[0] == 0xd41f8a37 && (blockHeader[1] & 0xffff) == 1210 && (blockHeader[3] & 0x3fffffff) == 16) { // this seams to be a RUNE (1210) block fIn->seekg( -4 * (Int_t)(sizeof(int)), ios::cur); fRunePos = fIn->tellg(); return kTRUE; } // we do not find a RUNE block at the end of the file // we have to search in the file fIn->seekg(currentPos, ios::beg); if (SeekNextBlock("RUNE", 1210)) { UnreadLastHeader(); fRunePos = fIn->tellg(); 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) return kTRUE; if (fileType == 1209) { // we found an eventEnd block, reset file pointer fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur); length = 0; return kFALSE; } // a unknown block, we jump to the next one 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; }