Ignore:
Timestamp:
11/29/10 14:24:23 (14 years ago)
Author:
rohlfs
Message:
two new command line arguments of readcorsika: -A=arrayNo and -T=telescopeNo. New design of program flow in MCorsikaRead: It is now determined by the order of the data blocks in the input file.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/Mars/msim/MPhotonEvent.cc

    r9949 r10060  
    507507        operator[](i).SimWavelength(wmin, wmax);
    508508}
    509 
    510 // --------------------------------------------------------------------------
    511 //
    512 // Read the Event section from the file
    513 //
    514 Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t id)
    515 {
    516     Int_t n = 0;
    517 
    518     // --- old I/O ---
    519     // Read only + Reflector (no absorption)
    520     // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
    521     // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
    522 
    523     // Read only:
    524     // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
    525     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    526 
    527     // --- new I/O ---
    528     // Read only (don't allocate storage space):
    529     // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
    530     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    531 
    532     // Read only in blocks (with storage allocation):
    533     // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
    534     // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
    535 
    536     // Read only in blocks (without storage allocation):
    537     // similar to just copy
    538 
    539     // Copy with cp
    540     // 1.57GB/ 5s   CPU
    541     // 1.57GB/28s   REAL
    542     // 1.06GB/ 3s   CPU
    543     // 1.06GB/22s   REAL
    544     Float_t *buffer = 0;
    545 
    546     if (fInFormat->IsEventioFormat())
    547     {
    548         while (1)
    549         {
    550             const Int_t rc = fInFormat->GetNextEvent(&buffer, id);
    551             if (rc==kERROR)
    552                 return kERROR;
    553             if (rc==kFALSE)
    554                 break;
    555 
    556             // Loop over number of photons in bunch
    557             while (Add(n).FillEventIO(buffer))
    558                 n++;
    559         }
    560     }
    561     else
    562     {
    563         while (1)
    564         {
    565             const Int_t rc1 = fInFormat->GetNextEvent(&buffer);
    566             if (rc1==kERROR)
    567                 return kERROR;
    568             if (rc1==kFALSE)
    569                 break;
    570 
    571             const Int_t rc2 = Add(n).FillCorsika(buffer, id);
    572             switch (rc2)
    573             {
    574             case kCONTINUE:  continue;        // No data in this bunch... skip it.
    575             case kERROR:     return kERROR;   // Error occured
    576             //case kFALSE:     return kFALSE;   // End of stream
    577             }
    578 
    579             // This is a photon we would like to keep later.
    580             // Increase the counter by one
    581             n++;
    582         }
    583     }
    584 
    585 /*
    586 
    587     while (1)
    588     {
    589         Float_t buffer[273];
    590         Float_t * ptr = buffer;
    591 
    592 
    593         if (!fInFormat->ReadData(273, buffer))
    594             return kFALSE;
    595 
    596         if (!memcmp(ptr, "EVTE", 4))
    597             {
    598 
    599             fInFormat->UnreadLastData();
    600             break;
    601             }
    602 
    603         Float_t *end = ptr + 273;
    604 
    605         // Loop over all sub-blocks (photons) in the block and if
    606         // they contain valid data add them to the array
    607         while (ptr<end)
    608         {
    609             // Get/Add the n-th entry from the array and
    610             // fill it with the current 7 floats
    611             const Int_t rc = Add(n).FillCorsika(ptr);
    612             ptr += 7;
    613 
    614             switch (rc)
    615             {
    616             case kCONTINUE:  continue;        // No data in this bunch... skip it.
    617             case kERROR:     return kERROR;   // Error occured
    618             //case kFALSE:     return kFALSE;   // End of stream
    619             }
    620 
    621             // This is a photon we would like to keep later.
    622             // Increase the counter by one
    623             n++;
    624         }
    625     }
    626 
    627 */
    628     Resize(n);
    629     fData.UnSort();
    630 
    631     SetReadyToSave();
    632 
    633     //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    634     return kTRUE;
    635 
    636 }
    637 
    638 Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin, Int_t i)
    639 {
    640     Int_t n = 0;
    641 
    642     // --- old I/O ---
    643     // Read only + Reflector (no absorption)
    644     // Muons:   1.06GB/115s =  9.2MB/s (100kEvs)
    645     // Gammas:  1.57GB/275s =  5.7MB/s (  1kEvs)
    646 
    647     // Read only:
    648     // Gammas:  1.57GB/158s =  9.9MB/s (  1kEvs)
    649     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    650 
    651     // --- new I/O ---
    652     // Read only (don't allocate storage space):
    653     // Gammas:  1.57GB/143s = 11.0MB/s (  1kEvs)
    654     // Muons:   1.06GB/ 77s = 13.8MB/s (100kEvs)
    655 
    656     // Read only in blocks (with storage allocation):
    657     // Gammas:  1.57GB/28s  =  56MB/s (  1kEvs)
    658     // Muons:   1.06GB/5.2s = 204MB/s (100kEvs)
    659 
    660     // Read only in blocks (without storage allocation):
    661     // similar to just copy
    662 
    663     // Copy with cp
    664     // 1.57GB/ 5s   CPU
    665     // 1.57GB/28s   REAL
    666     // 1.06GB/ 3s   CPU
    667     // 1.06GB/22s   REAL
    668 
    669     while (1)
    670     {
    671         // Stprage for one block
    672         char c[273*4];
    673 
    674         // Read the first four byte to check whether the next block
    675         // doen't belong to the event anymore
    676         fin.read(c, 4);
    677         if (!fin)
    678             return kFALSE;
    679 
    680         // Check if the event is finished
    681         if (!memcmp(c, "EVTE", 4))
    682             break;
    683 
    684         // Now read the rest of the data
    685         fin.read(c+4, 272*4);
    686 
    687         Float_t *ptr = reinterpret_cast<Float_t*>(c);
    688         Float_t *end = ptr + 273;
    689 
    690         // Loop over all sub-blocks (photons) in the block and if
    691         // they contain valid data add them to the array
    692         while (ptr<end)
    693         {
    694             // Get/Add the n-th entry from the array and
    695             // fill it with the current 7 floats
    696             const Int_t rc = Add(n).FillCorsika(ptr, i);
    697             ptr += 7;
    698 
    699             switch (rc)
    700             {
    701             case kCONTINUE:  continue;        // No data in this bunch... skip it.
    702             case kERROR:     return kERROR;   // Error occured
    703             //case kFALSE:     return kFALSE;   // End of stream
    704             }
    705 
    706             // This is a photon we would like to keep later.
    707             // Increase the counter by one
    708             n++;
    709         }
    710     }
    711 /*
    712     while (1)
    713     {
    714         // Check the first four bytes
    715         char c[4];
    716         fin.read(c, 4);
    717 
    718         // End of stream
    719         if (!fin)
    720             return kFALSE;
    721 
    722         // Check if we found the end of the event
    723         if (!memcmp(c, "EVTE", 4))
    724             break;
    725 
    726         // The first for byte contained data already --> go back
    727         fin.seekg(-4, ios::cur);
    728 
    729         // Do not modify this. It is optimized for execution
    730         // speed and flexibility!
    731         MPhotonData &ph = Add(n);
    732         // It checks how many entries the lookup table has. If it has enough
    733         // entries and the entry was already allocated, we can re-use it,
    734         // otherwise we have to allocate it.
    735 
    736         // Now we read a single cherenkov bunch. Note that for speed reason we have not
    737         // called the constructor if the event was already constructed (virtual table
    738         // set), consequently we must make sure that ReadCorsikaEvent does reset
    739         // all data mebers no matter whether they are read or not.
    740         const Int_t rc = ph.ReadCorsikaEvt(fin);
    741 
    742         // Evaluate result from reading event
    743         switch (rc)
    744         {
    745         case kCONTINUE:  continue;        // No data in this bunch... skip it.
    746         case kFALSE:     return kFALSE;   // End of stream
    747         case kERROR:     return kERROR;   // Error occured
    748         }
    749 
    750         // FIXME: If fNumPhotons!=1 add the photon more than once
    751 
    752         // Now increase the number of entries which are kept,
    753         // i.e. keep this photon(s)
    754         n++;
    755     }
    756   */
    757 
    758     Resize(n);
    759     fData.UnSort();
    760 
    761     SetReadyToSave();
    762 
    763     //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    764     return kTRUE;
    765 }
    766 
    767 // --------------------------------------------------------------------------
    768 /*
    769 Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i)
    770 {
    771     Int_t n = 0;
    772 
    773     while (1)
    774     {
    775         // Check the first four bytes
    776         char c[13];
    777         fin.read(c, 13);
    778 
    779         // End of stream
    780         if (!fin)
    781             return kFALSE;
    782 
    783         // Check if we found the end of the event
    784         if (!memcmp(c, "\nEND---EVENT\n", 13))
    785             break;
    786 
    787         // The first for byte contained data already --> go back
    788         fin.seekg(-13, ios::cur);
    789 
    790         // Do not modify this. It is optimized for execution
    791         // speed and flexibility!
    792         //TObject *o = n<fData.GetSize() && fData.UncheckedAt(n) ? fData.UncheckedAt(n) : fData.New(n);
    793 
    794         // Now we read a single cherenkov bunch
    795         //const Int_t rc = static_cast<MPhotonData*>(o)->ReadRflEvt(fin);
    796         const Int_t rc = Add(n).ReadRflEvt(fin, i);
    797 
    798         // Evaluate result from reading event
    799         switch (rc)
    800         {
    801         case kCONTINUE:  continue;        // No data in this bunch... skip it.
    802         case kFALSE:     return kFALSE;   // End of stream
    803         case kERROR:     return kERROR;   // Error occured
    804         }
    805 
    806         // Now increase the number of entries which are kept,
    807         // i.e. keep this photon(s)
    808         n++;
    809     }
    810 
    811     Shrink(n);
    812 
    813     SetReadyToSave();
    814 
    815     // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    816     return kTRUE;
    817 }
    818 */
     509Int_t MPhotonEvent::ReadEventIoEvt(MCorsikaFormat *fInFormat)
     510{
     511   Int_t  bunchHeader[3];
     512   fInFormat->Read(bunchHeader, 3 * sizeof(Int_t));
     513
     514   Int_t n = 0;
     515
     516   for (int bunch = 0; bunch < bunchHeader[2]; bunch++)
     517      {
     518      Float_t buffer[8];
     519      fInFormat->Read(buffer, 8 * sizeof(Float_t));
     520
     521      if (Add(n).FillEventIO(buffer))
     522         n++;
     523      }
     524
     525   Resize(n);
     526   fData.UnSort();
     527
     528   SetReadyToSave();
     529
     530   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     531   return kTRUE;
     532
     533}
     534
     535Int_t MPhotonEvent::ReadCorsikaEvt(Float_t * data, Int_t numEvents, Int_t arrayIdx)
     536{
     537   Int_t n = 0;
     538
     539   for (Int_t event = 0; event < numEvents; event++)
     540      {
     541      const Int_t rc = Add(n).FillCorsika(data + 7 * event, arrayIdx);
     542
     543      switch (rc)
     544      {
     545      case kCONTINUE:  continue;        // No data in this bunch... skip it.
     546      case kERROR:     return kERROR;   // Error occured
     547      }
     548
     549      // This is a photon we would like to keep later.
     550      // Increase the counter by one
     551      n++;
     552      }
     553
     554   Resize(n);
     555   fData.UnSort();
     556
     557   SetReadyToSave();
     558
     559   //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     560   return kTRUE;
     561
     562}
    819563
    820564// --------------------------------------------------------------------------
Note: See TracChangeset for help on using the changeset viewer.