Changeset 9616 for trunk/MagicSoft


Ignore:
Timestamp:
07/23/10 10:48:12 (14 years ago)
Author:
rohlfs
Message:
can read also EventIO files
Location:
trunk/MagicSoft/Mars
Files:
2 added
11 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.cc

    r9571 r9616  
    3434/////////////////////////////////////////////////////////////////////////////
    3535#include "MCorsikaEvtHeader.h"
     36#include "MCorsikaFormat.h"
    3637
    3738#include <iomanip>
    … …  
    121122// return FALSE if there is no  header anymore, else TRUE
    122123//
    123 Int_t MCorsikaEvtHeader::ReadEvt(std::istream &fin)
     124Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat * fInFormat)
    124125{
    125     char evth[4];
    126     fin.read(evth, 4);
    127     if (memcmp(evth, "EVTH", 4))
    128     {
    129         fin.seekg(-4, ios::cur);
     126
     127    if (!fInFormat->SeekNextBlock("EVTH", 1202))
    130128        return kFALSE;
    131     }
    132129
    133130    Float_t f[273];
    134     fin.read((char*)&f, 273*4);
     131    if (!fInFormat->ReadData(272, f))
     132        return kFALSE;
    135133
    136134    fEvtNumber  = TMath::Nint(f[0]);
    … …  
    160158    if (n!=1)
    161159    {
    162         *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl;
    163         return kFALSE;
     160        *fLog << err  << "ERROR  - Currently only one impact parameter per event is supported." << endl;
     161        *fLog << warn << "WARNING - This error is replaced by a warning." << endl;
    164162    }
    165163
    … …  
    167165    fY = -f[97];    //fY = f[117];
    168166
    169     fin.seekg(1088-273*4, ios::cur);
    170 
    171     return !fin.eof();
     167    return !fInFormat->Eof();
    172168}
    173169
    … …  
    175171// this member function is for reading the event end block
    176172
    177 Bool_t MCorsikaEvtHeader::ReadEvtEnd(std::istream &fin)
     173Bool_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
    178174{
     175
     176    if (!fInFormat->SeekNextBlock("EVTE", 1209))
     177        return kFALSE;
     178
    179179    //fin.seekg(-1088,ios::cur);
    180180
    181181    Float_t f[2];
    182     fin.read((char*)&f, 2*4);
     182
     183    if (!fInFormat->ReadData(2, f))
     184        return kFALSE;
    183185
    184186    const UInt_t evtnum = TMath::Nint(f[0]);
    … …  
    192194    fWeightedNumPhotons = f[1];
    193195
    194     fin.seekg(1080,ios::cur);
    195 
    196     return !fin.eof();
     196    return !fInFormat->Eof();
    197197}
    198198

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaEvtHeader.h

    r9331 r9616  
    1313//class ifstream;
    1414#include <iosfwd>
     15
     16class MCorsikaFormat;
    1517
    1618class MCorsikaEvtHeader : public MParContainer
    … …  
    6062    Double_t GetImpact() const;
    6163
    62     Int_t  ReadEvt(istream& fin);    // read in event header block
    63     Bool_t ReadEvtEnd(istream& fin); // read in event end block
     64    Int_t  ReadEvt(MCorsikaFormat * fInFormat);    // read in event header block
     65    Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); // read in event end block
    6466
    6567    ClassDef(MCorsikaEvtHeader, 1) // Parameter Conatiner for raw EVENT HEADER

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.cc

    r9314 r9616  
    4949#include "MStatusDisplay.h"
    5050
     51#include "MCorsikaFormat.h"
    5152#include "MCorsikaRunHeader.h"
    5253#include "MCorsikaEvtHeader.h"
    … …  
    9798    : fRunHeader(0), fEvtHeader(0), fEvent(0), /*fEvtData(0),*/ fForceMode(kFALSE),
    9899    fFileNames(0), fNumFile(0), fNumEvents(0), fNumTotalEvents(0),
    99     fIn(0), fParList(0)
     100    fIn(0),  fInFormat(0), fParList(0)
    100101{
    101102    fName  = name  ? name  : "MRead";
    … …  
    118119    if (fIn)
    119120        delete fIn;
     121    if (fInFormat)
     122        delete fInFormat;
    120123}
    121124
    … …  
    161164Bool_t MCorsikaRead::ReadEvtEnd()
    162165{
    163     if (!fRunHeader->SeekEvtEnd(*fIn))
     166    if (!fInFormat->SeekEvtEnd())
    164167    {
    165168        *fLog << (fForceMode?warn:err) << "Error: RUNE section not found in file." << endl;
    … …  
    168171    }
    169172
    170     if (!fRunHeader->ReadEvtEnd(*fIn))
     173    if (!fRunHeader->ReadEvtEnd(fInFormat))
    171174    {
    172175        *fLog << (fForceMode?warn:err) << "Error: Reading RUNE section failed." << endl;
    … …  
    184187Int_t MCorsikaRead::OpenNextFile(Bool_t print)
    185188{
     189
    186190    //
    187191    // open the input stream and check if it is really open (file exists?)
    … …  
    191195    fIn = NULL;
    192196
     197    if (fInFormat)
     198       delete fInFormat;
     199    fInFormat = NULL;
     200
    193201    //
    194202    // Check for the existance of a next file to read
    … …  
    204212
    205213    const char *expname = gSystem->ExpandPathName(name);
    206     fIn = new ifstream(expname);
    207 
    208     const Bool_t noexist = !(*fIn);
    209     if (noexist)
    210     {
    211         *fLog << err << "Cannot open file " << expname << ": ";
    212         *fLog << (errno!=0?strerror(errno):"Insufficient memory for decompression") << endl;
    213     }
    214     else
    215     {
    216         *fLog << inf << "Open file: '" << name << "'" << endl;
    217 
    218         if (fDisplay)
    219         {
    220             // Show the number of the last event after
    221             // which we now open a new file
    222             TString txt = GetFileName();
    223             txt += " @ ";
    224             txt += GetNumExecutions()-1;
    225             fDisplay->SetStatusLine2(txt);
    226         }
    227     }
    228 
     214    fInFormat = CorsikaFormatFactory(fLog, expname);
    229215    delete [] expname;
    230216
    231     if (noexist)
     217    if (fInFormat == NULL)
    232218        return kERROR;
    233219
     220    *fLog << inf << "Open file: '" << name << "'" << endl;
     221
     222    if (fDisplay)
     223    {
     224       // Show the number of the last event after
     225       // which we now open a new file
     226       TString txt = GetFileName();
     227       txt += " @ ";
     228       txt += GetNumExecutions()-1;
     229       fDisplay->SetStatusLine2(txt);
     230    }
     231
    234232    fNumFile++;
    235233
    … …  
    237235    // Read RUN HEADER (see specification) from input stream
    238236    //
    239     if (!fRunHeader->ReadEvt(*fIn))
     237    if (!fRunHeader->ReadEvt(fInFormat))
    240238        return kERROR;
    241239//    if (!fEvtHeader->ReadRunHeader(*fIn, *fRunHeader))
    242240//        return kERROR;
    243241
    244     const streampos pos = fIn->tellg();
     242    fInFormat->StorePos();
    245243    if (!ReadEvtEnd())
    246244        return kERROR;
    247     fIn->seekg(pos, ios::beg);
     245    fInFormat->ResetPos();
    248246
    249247
    … …  
    395393// Read a single event from the stream
    396394//
    397 Bool_t MCorsikaRead::ReadEvent(istream &fin)
     395Bool_t MCorsikaRead::ReadEvent()
    398396{
    399397    //
    … …  
    401399    // if there is no next event anymore stop eventloop
    402400    //
    403     Int_t rc = fEvtHeader->ReadEvt(fin); //read event header block
     401    Int_t rc = fEvtHeader->ReadEvt(fInFormat); //read event header block
    404402    if (!rc)
    405403        return kFALSE;
    406404
    407     rc = fEvent->ReadCorsikaEvt(fin);
     405    rc = fEvent->ReadCorsikaEvt(fInFormat);
    408406
    409407    /*
    … …  
    418416    */
    419417
    420     return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fin) : rc;
     418    return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fInFormat) : rc;
    421419}
    422420
    … …  
    433431    while (1)
    434432    {
    435         if (fIn)
     433        if (fInFormat)
    436434        {
    437435            // Read a single event from file
    438             const Bool_t rc = ReadEvent(*fIn);
     436            const Bool_t rc = ReadEvent();
    439437
    440438            // kFALSE means: end of file (try next one)
    … …  
    442440                return rc;
    443441
    444             if (!fRunHeader->ReadEvtEnd(*fIn))
     442
     443            fInFormat->UnreadLastHeader();
     444            if (!fRunHeader->ReadEvtEnd(fInFormat))
    445445                if (!fForceMode)
    446446                    return kERROR;

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaRead.h

    r9518 r9616  
    1313class MCorsikaEvtHeader;
    1414class MPhotonEvent;
     15class MCorsikaFormat;
    1516
    1617class MCorsikaRead : public MRead
    … …  
    2930
    3031    ifstream *fIn;             //! input stream (file to read from)
     32         MCorsikaFormat * fInFormat; //! access to input corsika data
    3133
    3234    MParList *fParList;        //! tasklist to call ReInit from
    … …  
    4042    Int_t  OpenNextFile(Bool_t print=kTRUE);
    4143    Bool_t CalcNumTotalEvents();
    42     Bool_t ReadEvent(istream &fin);
     44    Bool_t ReadEvent();
    4345
    4446    Int_t PreProcess(MParList *pList);

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.cc

    r9378 r9616  
    4646
    4747#include "MCorsikaRunHeader.h"
     48#include "MCorsikaFormat.h"
    4849
    4950#include <fstream>
    … …  
    7879// Read in one run header from the binary file
    7980//
    80 Bool_t MCorsikaRunHeader::ReadEvt(istream& fin)
    81 {
    82     char runh[4];
    83     fin.read(runh, 4);
    84     if (memcmp(runh, "RUNH", 4))
    85     {
    86         *fLog << err << "ERROR - Wrong identifier: RUNH expected." << endl;
    87         return kFALSE;
    88     }
    89 
    90     Float_t f[272*4];
    91     fin.read((char*)f, 272*4);
     81Bool_t MCorsikaRunHeader::ReadEvt(MCorsikaFormat * fInFormat)
     82{
     83    if (!fInFormat->SeekNextBlock("RUNH", 1200))
     84        return kFALSE;
     85
     86    Float_t f[272];
     87    if (!fInFormat->ReadData(272, f))
     88        return kFALSE;
    9289
    9390    fRunNumber = TMath::Nint(f[0]);
    … …  
    157154    // f[145] Muon multiple scattering flag
    158155
    159     char evth[4];
    160     fin.read(evth, 4);
    161     if (memcmp(evth, "EVTH", 4))
    162     {
    163         *fLog << err << "ERROR - Wrong identifier: EVTH expected." << endl;
    164         return kFALSE;
    165     }
     156    if (!fInFormat->SeekNextBlock("EVTH", 1202))
     157        return kFALSE;
    166158
    167159    Float_t g[273];
    168     fin.read((char*)&g, 273*4);
    169     if (fin.eof())
    170         return kFALSE;
    171 
    172     fin.seekg(-274*4, ios::cur);
     160    if (!fInFormat->ReadData(272, g))
     161        return kFALSE;
     162
     163    fInFormat->UnreadLastData();
     164    fInFormat->UnreadLastHeader();
    173165
    174166    const Int_t n = TMath::Nint(g[96]);  // Number i of uses of each cherenkov event
    175167    if (n!=1)
    176168    {
    177         *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl;
    178         return kFALSE;
     169        *fLog << err  << "ERROR  - Currently only one impact parameter per event is supported." << endl;
     170        *fLog << warn << "WARNING - This error is replaced by a warning." << endl;
    179171    }
    180172
    … …  
    223215}
    224216
    225 Bool_t MCorsikaRunHeader::ReadEvtEnd(istream& fin)
    226 {
    227     char runh[4];
    228     fin.read(runh, 4);
    229     if (memcmp(runh, "RUNE", 4))
    230     {
    231         *fLog << err << "ERROR - Wrong identifier: RUNE expected." << endl;
    232         return kFALSE;
    233     }
     217Bool_t MCorsikaRunHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
     218{
     219
     220    if (!fInFormat->SeekNextBlock("RUNE", 1210))
     221        return kFALSE;
    234222
    235223    Float_t f[2];
    236     fin.read((char*)f, 2*4);
     224    if (!fInFormat->ReadData(2, f))
     225        return kFALSE;
    237226
    238227    const UInt_t runnum = TMath::Nint(f[0]);
    … …  
    245234
    246235    fNumEvents = TMath::Nint(f[1]);
    247 
    248     fin.seekg(270*4, ios::cur);     // skip the remaining data of this block
    249236
    250237    return kTRUE;

  • trunk/MagicSoft/Mars/mcorsika/MCorsikaRunHeader.h

    r9595 r9616  
    1010#include "MTime.h"
    1111#endif
     12
     13class MCorsikaFormat;
    1214
    1315class MCorsikaRunHeader : public MParContainer
    … …  
    119121
    120122    // I/O
    121     Bool_t ReadEvt(istream& fin);
    122     Bool_t ReadEvtEnd(istream& fin);
     123    Bool_t ReadEvt(MCorsikaFormat * fInFormat);
     124    Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat);
    123125    Bool_t SeekEvtEnd(istream &fin);
    124126

  • trunk/MagicSoft/Mars/mcorsika/Makefile

    r9518 r9616  
    2222
    2323SRCFILES = MCorsikaRunHeader.cc \
     24           MCorsikaFormat.cc \
    2425           MCorsikaEvtHeader.cc \
    2526           MCorsikaRead.cc

  • trunk/MagicSoft/Mars/msim/MPhotonData.cc

    r9348 r9616  
    245245// --------------------------------------------------------------------------
    246246//
     247// Set the data member according to the 8 floats read from a eventio-file.
     248// This function MUST reset all data-members, no matter whether these are
     249// contained in the input stream.
     250//
     251// Currently we exchange x and y and set y=-y to convert Corsikas coordinate
     252// system into our own.
     253//
     254Int_t MPhotonData::FillEventIO(Float_t f[8])
     255{
     256    fPosX             =  f[1]; // xpos relative to telescope [cm]
     257    fPosY             = -f[0]; // ypos relative to telescope [cm]
     258    fCosU             =  f[3]; // cos to x
     259    fCosV             = -f[2]; // cos to y
     260    fTime             =  f[4]; // a relative arival time [ns]
     261    fProductionHeight =  f[5]; // altitude of emission [cm]
     262    fNumPhotons       =  f[6]; // photons in this bunch
     263    fWavelength       =  f[7]; // so far always zeor = unspec. [nm]
     264
     265
     266    // Now reset all data members which are not in the stream
     267    fPrimary    = MMcEvtBasic::kUNDEFINED;
     268    fTag        = -1;
     269    fWeight     =  1;
     270
     271    return kTRUE;
     272}
     273
     274// --------------------------------------------------------------------------
     275//
    247276// Read seven floats from the stream and call FillCorsika for them.
    248277//

  • trunk/MagicSoft/Mars/msim/MPhotonData.h

    r9370 r9616  
    133133
    134134    Int_t FillCorsika(Float_t f[7]);
     135    Int_t FillEventIO(Float_t f[7]);
    135136    Int_t FillRfl(Float_t f[8]);
    136137

  • trunk/MagicSoft/Mars/msim/MPhotonEvent.cc

    r9349 r9616  
    121121
    122122#include "MPhotonData.h"
     123#include "MCorsikaFormat.h"
    123124
    124125ClassImp(MPhotonEvent);
    … …  
    452453// Read the Event section from the file
    453454//
    454 Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin)
     455Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat * fInFormat)
    455456{
    456457    Int_t n = 0;
    … …  
    482483    // 1.06GB/ 3s   CPU
    483484    // 1.06GB/22s   REAL
     485    Bool_t readError = kFALSE;
     486    Float_t * buffer;
     487
     488    if ( fInFormat->IsEventioFormat() )
     489        {
     490        while (fInFormat->GetNextEvent(&buffer, readError))
     491            {
     492
     493            const Int_t rc = Add(n).FillEventIO(buffer);
     494            switch (rc)
     495            {
     496            case kCONTINUE:  continue;        // No data in this bunch... skip it.
     497            case kERROR:     return kERROR;   // Error occured
     498            //case kFALSE:     return kFALSE;   // End of stream
     499            }
     500
     501            // This is a photon we would like to keep later.
     502            // Increase the counter by one
     503            n++;
     504            }
     505        }
     506    else
     507        {
     508        while (fInFormat->GetNextEvent(&buffer, readError))
     509            {
     510
     511            const Int_t rc = Add(n).FillCorsika(buffer);
     512            switch (rc)
     513            {
     514            case kCONTINUE:  continue;        // No data in this bunch... skip it.
     515            case kERROR:     return kERROR;   // Error occured
     516            //case kFALSE:     return kFALSE;   // End of stream
     517            }
     518
     519            // This is a photon we would like to keep later.
     520            // Increase the counter by one
     521            n++;
     522            }
     523        }
     524     if (readError)      return kFALSE;
     525
     526/*
    484527
    485528    while (1)
    486529    {
    487         // Stprage for one block
    488         char c[273*4];
    489 
    490         // Read the first four byte to check whether the next block
    491         // doen't belong to the event anymore
    492         fin.read(c, 4);
    493         if (!fin)
     530        Float_t buffer[273];
     531        Float_t * ptr = buffer;
     532
     533
     534        if (!fInFormat->ReadData(273, buffer))
    494535            return kFALSE;
    495536
    496         // Check if the event is finished
    497         if (!memcmp(c, "EVTE", 4))
     537        if (!memcmp(ptr, "EVTE", 4))
     538            {
     539
     540            fInFormat->UnreadLastData();
    498541            break;
    499 
    500         // Now read the rest of the data
    501         fin.read(c+4, 272*4);
    502 
    503         Float_t *ptr = reinterpret_cast<Float_t*>(c);
     542            }
     543
    504544        Float_t *end = ptr + 273;
    505545
    … …  
    525565        }
    526566    }
     567
     568*/
     569    Resize(n);
     570    fData.UnSort();
     571
     572    SetReadyToSave();
     573
     574    //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     575    return kTRUE;
     576
     577}
     578
     579Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin)
     580{
     581    Int_t n = 0;
     582
     583    // --- old I/O ---
     584    // Read only + Reflector (no absorption)
     585    // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
     586    // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
     587
     588    // Read only:
     589    // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
     590    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
     591
     592    // --- new I/O ---
     593    // Read only (don't allocate storage space):
     594    // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
     595    // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
     596
     597    // Read only in blocks (with storage allocation):
     598    // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
     599    // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
     600
     601    // Read only in blocks (without storage allocation):
     602    // similar to just copy
     603
     604    // Copy with cp
     605    // 1.57GB/ 5s   CPU
     606    // 1.57GB/28s   REAL
     607    // 1.06GB/ 3s   CPU
     608    // 1.06GB/22s   REAL
     609
     610    while (1)
     611    {
     612        // Stprage for one block
     613        char c[273*4];
     614
     615        // Read the first four byte to check whether the next block
     616        // doen't belong to the event anymore
     617        fin.read(c, 4);
     618        if (!fin)
     619            return kFALSE;
     620
     621        // Check if the event is finished
     622        if (!memcmp(c, "EVTE", 4))
     623            break;
     624
     625        // Now read the rest of the data
     626        fin.read(c+4, 272*4);
     627
     628        Float_t *ptr = reinterpret_cast<Float_t*>(c);
     629        Float_t *end = ptr + 273;
     630
     631        // Loop over all sub-blocks (photons) in the block and if
     632        // they contain valid data add them to the array
     633        while (ptr<end)
     634        {
     635            // Get/Add the n-th entry from the array and
     636            // fill it with the current 7 floats
     637            const Int_t rc = Add(n).FillCorsika(ptr);
     638            ptr += 7;
     639
     640            switch (rc)
     641            {
     642            case kCONTINUE:  continue;        // No data in this bunch... skip it.
     643            case kERROR:     return kERROR;   // Error occured
     644            //case kFALSE:     return kFALSE;   // End of stream
     645            }
     646
     647            // This is a photon we would like to keep later.
     648            // Increase the counter by one
     649            n++;
     650        }
     651    }
    527652/*
    528653    while (1)

  • trunk/MagicSoft/Mars/msim/MPhotonEvent.h

    r9349 r9616  
    1616class MPhotonData;
    1717class MCorsikaRunHeader;
     18class MCorsikaFormat;
    1819
    1920class MPhotonEvent : public MParContainer
    … …  
    5253
    5354    // I/O
     55    Int_t ReadCorsikaEvt(MCorsikaFormat * fInFormat);
    5456    Int_t ReadCorsikaEvt(istream &fin);
    5557    Int_t ReadRflEvt(istream &fin);
Note: See TracChangeset for help on using the changeset viewer.