Index: trunk/Mars/msim/MPhotonData.cc =================================================================== --- trunk/Mars/msim/MPhotonData.cc (revision 9949) +++ 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 9949) +++ 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 9949) +++ 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