Index: /trunk/Mars/Changelog =================================================================== --- /trunk/Mars/Changelog (revision 9940) +++ /trunk/Mars/Changelog (revision 9941) @@ -18,4 +18,29 @@ -*-*- END OF LINE -*-*- + + 2010/09/24 Thomas Bretz + + * mcorsika/MCorsikaEvtHeader.cc: + - store impact parameters persistent for further use (because in + eventio format the header is not repeated this is needed) + + * mcorsika/MCorsikaFormat.[h,cc]: + - adapted to be able to read all telescopes and all reused shower + events in eventio + - a lot of cosmetics + - many more comments + - PRELIMINARY + + * mcorsika/MCorsikaRead.[h,cc]: + - implemented variable for selection of the telescope in eventio format + - implemented eventio format (PRELIMINARY) + + * msim/MPhotonData.[h,cc], msim/MPhotonEvent.[h,cc]: + - implemented function to calculate mean event time + - implemented function to simulate the wavelength (lambda^-2) + - implemented function to shift all photons by a certain xy + - adapted ReadCorsikaEvt to changes in MCorsikaFormat + + 2010/09/24 Daniela Dorner Index: /trunk/Mars/mcorsika/MCorsikaEvtHeader.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9940) +++ /trunk/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9941) @@ -160,7 +160,10 @@ // f[96] // Number of reuse of event [max=20] + memcpy(fTempX, f +97, 20*sizeof(Float_t)); + memcpy(fTempY, f+117, 20*sizeof(Float_t)); + // FIXME: Check fNumReuse<20 - fX = f[117 + fNumReuse]; - fY = -f[ 97 + fNumReuse]; +// fX = f[117 + fNumReuse]; +// fY = -f[ 97 + fNumReuse]; fWeightedNumPhotons = 0; Index: /trunk/Mars/mcorsika/MCorsikaEvtHeader.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9940) +++ /trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9941) @@ -40,4 +40,7 @@ Float_t fWeightedNumPhotons; // weighted number of photons arriving at observation level + Float_t fTempX[20]; //! Temporary storage for impact parameter + Float_t fTempY[20]; //! Temporary storage for impact parameter + public: MCorsikaEvtHeader(const char *name=NULL, const char *title=NULL); @@ -67,4 +70,6 @@ void ResetNumReuse() { fNumReuse=(UInt_t)-1; } + void InitXY() { fX=fTempY[fNumReuse]; fY=-fTempX[fNumReuse]; } + Int_t ReadEvt(MCorsikaFormat *informat); // 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 9940) +++ /trunk/Mars/mcorsika/MCorsikaFormat.cc (revision 9941) @@ -198,5 +198,5 @@ // returns kFALSE; // -Bool_t MCorsikaFormatRaw::GetNextEvent(Float_t ** buffer, Bool_t & readError) +Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, UInt_t) { static Float_t data[273]; @@ -204,21 +204,18 @@ if (next == data + 273) - { + { // read next block of events if (!ReadData(273, data)) - { - readError = kTRUE; - return kFALSE; - } + return kERROR; if (!memcmp(data, "EVTE", 4)) - { + { // we found the end of list of events UnreadLastData(); return kFALSE; - } + } next = data; - } + } *buffer = next; @@ -226,5 +223,4 @@ return kTRUE; - } @@ -259,26 +255,23 @@ if (fIn->eof()) - { - gLog << err << "ERROR - Missing identifier: " << id << - " type: " << type << endl; + { + gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - unexpected end-of-file." << endl; return kFALSE; - } - - unsigned short fileType = blockHeader[1] & 0xFFFF; - if (type == fileType) + } + + const unsigned short objType = blockHeader[1] & 0xFFFF; + if (type == objType) 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 + // 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 -//gLog << "unknown: " << id << " type: " << fileType << " sub-blocks: " << (blockHeader[3]>>29); -//gLog << " length: " << (blockHeader[3] & 0x3fffffff) << " pos: " << fIn->tellg() << endl; const int length = blockHeader[3] & 0x3fffffff; - - fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur); + + fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur); } @@ -307,5 +300,5 @@ // is the RUNE block at the very end of the file? - std::streampos currentPos = fIn->tellg(); + const streampos currentPos = fIn->tellg(); fIn->seekg(-32, ios::end); @@ -332,4 +325,5 @@ UnreadLastHeader(); fRunePos = fIn->tellg(); + return kTRUE; } @@ -342,46 +336,75 @@ // 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) - { +Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, UInt_t telescope) +{ + static Float_t *data = NULL; + static Int_t lengthEvent = -1; + static Int_t lengthPhotonData = 0; + + while (lengthPhotonData == 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)) - { + // If we are at the end of the event signal this and + // start next time with a new event + if (lengthEvent==0) + { + lengthEvent=-1; + return kFALSE; // Signal to start with a new event (process task list first) + } + + // Look for the next top level block (an "event", object type 1204) + // A top level block contains the data of a full array simulation + // (a single re-usage of the shower) + if (lengthEvent<0) + { + // 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!=kTRUE) + return rc; + } + + // Look for the next event photon bunch (object type 1205) + const Int_t tel = NextPhotonObject(lengthPhotonData); + if (tel<0) + return kERROR; + + // Decrease the distance to the next "event" by the current photon bunches + lengthEvent -= lengthPhotonData + 3 * sizeof(int); + + // If the telescope is not the one we want to read skip the photon data + if (UInt_t(tel)!=telescope) + { + fPrevPos = fIn->tellg(); + fIn->seekg(lengthPhotonData, ios::cur); + lengthPhotonData=0; + continue; + } + + 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; - readError = kTRUE; - return kFALSE; - } - next = data + 3; - eventLength -= 3 * sizeof(Float_t); - } - - eventLength -= 8 * sizeof(Float_t); - *buffer = next; - next += 8; - + return kERROR; + } + + // The object containing the photon bunches starts with: + // Array[short], Telescope[Short], NumPhotons[Real], NumBunches[Long] + // eventLength contains now the number of bytes with photon data not yet evaluated + lengthPhotonData -= 3 * sizeof(Float_t); + } + + // The photon data itself is 8 floats long. When we return we expect this to be + // evaluated already when the function is called the next time. + lengthPhotonData -= 8 * sizeof(Float_t); + + // Return the pointer to the start of the current photon data + *buffer = data + 3; return kTRUE; @@ -394,9 +417,6 @@ // EventEnd block (1209). In this case kFALSE is returned // -Bool_t MCorsikaFormatEventIO::NextTopLevelBlock(Int_t & length, - Bool_t & readError) const -{ - Int_t blockHeader[4]; - +Int_t MCorsikaFormatEventIO::NextEventObject(Int_t &length) const +{ while (1) { @@ -405,40 +425,46 @@ // - identification field // - length + + UInt_t blockHeader[4]; fIn->read((char*)blockHeader, 4 * sizeof(Int_t)); if (fIn->eof()) - { - gLog << err << "ERROR - Missing identifier: 1204 or 1209" << endl; - readError = kTRUE; - return kFALSE; - } - + { + 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; + + // Decode length of block length = blockHeader[3] & 0x3fffffff; - const unsigned short fileType = blockHeader[1] & 0xFFFF; - if (fileType == 1204) + + // Check if we found the next array (reuse / impact parameter) + // blockHeader[2] == array number (reuse) + if (objType==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; - } - + // 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); - } - - return kTRUE; -} - -// -------------------------------------------------------------------------- -// -Bool_t MCorsikaFormatEventIO::NextEventBlock(Int_t & length, - Bool_t & readError) const -{ - Int_t blockHeader[3]; + 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 @@ -446,32 +472,29 @@ // - identification field // - length - fIn->read((char*)blockHeader, 3 * sizeof(Int_t)); + fIn->read((char*)blockHeader, 3 * sizeof(UInt_t)); if (fIn->eof()) - { - gLog << err << "ERROR - Missing identifier: 1205" << endl; - readError = kTRUE; - return kFALSE; - } - - const unsigned short fileType = blockHeader[0] & 0xFFFF; - if (fileType != 1205) - { - gLog << err << "ERROR - Unexpected type: " << fileType << "expected 1205" << endl; - readError = kTRUE; - return kFALSE; - } + { + 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 << "ERROR - Unexpected version: " << version << "expected: 0" << endl; - readError = kTRUE; - return kFALSE; - } + { + gLog << err << "MCorsikaFormatEventIO::NextPhotonObject: ERROR - Unexpected EventIO object version " << version << " (expected 0 for object type 1205)." << endl; + return -1; + } length = blockHeader[2] & 0x3fffffff; - return kTRUE; - -} + // 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 9940) +++ /trunk/Mars/mcorsika/MCorsikaFormat.h (revision 9941) @@ -2,4 +2,7 @@ #define MARS_MDataFormat +#ifndef MARS_MAGIC +#include "MAGIC.h" +#endif #ifndef ROOT_Rtypes @@ -35,5 +38,5 @@ virtual void ResetPos() const; - virtual Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError) = 0; + virtual Int_t GetNextEvent(Float_t **buffer, UInt_t tel=0) = 0; virtual Bool_t IsEventioFormat() const = 0; @@ -59,5 +62,5 @@ Bool_t SeekEvtEnd(); - Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError); + Int_t GetNextEvent(Float_t **buffer, UInt_t); Bool_t IsEventioFormat() const {return kFALSE;} }; @@ -67,22 +70,21 @@ { private: - std::streampos fRunePos; // file position of the RUNE block + std::streampos fRunePos; // file position of the RUNE block public: - MCorsikaFormatEventIO(std::istream * in) + MCorsikaFormatEventIO(std::istream *in) : MCorsikaFormat(in) {fRunePos = std::streampos(0);} - Bool_t SeekNextBlock(const char * id, unsigned short type) const; - void UnreadLastHeader() const; + Bool_t SeekNextBlock(const char *id, unsigned short type) const; + void UnreadLastHeader() const; - Bool_t SeekEvtEnd(); + Bool_t SeekEvtEnd(); - Bool_t GetNextEvent(Float_t ** buffer, Bool_t & readError); - Bool_t IsEventioFormat() const {return kTRUE;} + Int_t GetNextEvent(Float_t **buffer, UInt_t tel); + Bool_t IsEventioFormat() const { return kTRUE; } private: - Bool_t NextTopLevelBlock(Int_t & length, Bool_t & readError) const; - Bool_t NextEventBlock(Int_t & length, Bool_t & readError) const; - + 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 9940) +++ /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 9941) @@ -96,5 +96,5 @@ // MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title) - : fRunHeader(0), fEvtHeader(0), fEvent(0), /*fEvtData(0),*/ fForceMode(kFALSE), + : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(0), fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0), fIn(0), fInFormat(0), fParList(0) @@ -212,5 +212,5 @@ const char *expname = gSystem->ExpandPathName(name); - fInFormat = CorsikaFormatFactory(fLog, expname); + fInFormat = MCorsikaFormat::CorsikaFormatFactory(expname); delete [] expname; @@ -401,22 +401,54 @@ fInFormat->StorePos(); fEvtHeader->ResetNumReuse(); + + if (fInFormat->IsEventioFormat()) + { + // Read the event header + const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); + if (rc1!=kTRUE) + return rc1; + } + } else - fInFormat->ResetPos(); + { + if (!fInFormat->IsEventioFormat()) + fInFormat->ResetPos(); + } + + if (!fInFormat->IsEventioFormat()) + { + // Read the event header + const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); + if (rc1!=kTRUE) + return rc1; + + } fEvtHeader->IncNumReuse(); - - // Read the event header - const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); - if (rc1!=kTRUE) - return rc1; - - // Check if reusage number should be reset and increase it - // Note, that the trick here is that after reset it is set to -1 + 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 + const Int_t id = fInFormat->IsEventioFormat() ? fTelescopeIdx : fEvtHeader->GetNumReuse()+1; // Read the photons corresponding to the i-th core location - const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1); + // EventIO: Number of telescope + // Raw: Number of re-use + const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, id); if (rc2!=kTRUE) return rc2; + + //---- + if (fInFormat->IsEventioFormat()) + { + fEvent->AddXY(fEvtHeader->GetX(), fEvtHeader->GetY()); + fEvent->SimWavelength(fRunHeader->GetWavelengthMin(), fRunHeader->GetWavelengthMax()); + return kTRUE; + } + //---- + // read event end Index: /trunk/Mars/mcorsika/MCorsikaRead.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRead.h (revision 9940) +++ /trunk/Mars/mcorsika/MCorsikaRead.h (revision 9941) @@ -22,4 +22,5 @@ MPhotonEvent *fEvent; //! event information + UInt_t fTelescopeIdx; // Index of telescope to be extracted Bool_t fForceMode; // Force mode skipping defect RUNE @@ -52,8 +53,6 @@ ~MCorsikaRead(); - //static Byte_t IsFileValid(const char *name); - - //void SetInterleave(UInt_t i) { fInterleave = i; } void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; } + void SetTelescopeIdx(UInt_t idx=0) { fTelescopeIdx = idx; } TString GetFullFileName() const; Index: /trunk/Mars/msim/MPhotonData.cc =================================================================== --- /trunk/Mars/msim/MPhotonData.cc (revision 9940) +++ /trunk/Mars/msim/MPhotonData.cc (revision 9941) @@ -47,4 +47,5 @@ #include +#include #include "MLog.h" @@ -171,4 +172,17 @@ // -------------------------------------------------------------------------- // +// Set the wavelength to a random lambda^-2 distributed value +// between wmin and wmax. +// +void MPhotonData::SimWavelength(Float_t wmin, Float_t wmax) +{ + const Double_t w = gRandom->Uniform(wmin, wmax); + + fWavelength = TMath::Nint(wmin*wmax / w); +} + + +// -------------------------------------------------------------------------- +// // Set the data member according to the 8 floats read from a reflector-file. // This function MUST reset all data-members, no matter whether these are @@ -255,9 +269,10 @@ Int_t MPhotonData::FillEventIO(Float_t f[8]) { - if (TMath::Nint(f[6])!=1) - { - gLog << err << "ERROR - Bunch sizes != 1 are not supported." << endl; - return kFALSE; - } + // photons in this bunch + const UInt_t n = TMath::Nint(f[6]); + if (n==0) + return 0; + + f[6] = n-1; fPosX = f[1]; // xpos relative to telescope [cm] @@ -267,5 +282,4 @@ fTime = f[4]; // a relative arival time [ns] fProductionHeight = f[5]; // altitude of emission [cm] -// fNumPhotons = TMath::Nint(f[6]); // photons in this bunch fWavelength = TMath::Nint(f[7]); // so far always zeor = unspec. [nm] @@ -275,6 +289,7 @@ fWeight = 1; - return kTRUE; -} + return n; +} + /* // -------------------------------------------------------------------------- Index: /trunk/Mars/msim/MPhotonData.h =================================================================== --- /trunk/Mars/msim/MPhotonData.h (revision 9940) +++ /trunk/Mars/msim/MPhotonData.h (revision 9941) @@ -101,4 +101,6 @@ void SetTime(Double_t t) { fTime = t; } + void SimWavelength(Float_t wmin, Float_t wmax); + void Copy(TObject &obj) const; Index: /trunk/Mars/msim/MPhotonEvent.h =================================================================== --- /trunk/Mars/msim/MPhotonEvent.h (revision 9940) +++ /trunk/Mars/msim/MPhotonEvent.h (revision 9941) @@ -39,4 +39,5 @@ Double_t GetMeanX() const; Double_t GetMeanY() const; + Double_t GetMeanT() const; TClonesArray &GetArray() { return fData; } @@ -54,4 +55,7 @@ Int_t Shrink(Int_t n); void Resize(Int_t n); + + void AddXY(Double_t x, Double_t y); + void SimWavelength(Float_t wmin, Float_t wmax); // I/O