Index: /trunk/Mars/mcorsika/MCorsikaEvtHeader.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9942) +++ /trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9943) @@ -71,4 +71,5 @@ 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 Index: /trunk/Mars/mcorsika/MCorsikaFormat.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaFormat.cc (revision 9942) +++ /trunk/Mars/mcorsika/MCorsikaFormat.cc (revision 9943) @@ -72,4 +72,9 @@ } +void MCorsikaFormat::Read(void *ptr, Int_t i) const +{ + fIn->read((char*)ptr, i); +} + // -------------------------------------------------------------------------- // @@ -122,4 +127,9 @@ { fIn->seekg(fPos, ios::beg); +} + +void MCorsikaFormat::Rewind() const +{ + fIn->seekg(0, ios::beg); } @@ -198,5 +208,5 @@ // returns kFALSE; // -Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, UInt_t) +Int_t MCorsikaFormatRaw::GetNextEvent(Float_t **buffer, Int_t) { static Float_t data[273]; @@ -242,6 +252,5 @@ Bool_t MCorsikaFormatEventIO::SeekNextBlock(const char * id, unsigned short type) const { - int blockHeader[7]; - + cout << "Seek " << type << endl; while (1) { @@ -252,15 +261,30 @@ // - unknown field // - id (first 4 bytes of data field) - fIn->read((char*)blockHeader, 6 * sizeof(int)); + int blockHeader[4]; + fIn->read((char*)blockHeader, 4 * sizeof(int)); if (fIn->eof()) { - gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - unexpected end-of-file." << endl; + gLog << err << "MCorsikaFormatEventIO::SeekNextBlock: ERROR - Unexpected end-of-file seeking " << id << " (" << type << ")." << endl; return kFALSE; } const unsigned short objType = blockHeader[1] & 0xFFFF; + cout << "objType=" << objType << endl; if (type == objType) - break; + { + 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. @@ -272,6 +296,5 @@ // a unknown block, we jump to the next one const int length = blockHeader[3] & 0x3fffffff; - - fIn->seekg(length - 2 * (Int_t)(sizeof(int)), ios::cur); + fIn->seekg(length, ios::cur); } @@ -336,77 +359,82 @@ // returns kFALSE; // -Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, UInt_t telescope) +Int_t MCorsikaFormatEventIO::GetNextEvent(Float_t **buffer, Int_t telescope) { static Float_t *data = NULL; - static Int_t lengthEvent = -1; - static Int_t lengthPhotonData = 0; - - while (lengthPhotonData == 0) + + 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; - - // If we are at the end of the event signal this and - // start next time with a new event - if (lengthEvent==0) + 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) { - lengthEvent=-1; - return kFALSE; // Signal to start with a new event (process task list first) + cout << "skip EVTE." << endl; + return kFALSE; } - // 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; - 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. + 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); - - // Return the pointer to the start of the current photon data - *buffer = data + 3; + *buffer = data; // Return the pointer to the start of the current photon data return kTRUE; +} } @@ -425,5 +453,4 @@ // - identification field // - length - UInt_t blockHeader[4]; fIn->read((char*)blockHeader, 4 * sizeof(Int_t)); @@ -440,4 +467,6 @@ const unsigned short objType = blockHeader[1] & 0xFFFF; + cout << "objType=" << objType << endl; + // Decode length of block length = blockHeader[3] & 0x3fffffff; Index: /trunk/Mars/mcorsika/MCorsikaFormat.h =================================================================== --- /trunk/Mars/mcorsika/MCorsikaFormat.h (revision 9942) +++ /trunk/Mars/mcorsika/MCorsikaFormat.h (revision 9943) @@ -37,6 +37,7 @@ virtual void StorePos(); virtual void ResetPos() const; + virtual void Rewind() const; - virtual Int_t GetNextEvent(Float_t **buffer, UInt_t tel=0) = 0; + virtual Int_t GetNextEvent(Float_t **buffer, Int_t tel=0) = 0; virtual Bool_t IsEventioFormat() const = 0; @@ -44,4 +45,6 @@ std::streampos GetCurrPos() const; + + void Read(void *ptr, Int_t i=4) const; static MCorsikaFormat *CorsikaFormatFactory(const char *fileName); @@ -62,5 +65,5 @@ Bool_t SeekEvtEnd(); - Int_t GetNextEvent(Float_t **buffer, UInt_t); + Int_t GetNextEvent(Float_t **buffer, Int_t); Bool_t IsEventioFormat() const {return kFALSE;} }; @@ -81,5 +84,5 @@ Bool_t SeekEvtEnd(); - Int_t GetNextEvent(Float_t **buffer, UInt_t tel); + Int_t GetNextEvent(Float_t **buffer, Int_t tel); Bool_t IsEventioFormat() const { return kTRUE; } Index: /trunk/Mars/mcorsika/MCorsikaRead.cc =================================================================== --- /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 9942) +++ /trunk/Mars/mcorsika/MCorsikaRead.cc (revision 9943) @@ -96,7 +96,7 @@ // MCorsikaRead::MCorsikaRead(const char *fname, const char *name, const char *title) - : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(0), fForceMode(kFALSE), + : fRunHeader(0), fEvtHeader(0), fEvent(0), fTelescopeIdx(-1), fForceMode(kFALSE), fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0), - fIn(0), fInFormat(0), fParList(0) + /*fIn(0),*/ fInFormat(0), fParList(0), fNumTelescopes(1) { fName = name ? name : "MRead"; @@ -117,6 +117,6 @@ { delete fFileNames; - if (fIn) - delete fIn; +// if (fIn) +// delete fIn; if (fInFormat) delete fInFormat; @@ -191,8 +191,8 @@ // open the input stream and check if it is really open (file exists?) // - if (fIn) +/* if (fIn) delete fIn; fIn = NULL; - +*/ if (fInFormat) delete fInFormat; @@ -239,4 +239,59 @@ // 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(); @@ -319,14 +374,15 @@ } - fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse(); + fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()* + (fTelescopeIdx<0 ? fNumTelescopes : 1); continue; } break; } - +/* if (fIn) delete fIn; fIn = NULL; - +*/ return rc; } @@ -356,8 +412,8 @@ // fParList = 0; - +/* if (!OpenStream()) return kFALSE; - +*/ // // check if all necessary containers exist in the Parameter list. @@ -395,28 +451,48 @@ Int_t MCorsikaRead::ReadEvent() { - // Store the position to re-read a single event several times - // FIXME: Does this work with EventIO, too? - if (fEvtHeader->GetNumReuse()>=fRunHeader->GetNumReuse()-1) - { - fInFormat->StorePos(); - fEvtHeader->ResetNumReuse(); - - if (fInFormat->IsEventioFormat()) + if (fInFormat->IsEventioFormat()) + { + if (fNumTelescope==0) { - // Read the event header 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 { - if (!fInFormat->IsEventioFormat()) + // Store the position to re-read a single event several times + if (fEvtHeader->GetNumReuse()GetNumReuse()-1) fInFormat->ResetPos(); - } - - if (!fInFormat->IsEventioFormat()) - { + else + { + fInFormat->StorePos(); + fEvtHeader->ResetNumReuse(); + } + // Read the event header const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); @@ -424,34 +500,22 @@ 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 - const Int_t id = fInFormat->IsEventioFormat() ? fTelescopeIdx : fEvtHeader->GetNumReuse()+1; - - // 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, 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 - return fEvtHeader->ReadEvtEnd(fInFormat); + 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); + } } @@ -502,5 +566,5 @@ Int_t MCorsikaRead::PostProcess() { - const UInt_t n = fNumEvents*fRunHeader->GetNumReuse(); + const UInt_t n = fNumEvents*fRunHeader->GetNumReuse()*(fTelescopeIdx<0?fNumTelescopes:1); // @@ -517,2 +581,21 @@ return kTRUE; } + +// -------------------------------------------------------------------------- +// +// 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 9942) +++ /trunk/Mars/mcorsika/MCorsikaRead.h (revision 9943) @@ -4,4 +4,8 @@ #ifndef MARS_MRead #include "MRead.h" +#endif + +#ifndef ROOT_TArrayF +#include #endif @@ -22,5 +26,5 @@ MPhotonEvent *fEvent; //! event information - UInt_t fTelescopeIdx; // Index of telescope to be extracted + Int_t fTelescopeIdx; // Index of telescope to be extracted Bool_t fForceMode; // Force mode skipping defect RUNE @@ -30,13 +34,20 @@ UInt_t fNumTotalEvents; //! total number of events in all files - ifstream *fIn; //! input stream (file to read from) +// ifstream *fIn; //! input stream (file to read from) MCorsikaFormat *fInFormat; //! access to input corsika data MParList *fParList; //! tasklist to call ReInit from + Int_t fNumTelescopes; //! Number of telescopes in array + Int_t fNumTelescope; //! Number of telescope currently being read + TArrayF fTelescopeX; //! x pos (measured towards north, unit: cm) + TArrayF fTelescopeY; //! y pos (measured towards west, unit: cm) + TArrayF fTelescopeZ; //! z pos (from detection level, unit: cm) + TArrayF fTelescopeR; //! Radii of spheres around tel. (unit: cm) + //UInt_t fInterleave; //Bool_t fForce; - virtual Bool_t OpenStream() { return kTRUE; } +// virtual Bool_t OpenStream() { return kTRUE; } Bool_t ReadEvtEnd(); @@ -45,7 +56,11 @@ Int_t ReadEvent(); + // MTask Int_t PreProcess(MParList *pList); Int_t Process(); Int_t PostProcess(); + + // MParContainer + Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print); public: @@ -54,5 +69,5 @@ void SetForceMode(Bool_t b=kTRUE) { fForceMode=b; } - void SetTelescopeIdx(UInt_t idx=0) { fTelescopeIdx = idx; } + void SetTelescopeIdx(Int_t idx=-1) { fTelescopeIdx = idx; } TString GetFullFileName() const;