/* ======================================================================== *\ ! ! * ! * 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): Reiner Rohlfs 2010 ! Author(s): Thomas Bretz 2010 ! ! Copyright: Software Development, 2000-2010 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MCorsikaFormat // ////////////////////////////////////////////////////////////////////////////// #include "MCorsikaFormat.h" #include #include #include "MLogManip.h" using namespace std; const unsigned int MCorsikaFormat::kSyncMarker = 0xd41f8a37; // -------------------------------------------------------------------------- // MCorsikaFormat *MCorsikaFormat::CorsikaFormatFactory(const char * fileName) { ifstream * fileIn = new ifstream(fileName); const Bool_t noexist = !(*fileIn); if (noexist) { gLog << err << "Cannot open file " << fileName << ": "; gLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl; delete fileIn; return NULL; } char buffer[5]="\0\0\0\0"; fileIn->read(buffer, 4); fileIn->seekg(-4, ios::cur); if (strcmp(buffer, "RUNH") == 0) return new MCorsikaFormatRaw(fileIn); if (*reinterpret_cast(buffer) == kSyncMarker) return new MCorsikaFormatEventIO(fileIn); gLog << err << "File " << fileName << " is neither a CORSIKA raw nor EventIO file" << endl; delete fileIn; return NULL; } void MCorsikaFormat::Read(void *ptr, Int_t i) const { fIn->read((char*)ptr, i); } // -------------------------------------------------------------------------- // Bool_t MCorsikaFormat::Eof() const { return fIn->eof(); } // -------------------------------------------------------------------------- // streampos MCorsikaFormat::GetCurrPos() const { return fIn->tellg(); } // -------------------------------------------------------------------------- // 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() const { fIn->seekg(fPrevPos, ios::beg); } // -------------------------------------------------------------------------- // void MCorsikaFormat::StorePos() { fPos = fIn->tellg(); } // -------------------------------------------------------------------------- // void MCorsikaFormat::ResetPos() const { fIn->seekg(fPos, ios::beg); } void MCorsikaFormat::Rewind() const { fIn->seekg(0, 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. // Int_t MCorsikaFormatRaw::SeekNextBlock(const char * id, unsigned short type) const { char blockHeader[5]="\0\0\0\0"; fIn->read(blockHeader, 4); if (strcmp(blockHeader, id)==0) return kTRUE; // 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. if (strcmp(blockHeader, "RUNE")==0) return kFALSE; gLog << err << "ERROR - Wrong identifier: " << id << " expected."; gLog << " But read " << blockHeader << " from file." << endl; return kERROR; } // -------------------------------------------------------------------------- // void MCorsikaFormatRaw::UnreadLastHeader() const { 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[5]="\0\0\0\0"; fIn->read(runh, 4); if (!strcmp(runh, "RUNE")) { 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; // Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, Int_t) { static Float_t data[273]; static Float_t * next = data + 273; if (next == data + 273) { // read next block of events if (!ReadData(273, data)) return kERROR; 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. // Int_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) const { cout << "Seek " << type << endl; while (1) { // we read - synchronisation marker // - type / version field // - identification field // - length // - unknown field // - id (first 4 bytes of data field) int blockHeader[4]; fIn->read((char*)blockHeader, 4 * sizeof(int)); if (fIn->eof()) { gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - Unexpected end-of-file seeking " << id << " (" << type << ")." << endl; return kERROR; } const unsigned short objType = blockHeader[1] & 0xFFFF; cout << "objType=" << objType << endl; if (type == objType) { if (!id) break; char c[9]; fIn->read(c, 8); if (memcmp(c+4, id, 4)==0) break; c[8] = 0; gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - '" << c+4 << "' doesn't match expected one '" << id << "' in object " << type << endl; return kFALSE; } // 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 if (type == 1202 && objType == 1210) return kFALSE; // a unknown block, we jump to the next one const int length = blockHeader[3] & 0x3fffffff; fIn->seekg(length, ios::cur); } return kTRUE; } // -------------------------------------------------------------------------- // void MCorsikaFormatEventIO::UnreadLastHeader() const { fIn->seekg( -6 * (Int_t)(sizeof(int)), ios::cur); } // -------------------------------------------------------------------------- // Bool_t MCorsikaFormatEventIO::SeekEvtEnd() { if (fRunePos != streampos(0)) { fIn->seekg(fRunePos, ios::beg); return kTRUE; } // it is the first time we are looking for the RUNE block // is the RUNE block at the very end of the file? const streampos currentPos = fIn->tellg(); fIn->seekg(-32, ios::end); unsigned int blockHeader[4]; fIn->read((char*)blockHeader, 4 * sizeof(int)); if ( blockHeader[0] == kSyncMarker && (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)!=kTRUE) return kFALSE; UnreadLastHeader(); fRunePos = fIn->tellg(); 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; // Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, Int_t telescope) { static Float_t *data = NULL; static int lengthPhotonData = 0; static int lengthEvent = 0; if (lengthPhotonData>0) { cout << "Return Bunch l2=" << lengthPhotonData << endl; lengthPhotonData -= 8 * sizeof(Float_t); *buffer += 8; // Return the pointer to the start of the current photon data return kTRUE; } if (data) { delete [] data; data = NULL; cout << "Return end-of-tel LE=" << lengthEvent << endl; return kFALSE; } if (lengthEvent==0) { cout << "Searching 1204... " << flush; // The length of the block is stored and we use it to determine // when the next top level block is expected const Int_t rc = NextEventObject(lengthEvent); if (rc==kERROR) return kERROR; if (!rc) { cout << "skip EVTE." << endl; return kFALSE; } cout << "found new array." << endl; } cout << "Searching 1205..." << flush; // Look for the next event photon bunch (object type 1205) const Int_t tel = NextPhotonObject(lengthPhotonData); if (tel<0) return kERROR; lengthEvent -= lengthPhotonData+12; // Length of data + Length of header lengthPhotonData -= 12; // Length of data before the photon bunches fIn->seekg(12, ios::cur); // Skip this data cout << " found (tel=" << tel << ",LEN=" << lengthEvent << ",len=" << lengthPhotonData << ")... " << flush; if (lengthPhotonData==0) { cout << "empty!" << endl; return kFALSE; } const UInt_t cnt = lengthPhotonData / sizeof(Float_t); // Read next object (photon bunch) data = new Float_t[cnt]; if (!ReadData(cnt, data, 0)) { delete [] data; data = NULL; return kERROR; } cout << "return!" << endl; lengthPhotonData -= 8 * sizeof(Float_t); *buffer = data; // Return the pointer to the start of the current photon data 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 // Int_t MCorsikaFormatEventIO::NextEventObject(Int_t &length) const { while (1) { // we read - synchronisation marker // - type / version field // - identification field // - length UInt_t blockHeader[4]; fIn->read((char*)blockHeader, 4 * sizeof(Int_t)); if (fIn->eof()) { gLog << err << "MCorsikaFormatEventIO::NextEventObject: ERROR - Unexpected end-of-file." << endl; return kERROR; } // For speed reasons we skip the check of the synchronization marker // Decode the object type const unsigned short objType = blockHeader[1] & 0xFFFF; cout << "objType=" << objType << endl; // Decode length of block length = blockHeader[3] & 0x3fffffff; // Check if we found the next array (reuse / impact parameter) // blockHeader[2] == array number (reuse) if (objType==1204) return kTRUE; // we found an eventEnd block, reset file pointer if (objType==1209) break; // a unknown block, we jump to the next one fIn->seekg(length, ios::cur); } fIn->seekg( -4 * (Int_t)(sizeof(Int_t)), ios::cur); length = 0; return kFALSE; } // -------------------------------------------------------------------------- // Int_t MCorsikaFormatEventIO::NextPhotonObject(Int_t &length) const { UInt_t blockHeader[3]; // we read - synchronisation marker // - type / version field // - identification field // - length fIn->read((char*)blockHeader, 3 * sizeof(UInt_t)); if (fIn->eof()) { gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected end-of-file.." << endl; return -1; } const unsigned short objType = blockHeader[0] & 0xFFFF; if (objType != 1205) { gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object type " << objType << " (expected 1205)." << endl; return -1; } const unsigned short version = (blockHeader[0] >> 20) & 0x0FFF; if (version != 0) { gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object version " << version << " (expected 0 for object type 1205)." << endl; return -1; } length = blockHeader[2] & 0x3fffffff; // blockHeader[1] == 1000 x array number (reuse) + telescope number return blockHeader[1] % 1000; }