Index: /trunk/Mars/Changelog =================================================================== --- /trunk/Mars/Changelog (revision 10059) +++ /trunk/Mars/Changelog (revision 10060) @@ -18,4 +18,48 @@ -*-*- END OF LINE -*-*- + + 2010/11/29 Reiner Rohlfs + + * msim/MPhotonData.cc: + - Use all data if telescope array is not defined + + * msim/MPhotonEvent.[h,cc]: + - remove following methods: + Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i); + Int_t ReadCorsikaEvt(istream &fin, Int_t i); + - replace those methods by + Int_t ReadEventIoEvt(MCorsikaFormat *fInFormat); + Int_t ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx); + + * mcorsika/MCorsikaRunHeader.h mcorsika/MCorsikaRunHeader.cc + - Split method ReadEvtEnd() into two functions: + ReadEvtEnd() and ReadEventHeader() + + * mcorsika/MCorsikaEvtHeader.[h,cc]: + - method ReadEvt does not read data from file by itself, but gets + buffer as input + - test that number of reuse of same shour is not greater than 20 + - method ReadEvtEnd does not seek to the EVTE block and does not + read the Bloch header any more. + + * mcorsika/MCorsikaFormat.[h,cc]: + - remove following methods: + SeekNextBlock(), UnreadLastHeader(), UnreadLastData(), StorePos(), + ResetPos(), Rewind(), GetNextEvent(), GetCurrPos(), + NextEventObject() and NextPhotonObject() + - new method: NextBlock() + + * mcorsika/MCorsikaRead.[h,cc]: + - new design of program flow. It is now determined by the order of + the data blocks in the input file. + - Depending on the variables fTelescopeIdx and fArrayIdx, only data + of the requested telescope and array are stored in the output + file. + + * readcorsika.cc: + - two new command line arguments: -A=arrayNo and -T=telescopeNo + The values of these parameters are store in the MCorsikaRead - + class. + 2010/11/22 Thomas Bretz Index: /trunk/Mars/mcorsika/MCorsikaEvtHeader.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 10060) @@ -124,16 +124,8 @@ // -------------------------------------------------------------------------- // -// read the EVTH information from the input stream -// return FALSE if there is no header anymore, else TRUE +// get the EVTH information from the input parameter f // -Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat *fInFormat) +Int_t MCorsikaEvtHeader::ReadEvt(Float_t * f) { - const Int_t rc=fInFormat->SeekNextBlock("EVTH", 1202); - if (rc!=kTRUE) - return rc; - - Float_t f[273]; - if (!fInFormat->ReadData(272, f)) - return kERROR; fEvtNumber = TMath::Nint(f[0]); @@ -160,7 +152,15 @@ // f[96] // Number of reuse of event [max=20] + fTotReuse = f[96]; + if (fTotReuse > 20) + { + *fLog << err << "Number of reuse of shower is " << fTotReuse + << " But maximum implemented: 20" << endl; + return kFALSE; + } memcpy(fTempX, f +97, 20*sizeof(Float_t)); memcpy(fTempY, f+117, 20*sizeof(Float_t)); + // FIXME: Check fNumReuse<20 @@ -170,5 +170,5 @@ fWeightedNumPhotons = 0; - return fInFormat->Eof() ? kERROR : kTRUE; + return kTRUE; } @@ -178,11 +178,6 @@ Int_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) { - if (fInFormat->SeekNextBlock("EVTE", 1209)!=kTRUE) - return kERROR; - - //fin.seekg(-1088,ios::cur); - - Float_t f[2]; - if (!fInFormat->ReadData(2, f)) + Float_t f[272]; + if (!fInFormat->Read(f, 272 * sizeof(Float_t))) return kERROR; Index: /trunk/Mars/mcorsika/MCorsikaEvtHeader.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 10060) @@ -21,4 +21,5 @@ UInt_t fEvtNumber; // Event number UInt_t fNumReuse; // Counter of the reuse of the same shower + UInt_t fTotReuse; //! number of reuse of the same shower // UInt_t fParticleID; // Particle ID (see MMcEvtBasic or CORSIKA manual) Float_t fTotalEnergy; // [GeV] @@ -51,4 +52,5 @@ UInt_t GetEvtNumber() const { return fEvtNumber; } UInt_t GetNumReuse() const { return fNumReuse; } + UInt_t GetTotReuse() const { return fTotReuse; } // UInt_t GetParticleID() const { return fParticleID; } @@ -67,11 +69,16 @@ Double_t GetImpact() const; + void GetArrayOffset(Int_t arrayIdx, Float_t & xArrOff, Float_t & yArrOff) + {xArrOff = fTempY[arrayIdx]; yArrOff=-fTempX[arrayIdx]; } + void SetTelescopeOffset(Int_t arrayIdx, Float_t xTelOff, Float_t yTelOff) + {fNumReuse = arrayIdx; fX = xTelOff; fY = yTelOff;} + void IncNumReuse() { fNumReuse++; } - void ResetNumReuse() { fNumReuse=(UInt_t)-1; } + void ResetNumReuse() { fNumReuse=0; } void InitXY() { fX=fTempY[fNumReuse]; fY=-fTempX[fNumReuse]; } void AddXY(Float_t x, Float_t y) { fX+=x; fY+=y; } - Int_t ReadEvt(MCorsikaFormat *informat); // read in event header block + Int_t ReadEvt(Float_t * f); // read in event header block Int_t ReadEvtEnd(MCorsikaFormat *informat); // read in event end block Index: /trunk/Mars/mcorsika/MCorsikaFormat.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaFormat.cc (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaFormat.cc (revision 10060) @@ -36,5 +36,7 @@ #include "MLogManip.h" + using namespace std; + const unsigned int MCorsikaFormat::kSyncMarker = 0xd41f8a37; @@ -55,26 +57,35 @@ } - char buffer[5]="\0\0\0\0"; + char *buffer = new char[5]; + memset(buffer, 0, 5); fileIn->read(buffer, 4); fileIn->seekg(-4, ios::cur); if (strcmp(buffer, "RUNH") == 0) + { + delete [] buffer; return new MCorsikaFormatRaw(fileIn); + } if (*reinterpret_cast(buffer) == kSyncMarker) + { + delete [] buffer; return new MCorsikaFormatEventIO(fileIn); + } gLog << err << "File " << fileName << " is neither a CORSIKA raw nor EventIO file" << endl; delete fileIn; + delete [] buffer; return NULL; } -void MCorsikaFormat::Read(void *ptr, Int_t i) const -{ - fIn->read((char*)ptr, i); -} - +Bool_t MCorsikaFormat::Read(void *ptr, Int_t i) const +{ + fIn->read((char*)ptr, i); + return !fIn->fail(); + +} // -------------------------------------------------------------------------- // @@ -86,54 +97,4 @@ // -------------------------------------------------------------------------- // -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() { @@ -141,41 +102,43 @@ } -// -------------------------------------------------------------------------- -// -// 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"; + +// -------------------------------------------------------------------------- +// +// After a call to this function, the file pointer is located after the +// header of the block. As the event block has no header it is located +// at the beginning of this block, which is as the beginning of the data +Bool_t MCorsikaFormatRaw::NextBlock(Bool_t subBlock, + Int_t & blockType, + Int_t & blockVersion, + Int_t & blockIdentifier, + Int_t & blockLength) const +{ + char blockHeader[5]; + blockHeader[4] = 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) + if (fIn->eof()) 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); -} - + blockVersion = 0; + blockIdentifier = 0; + blockLength = 272 * 4; + + if (strcmp(blockHeader, "RUNH") == 0) + blockType = 1200; + else if (strcmp(blockHeader, "RUNE") == 0) + blockType = 1210; + else if (strcmp(blockHeader, "EVTH") == 0) + blockType = 1202; + else if (strcmp(blockHeader, "EVTE") == 0) + blockType = 1209; + else // the events, they don't have a specific header + { + blockType = 1105; + fIn->seekg(-4, ios::cur); + blockLength += 4; + } + + return kTRUE; +} // -------------------------------------------------------------------------- // @@ -192,5 +155,5 @@ if (!strcmp(runh, "RUNE")) { - fIn->seekg(-4, ios::cur); +// fIn->seekg(-4, ios::cur); return kTRUE; } @@ -199,104 +162,58 @@ 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); - } - - +Bool_t MCorsikaFormatEventIO::NextBlock(Bool_t subBlock, + Int_t & blockType, + Int_t & blockVersion, + Int_t & blockIdentifier, + Int_t & blockLength) const +{ +// we read - synchronisation markerif not subBlock +// - type / version field +// - identification field +// - length +// - unknown field +// - id (first 4 bytes of data field) + + if (subBlock) + { + int blockHeader[3]; + fIn->read((char*)blockHeader, 3 * sizeof(int)); + + if (fIn->eof()) + return kFALSE; + + + blockType = blockHeader[0] & 0xFFFF; + blockVersion = (blockHeader[0] & 0xFFF00000) >> 20; + blockIdentifier = blockHeader[1]; + blockLength = blockHeader[2] & 0x3FFFFFFF; + } + else + { + int blockHeader[4]; + fIn->read((char*)blockHeader, 4 * sizeof(int)); + + if (fIn->eof()) + return kFALSE; + + + blockType = blockHeader[1] & 0xFFFF; + blockVersion = (blockHeader[1] & 0xFFF00000) >> 20; + blockIdentifier = blockHeader[2]; + blockLength = blockHeader[3] & 0x3FFFFFFF; + + if (blockType == 1200 || blockType == 1210 || + blockType == 1202 || blockType == 1209 ) + { + // read the "old" corsika header plus additional 4 unknown byte + char tmp[8]; + fIn->read(tmp, 8); + blockLength -= 8; + } + + } return kTRUE; } @@ -304,24 +221,8 @@ // -------------------------------------------------------------------------- // -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(); - + + // the RUNE block it at the very end of the file. fIn->seekg(-32, ios::end); @@ -334,194 +235,9 @@ { // this seams to be a RUNE (1210) block - fIn->seekg( -4 * (Int_t)(sizeof(int)), ios::cur); - fRunePos = fIn->tellg(); + fIn->seekg( 8, ios::cur); 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; -} Index: /trunk/Mars/mcorsika/MCorsikaFormat.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaFormat.h (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaFormat.h (revision 10060) @@ -17,7 +17,4 @@ std::istream *fIn; - std::streampos fPrevPos; // file position before previous read - std::streampos fPos; - public: static const unsigned int kSyncMarker; @@ -27,24 +24,16 @@ virtual ~MCorsikaFormat(); - virtual Int_t SeekNextBlock(const char * id, unsigned short type) const = 0; - virtual void UnreadLastHeader() const = 0; + virtual Bool_t NextBlock(Bool_t subBlock, Int_t & blockType, Int_t & blockVersion, + Int_t & blockIdentifier, Int_t & blockLength) const = 0; - virtual Bool_t ReadData(Int_t numValues, Float_t * buffer, - Int_t minSeekValues = 272); - virtual void UnreadLastData() const; + void Seek(std::streampos offset) {fIn->seekg(offset, std::ios::cur);} virtual Bool_t SeekEvtEnd() = 0; - virtual void StorePos(); - virtual void ResetPos() const; - virtual void Rewind() const; - virtual Int_t GetNextEvent(Float_t **buffer, Int_t tel=0) = 0; virtual Bool_t IsEventioFormat() const = 0; virtual Bool_t Eof() const; - std::streampos GetCurrPos() const; - - void Read(void *ptr, Int_t i=4) const; + Bool_t Read(void *ptr, Int_t i) const; static MCorsikaFormat *CorsikaFormatFactory(const char *fileName); @@ -60,10 +49,9 @@ : MCorsikaFormat(in) {} - Int_t SeekNextBlock(const char * id, unsigned short type) const; - void UnreadLastHeader() const; + Bool_t NextBlock(Bool_t subBlock, Int_t & blockType, Int_t & blockVersion, + Int_t & blockIdentifier, Int_t & blockLength) const; Bool_t SeekEvtEnd(); - Int_t GetNextEvent(Float_t **buffer, Int_t); Bool_t IsEventioFormat() const {return kFALSE;} }; @@ -72,22 +60,17 @@ class MCorsikaFormatEventIO : public MCorsikaFormat { -private: - std::streampos fRunePos; // file position of the RUNE block public: MCorsikaFormatEventIO(std::istream *in) - : MCorsikaFormat(in) {fRunePos = std::streampos(0);} + : MCorsikaFormat(in) {} - Int_t SeekNextBlock(const char *id, unsigned short type) const; - void UnreadLastHeader() const; + + Bool_t NextBlock(Bool_t subBlock, Int_t & blockType, Int_t & blockVersion, + Int_t & blockIdentifier, Int_t & blockLength) const; Bool_t SeekEvtEnd(); - Int_t GetNextEvent(Float_t **buffer, Int_t tel); Bool_t IsEventioFormat() const { return kTRUE; } -private: - Int_t NextEventObject(Int_t &length) const; - Int_t NextPhotonObject(Int_t &length) const; }; Index: /trunk/Mars/mcorsika/MCorsikaRead.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 10060) @@ -17,4 +17,5 @@ ! ! Author(s): Thomas Bretz 11/2008 + Reiner Rohlfs 11/2010 ! ! Copyright: Software Development, 2000-2008 @@ -58,4 +59,5 @@ using namespace std; + /* ----------- please don't delete and don't care about (Thomas) ------------ @@ -96,7 +98,7 @@ // MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title) - : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fForceMode(kFALSE), - fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0), - /*fIn(0),*/ fInFormat(0), fParList(0), fNumTelescopes(1) + : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fArrayIdx(-1), + fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0), + fNumTotalEvents(0), fInFormat(0), fParList(0), fNumTelescopes(1), fReadState(0) { fName = name ? name : "MRead"; @@ -117,6 +119,4 @@ { delete fFileNames; -// if (fIn) -// delete fIn; if (fInFormat) delete fInFormat; @@ -161,24 +161,4 @@ } - -Bool_t MCorsikaRead::ReadEvtEnd() -{ - if (!fInFormat->SeekEvtEnd()) - { - *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl; - if (!fForceMode) - return kFALSE; - } - - if (!fRunHeader->ReadEvtEnd(fInFormat)) - { - *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl; - if (!fForceMode) - return kFALSE; - } - - return kTRUE; -} - // -------------------------------------------------------------------------- // @@ -191,8 +171,4 @@ // open the input stream and check if it is really open (file exists?) // -/* if (fIn) - delete fIn; - fIn = NULL; -*/ if (fInFormat) delete fInFormat; @@ -231,83 +207,4 @@ fNumFile++; - - // - // Read RUN HEADER (see specification) from input stream - // - if (!fRunHeader->ReadEvt(fInFormat)) - return kERROR; -// if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader)) -// return kERROR; - - fNumTelescopes = 1; - if (fInFormat->IsEventioFormat()) - { - fInFormat->StorePos(); - fInFormat->Rewind(); - - if (!fInFormat->SeekNextBlock(0, 1201)) - { - *fLog << err << "ERROR - Object 1201 not found in file." << endl; - return kERROR; - } - - fInFormat->Read(&fNumTelescopes); - - if (fTelescopeIdx>=fNumTelescopes) - { - *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx << " exceeds number of telescopes " << fNumTelescopes << " in file." << endl; - return kERROR; - } - - fTelescopeX.Set(fNumTelescopes); - fTelescopeY.Set(fNumTelescopes); - fTelescopeZ.Set(fNumTelescopes); - fTelescopeR.Set(fNumTelescopes); - - fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes); - fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes); - fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes); - fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes); - - *fLog << all; - *fLog << "Number of telescopes: " << fNumTelescopes << " (using "; - if (fTelescopeIdx>=0) - *fLog << "telescope " << fTelescopeIdx; - else - *fLog << "all telescopes"; - *fLog << ")" << endl; - - const Int_t lo = fTelescopeIdx<0 ? 0 : fTelescopeIdx; - const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1; - - for (int i=lo; iResetPos(); - } - - - fInFormat->StorePos(); - if (!ReadEvtEnd()) - return kERROR; - fInFormat->ResetPos(); - - - 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) @@ -368,8 +265,45 @@ break; case kTRUE: - if (!ReadEvtEnd()) + + // read run header(1200), telescope position(1201) and + // first event header(1202) + Bool_t status = kTRUE; + Int_t blockType, blockVersion, blockIdentifier, blockLength; + while (status && + fInFormat->NextBlock(fReadState == 4, blockType, blockVersion, + blockIdentifier, blockLength)) + { + if (blockType == 1200) + status = fRunHeader->ReadEvt(fInFormat); + + else if(blockType == 1201) + status = ReadTelescopePosition(); + + else if (blockType == 1202) + { + Float_t buffer[272]; + status = fInFormat->Read(buffer, 272 * sizeof(Float_t)); + status = fRunHeader->ReadEventHeader(buffer); + break; + } + else + fInFormat->Seek(blockLength); + } + + if (status != kTRUE) + return status; + + if (!fInFormat->SeekEvtEnd()) { - rc = kFALSE; - break; + *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl; + if (!fForceMode) + return fForceMode ? kTRUE : kFALSE; + } + + if (!fRunHeader->ReadEvtEnd(fInFormat, kTRUE)) + { + *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl; + if (!fForceMode) + return kFALSE; } @@ -380,10 +314,51 @@ break; } -/* - if (fIn) - delete fIn; - fIn = NULL; -*/ + return rc; +} + +Int_t MCorsikaRead::ReadTelescopePosition() +{ + if (!fInFormat->Read(&fNumTelescopes, 4)) return kERROR; + + if (fTelescopeIdx>=fNumTelescopes) + { + *fLog << err << "ERROR - Requested telescope index " << fTelescopeIdx << + " exceeds number of telescopes " << fNumTelescopes << " in file." << endl; + return kERROR; + } + + fTelescopeX.Set(fNumTelescopes); + fTelescopeY.Set(fNumTelescopes); + fTelescopeZ.Set(fNumTelescopes); + fTelescopeR.Set(fNumTelescopes); + + if (!fInFormat->Read(fTelescopeX.GetArray(), 4*fNumTelescopes)) return kERROR; + if (!fInFormat->Read(fTelescopeY.GetArray(), 4*fNumTelescopes)) return kERROR; + if (!fInFormat->Read(fTelescopeZ.GetArray(), 4*fNumTelescopes)) return kERROR; + if (!fInFormat->Read(fTelescopeR.GetArray(), 4*fNumTelescopes)) return kERROR; + + *fLog << all; + *fLog << "Number of telescopes: " << fNumTelescopes << " (using "; + if (fTelescopeIdx>=0) + *fLog << "telescope " << fTelescopeIdx; + else + *fLog << "all telescopes"; + *fLog << ")" << endl; + + const Int_t lo = fTelescopeIdx<0 ? 0 : fTelescopeIdx; + const Int_t up = fTelescopeIdx<0 ? fNumTelescopes : fTelescopeIdx+1; + + for (int i=lo; iIsEventioFormat()) - { - if (fNumTelescope==0) - { - const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); - if (rc1!=kTRUE) - return rc1; - - if (fEvtHeader->GetNumReuse()==fRunHeader->GetNumReuse()-1) - fEvtHeader->ResetNumReuse(); - fEvtHeader->IncNumReuse(); - - cout << "===== Array " << fEvtHeader->GetNumReuse() << " =====" << endl; - } - - const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fTelescopeIdx); - if (rc2!=kTRUE) - return rc2; - - fEvtHeader->InitXY(); - - const Double_t x = fTelescopeX[fNumTelescope]; - const Double_t y = fTelescopeY[fNumTelescope]; - - fEvtHeader->AddXY(y, -x); - fEvent->AddXY(fEvtHeader->GetX(), fEvtHeader->GetY()); - fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), fRunHeader->GetWavelengthMax()); - - fNumTelescope++; - fNumTelescope%=fNumTelescopes; - - return kTRUE; - } - else - { - // Store the position to re-read a single event several times - if (fEvtHeader->GetNumReuse()GetNumReuse()-1) - fInFormat->ResetPos(); - else - { - fInFormat->StorePos(); - fEvtHeader->ResetNumReuse(); - } - - // Read the event header - const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); - if (rc1!=kTRUE) - return rc1; - - fEvtHeader->IncNumReuse(); - fEvtHeader->InitXY(); - - // In the case of corsika raw data (MMCS only) we have too loop over one event - // several times to read all impact parameters (reusage of events) - // In case of the eventio format we have to decide, the data of which telescope - // we want to extract - - // Read the photons corresponding to the i-th core location - // EventIO: Number of telescope - // Raw: Number of re-use - const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1); - if (rc2!=kTRUE) - return rc2; - - // read event end - return fEvtHeader->ReadEvtEnd(fInFormat); - } -} // -------------------------------------------------------------------------- @@ -530,30 +428,185 @@ Int_t MCorsikaRead::Process() { - while (1) - { - if (fInFormat) - { - // Read a single event from file - const Int_t rc = ReadEvent(); - - // kFALSE means: end of file (try next one) - if (rc!=kFALSE) - return rc; - - - fInFormat->UnreadLastHeader(); - if (!fRunHeader->ReadEvtEnd(fInFormat)) - 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; + + if (fReadState == 11) + { + // we are currently saving the events of the raw format in the root file + if (fEvtHeader->GetNumReuse() == fRunHeader->GetNumReuse()) + { + // all data are saved + fRawEvemtBuffer.resize(0); + fReadState = 3; + } + else + { + fEvtHeader->InitXY(); + Int_t rc = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0], + fRawEvemtBuffer.size() / 7, + fEvtHeader->GetNumReuse()+1); + fEvtHeader->IncNumReuse(); + return rc; + } + } + + while (1) + { + if (fInFormat) + { + Int_t blockType, blockVersion, blockIdentifier, blockLength; + + while (fInFormat->NextBlock(fReadState == 4, blockType, blockVersion, + blockIdentifier, blockLength)) + { +// gLog << "Next block: type=" << blockType << " version=" << blockVersion; +// gLog << " identifier=" << blockIdentifier << " length=" << blockLength << endl; + + + if (fReadState == 4) + { + fTopBlockLength -= blockLength + 12; + if (fTopBlockLength <= 0) + // all data of a top block are read, go back to normal state + fReadState = 3; + } + + Int_t status = kTRUE; + switch (blockType) + { + case 1200: // the run header + status = fRunHeader->ReadEvt(fInFormat); + fReadState = 1; // RUNH is read + break; + + case 1201: // telescope position + status = ReadTelescopePosition(); + break; + + case 1202: // the event header + Float_t buffer[272]; + if (!fInFormat->Read(buffer, 272 * sizeof(Float_t))) + return kFALSE; + + if (fReadState == 1) // first event after RUN header + { + fRunHeader->ReadEventHeader(buffer); + fRunHeader->Print(); + } + + status = fEvtHeader->ReadEvt(buffer); + if (fArrayIdx >= (Int_t)fEvtHeader->GetTotReuse()) + { + *fLog << err << "ERROR - Requested array index " << fArrayIdx << + " exceeds number of arrays " << fEvtHeader->GetTotReuse() << + " in file." << endl; + return kERROR; + } + + + fReadState = 2; + break; + + case 1204: + if (fArrayIdx < 0 || fArrayIdx == blockIdentifier) + { + fReadState = 4; + fTopBlockLength = blockLength; + } + else + // skip this array of telescopes + fInFormat->Seek(blockLength); + + break; + + case 1205: + { + Int_t telIdx = blockIdentifier % 1000; + if (blockVersion == 0 && + (fTelescopeIdx < 0 || fTelescopeIdx == telIdx) && + blockLength > 12) + { + status = fEvent->ReadEventIoEvt(fInFormat); + + Int_t arrayIdx = blockIdentifier / 1000; + Float_t xArrOff, yArrOff; + fEvtHeader->GetArrayOffset(arrayIdx, xArrOff, yArrOff); + fEvtHeader->SetTelescopeOffset(arrayIdx, + xArrOff + fTelescopeY[telIdx], + yArrOff - fTelescopeX[telIdx] ); + fEvent->AddXY(xArrOff + fTelescopeY[telIdx], + yArrOff - fTelescopeX[telIdx]); + fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), + fRunHeader->GetWavelengthMax()); + + if (status == kTRUE) + // end of reading for one telescope in one array == + // end of this Process - step + return status; + } + else + // skip this telescope + fInFormat->Seek(blockLength); + } + break; + + case 1209: // the event end + status = fEvtHeader->ReadEvtEnd(fInFormat); + if (fReadState == 10) + { + // all particles of this event are read, now save them + fReadState = 11; + fEvtHeader->ResetNumReuse(); + + fEvtHeader->InitXY(); + Int_t rc = fEvent->ReadCorsikaEvt(&fRawEvemtBuffer[0], + fRawEvemtBuffer.size() / 7, + fEvtHeader->GetNumReuse()+1); + fEvtHeader->IncNumReuse(); + return rc; + } + else + fReadState = 3; // event closed, run still open + break; + + case 1210: // the run end + status = fRunHeader->ReadEvtEnd(fInFormat, kTRUE); + fNumEvents += fRunHeader->GetNumEvents(); + fRunHeader->SetReadyToSave(); + fReadState = 5; // back to starting point + return status; + break; + + case 1105: // event block of raw format + if (fReadState == 2 || fReadState == 10) + { + Int_t oldSize = fRawEvemtBuffer.size(); + fRawEvemtBuffer.resize(oldSize + blockLength / sizeof(Float_t)); + status = fInFormat->Read(&fRawEvemtBuffer[oldSize], blockLength); + fReadState = 10; + } + else + fInFormat->Seek(blockLength); + break; + + default: + // unknown block, we skip it + fInFormat->Seek(blockLength); + + } + + if (status != kTRUE) + return status; + + } + + } + + // + // 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; } @@ -566,4 +619,5 @@ Int_t MCorsikaRead::PostProcess() { + const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1); @@ -575,5 +629,5 @@ *fLog << warn << dec; - *fLog << "Warning - number of read events (" << GetNumExecutions()-1; + *fLog << "Warning - number of read events (" << GetNumExecutions()-2; *fLog << ") doesn't match number in run header(s) ("; *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl; @@ -582,20 +636,2 @@ } -// -------------------------------------------------------------------------- -// -// Telescope: 1 -// -// In the case of eventio-format select which telescope to extract from the -// file. -1 will extract all telescopes -// -Int_t MCorsikaRead::ReadEnv(const TEnv &env, TString prefix, Bool_t print) -{ - Bool_t rc = kFALSE; - if (IsEnvDefined(env, prefix, "Telescope", print)) - { - rc = kTRUE; - fTelescopeIdx = GetEnvValue(env, prefix, "Telescope", fTelescopeIdx); - } - - return rc; -} Index: /trunk/Mars/mcorsika/MCorsikaRead.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRead.h (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaRead.h (revision 10060) @@ -1,4 +1,6 @@ #ifndef MARS_MCorsikaRead #define MARS_MCorsikaRead + +#include #ifndef MARS_MRead @@ -27,4 +29,5 @@ Int_t fTelescopeIdx; // Index of telescope to be extracted + Int_t fArrayIdx; // Index of telescope-array to be extracted Bool_t fForceMode; // Force mode skipping defect RUNE @@ -34,5 +37,4 @@ UInt_t fNumTotalEvents; //! total number of events in all files -// ifstream *fIn; //! input stream (file to read from) MCorsikaFormat *fInFormat; //! access to input corsika data @@ -46,13 +48,22 @@ TArrayF fTelescopeR; //! Radii of spheres around tel. (unit: cm) + Int_t fReadState; //! 0: nothing read yet + // 1: RUNH is already read + // 2: EVTH is already read + // 3: EVTE is already read, run still open + // 4: inside of a top level block (1205) + // 5: RUNE is already read + // 10: raw events are currently read + // 11: raw events are currently saved + Int_t fTopBlockLength; // remaining length of the current top-level block 1204 + + std::vector fRawEvemtBuffer; //! buffer of raw event data //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(); - Int_t ReadEvent(); + Int_t ReadTelescopePosition(); // MTask @@ -61,6 +72,4 @@ Int_t PostProcess(); - // MParContainer - Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); public: @@ -69,5 +78,6 @@ void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; } - void SetTelescopeIdx(Int_t idx=-1) { fTelescopeIdx = idx; } + void SetTelescopeIdx(Int_t idx) { fTelescopeIdx = idx; } + void SetArrayIdx(Int_t idx) { fArrayIdx = idx; } TString GetFullFileName() const; Index: /trunk/Mars/mcorsika/MCorsikaRunHeader.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRunHeader.cc (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaRunHeader.cc (revision 10060) @@ -85,9 +85,6 @@ Bool_t MCorsikaRunHeader::ReadEvt(MCorsikaFormat * fInFormat) { - if (fInFormat->SeekNextBlock("RUNH", 1200)!=kTRUE) - return kFALSE; - Float_t f[272]; - if (!fInFormat->ReadData(272, f)) + if (!fInFormat->Read(f, 272 * sizeof(Float_t))) return kFALSE; @@ -137,4 +134,15 @@ memcpy(fAtmosphericCoeffB, f+258, 5*4); memcpy(fAtmosphericCoeffC, f+263, 5*4); + + return kTRUE; +} + +// -------------------------------------------------------------------------- +// +// Read in one event header. It is called for the first event header after +// a run header +// +Bool_t MCorsikaRunHeader::ReadEventHeader(Float_t * g) +{ // -------------------- Read first event header ------------------- @@ -158,20 +166,6 @@ // f[145] Muon multiple scattering flag - if (fInFormat->SeekNextBlock("EVTH", 1202)!=kTRUE) - return kFALSE; - - Float_t g[273]; - if (!fInFormat->ReadData(272, g)) - return kFALSE; - - fInFormat->UnreadLastData(); - fInFormat->UnreadLastHeader(); fNumReuse = TMath::Nint(g[96]); // Number i of uses of each cherenkov event - if (fNumReuse!=1) - { - *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; - *fLog << warn << "WARNING - This error is replaced by a warning." << endl; - } fParticleID = TMath::Nint(g[1]); @@ -219,20 +213,20 @@ } -Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) -{ - if (fInFormat->SeekNextBlock("RUNE", 1210)!=kTRUE) +Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat, Bool_t runNumberVerify) +{ + Float_t f[2]; + if (!fInFormat->Read(f, 2 * sizeof(Float_t))) return kFALSE; - Float_t f[2]; - if (!fInFormat->ReadData(2, f)) - return kFALSE; - - 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; - } + if (runNumberVerify) + { + 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]); @@ -345,4 +339,5 @@ *fLog << endl; -} - + +} + Index: /trunk/Mars/mcorsika/MCorsikaRunHeader.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRunHeader.h (revision 10059) +++ /trunk/Mars/mcorsika/MCorsikaRunHeader.h (revision 10060) @@ -124,5 +124,6 @@ // I/O Bool_t ReadEvt(MCorsikaFormat * fInFormat); - Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); + Bool_t ReadEventHeader(Float_t * g); + Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat, Bool_t runNumberVerify); Bool_t SeekEvtEnd(istream &fin); Index: /trunk/Mars/msim/MPhotonData.cc =================================================================== --- /trunk/Mars/msim/MPhotonData.cc (revision 10059) +++ /trunk/Mars/msim/MPhotonData.cc (revision 10060) @@ -222,4 +222,5 @@ { const UInt_t n = TMath::Nint(f[0]); + if (n==0) // FIXME: Do we need to decode the rest anyway? @@ -227,7 +228,10 @@ // Check reuse - const Int_t reuse = (n/1000)%100; // Force this to be 1! - if (reuse!=i) - return kCONTINUE; + if (i >=0) + { + const Int_t reuse = (n/1000)%100; // Force this to be 1! + if (reuse!=i) + return kCONTINUE; + } // This seems to be special to mmcs Index: /trunk/Mars/msim/MPhotonEvent.cc =================================================================== --- /trunk/Mars/msim/MPhotonEvent.cc (revision 10059) +++ /trunk/Mars/msim/MPhotonEvent.cc (revision 10060) @@ -507,314 +507,58 @@ operator[](i).SimWavelength(wmin, wmax); } - -// -------------------------------------------------------------------------- -// -// Read the Event section from the file -// -Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t id) -{ - 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 - Float_t *buffer = 0; - - if (fInFormat->IsEventioFormat()) - { - while (1) - { - const Int_t rc = fInFormat->GetNextEvent(&buffer, id); - if (rc==kERROR) - return kERROR; - if (rc==kFALSE) - break; - - // Loop over number of photons in bunch - while (Add(n).FillEventIO(buffer)) - n++; - } - } - else - { - while (1) - { - const Int_t rc1 = fInFormat->GetNextEvent(&buffer); - if (rc1==kERROR) - return kERROR; - if (rc1==kFALSE) - break; - - const Int_t rc2 = Add(n).FillCorsika(buffer, id); - switch (rc2) - { - 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++; - } - } - -/* - - while (1) - { - Float_t buffer[273]; - Float_t * ptr = buffer; - - - if (!fInFormat->ReadData(273, buffer)) - return kFALSE; - - if (!memcmp(ptr, "EVTE", 4)) - { - - fInFormat->UnreadLastData(); - break; - } - - 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(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 go back - fin.seekg(-4, ios::cur); - - // Do not modify this. It is optimized for execution - // speed and flexibility! - MPhotonData &ph = Add(n); - // It checks how many entries the lookup table has. If it has enough - // entries and the entry was already allocated, we can re-use it, - // otherwise we have to allocate it. - - // Now we read a single cherenkov bunch. Note that for speed reason we have not - // called the constructor if the event was already constructed (virtual table - // set), consequently we must make sure that ReadCorsikaEvent does reset - // all data mebers no matter whether they are read or not. - const Int_t rc = ph.ReadCorsikaEvt(fin); - - // Evaluate result from reading event - switch (rc) - { - case kCONTINUE: continue; // No data in this bunch... skip it. - case kFALSE: return kFALSE; // End of stream - case kERROR: return kERROR; // Error occured - } - - // FIXME: If fNumPhotons!=1 add the photon more than once - - // Now increase the number of entries which are kept, - // i.e. keep this photon(s) - n++; - } - */ - - Resize(n); - fData.UnSort(); - - SetReadyToSave(); - - //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; - return kTRUE; -} - -// -------------------------------------------------------------------------- -/* -Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i) -{ - Int_t n = 0; - - while (1) - { - // Check the first four bytes - char c[13]; - fin.read(c, 13); - - // End of stream - if (!fin) - return kFALSE; - - // Check if we found the end of the event - if (!memcmp(c, "\nEND---EVENT\n", 13)) - break; - - // The first for byte contained data already --> go back - fin.seekg(-13, ios::cur); - - // Do not modify this. It is optimized for execution - // speed and flexibility! - //TObject *o = n(o)->ReadRflEvt(fin); - const Int_t rc = Add(n).ReadRflEvt(fin, i); - - // Evaluate result from reading event - switch (rc) - { - case kCONTINUE: continue; // No data in this bunch... skip it. - case kFALSE: return kFALSE; // End of stream - case kERROR: return kERROR; // Error occured - } - - // Now increase the number of entries which are kept, - // i.e. keep this photon(s) - n++; - } - - Shrink(n); - - SetReadyToSave(); - - // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; - return kTRUE; -} -*/ +Int_t MPhotonEvent::ReadEventIoEvt(MCorsikaFormat *fInFormat) +{ + Int_t bunchHeader[3]; + fInFormat->Read(bunchHeader, 3 * sizeof(Int_t)); + + Int_t n = 0; + + for (int bunch = 0; bunch < bunchHeader[2]; bunch++) + { + Float_t buffer[8]; + fInFormat->Read(buffer, 8 * sizeof(Float_t)); + + if (Add(n).FillEventIO(buffer)) + n++; + } + + Resize(n); + fData.UnSort(); + + SetReadyToSave(); + + //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; + return kTRUE; + +} + +Int_t MPhotonEvent::ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx) +{ + Int_t n = 0; + + for (Int_t event = 0; event < numEvents; event++) + { + const Int_t rc = Add(n).FillCorsika(data + 7 * event, arrayIdx); + + switch (rc) + { + case kCONTINUE: continue; // No data in this bunch... skip it. + case kERROR: return kERROR; // Error occured + } + + // This is a photon we would like to keep later. + // Increase the counter by one + n++; + } + + Resize(n); + fData.UnSort(); + + SetReadyToSave(); + + //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; + return kTRUE; + +} // -------------------------------------------------------------------------- Index: /trunk/Mars/msim/MPhotonEvent.h =================================================================== --- /trunk/Mars/msim/MPhotonEvent.h (revision 10059) +++ /trunk/Mars/msim/MPhotonEvent.h (revision 10060) @@ -60,7 +60,6 @@ // I/O - Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i); - Int_t ReadCorsikaEvt(istream &fin, Int_t i); - //Int_t ReadRflEvt(istream &fin, Int_t i); + Int_t ReadEventIoEvt(MCorsikaFormat *fInFormat); + Int_t ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx); // TObject Index: /trunk/Mars/readcorsika.cc =================================================================== --- /trunk/Mars/readcorsika.cc (revision 10059) +++ /trunk/Mars/readcorsika.cc (revision 10060) @@ -38,7 +38,9 @@ gLog << all << endl; gLog << "Sorry the usage is:" << endl; - gLog << " readcorsika [-h] [-?] [-vn] [-dec] [-a0] inputfile[.raw] [outputfile[.root]]" << endl << endl; - gLog << " inputfile: corsika raw file" << endl; - gLog << " outputfile: root file" << endl; + gLog << " readcorsika [-h] [-?] [-vn] [-dec] [-a0] [-A=arrayNo] [-T=telescopeNo] inputfile[.raw] [outputfile[.root]]" << endl << endl; + gLog << " inputfile: corsika raw file or eventio file" << endl; + gLog << " outputfile: root file" << endl; + gLog << " -A=arrayNo: use data only of array number arrayNo" << endl; + gLog << " -T=telescopeNo: use data only of telescope number telescopeNo (used only for eventio input file)" << endl; gLog << " -ff Force reading of file even if problems occur" << endl; gLog.Usage(); @@ -71,6 +73,9 @@ const Int_t kCompLvl = arg.GetIntAndRemove("--comp=", 1); + const Int_t kArrayNo = arg.GetIntAndRemove("-A=", -1); + const Int_t kTelNo = arg.GetIntAndRemove("-T=", -1); const Bool_t kForce = arg.HasOnlyAndRemove("-f"); const Bool_t kForceRd = arg.HasOnlyAndRemove("-ff"); + // @@ -149,4 +154,6 @@ MCorsikaRead read(kNamein); read.SetForceMode(kForceRd); + read.SetArrayIdx(kArrayNo); + read.SetTelescopeIdx(kTelNo); tasks.AddToList(&read);