Changeset 9937 for trunk/Mars/mcorsika


Ignore:
Timestamp:
09/22/10 15:25:55 (14 years ago)
Author:
tbretz
Message:
Added support for re-use of showers in MMCS.
Location:
trunk/Mars/mcorsika
Files:
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/Mars/mcorsika/MCorsikaEvtHeader.cc

    r9616 r9937  
    3232//
    3333//
     34// Class Version 3:
     35// ----------------
     36//  - UInt_t fNumReuse
     37//
     38//
    3439/////////////////////////////////////////////////////////////////////////////
    3540#include "MCorsikaEvtHeader.h"
    … …  
    5560//
    5661MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title)
    57     : fX(0), fY(0)
     62    : fNumReuse((UInt_t)-1), fX(0), fY(0)
    5863{
    5964    fName  = name  ? name  : "MCorsikaEvtHeader";
    … …  
    104109{
    105110    *fLog << all;
    106     *fLog << "Event Number:              " << dec << fEvtNumber << endl;
     111    *fLog << "Event Number:              " << dec << fEvtNumber << " (reused=" << fNumReuse << ")" << endl;
    107112//    *fLog << "Particle ID:               " << MMcEvt::GetParticleName(fParticleID) << endl;
    108113    *fLog << "Energy:                    " << fTotalEnergy << "GeV" << endl;
    … …  
    122127// return FALSE if there is no  header anymore, else TRUE
    123128//
    124 Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat * fInFormat)
     129Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat *fInFormat)
    125130{
    126 
    127131    if (!fInFormat->SeekNextBlock("EVTH", 1202))
    128132        return kFALSE;
    … …  
    151155    fMomentumZ = -f[8];  // Does this minus make sense?!
    152156
    153     // FIXME: Correct for direction of magnetic field!
    154157    fZd        = f[9];
    155158    fAz        = TMath::Pi()-f[10];
    156159
    157     const Int_t n = TMath::Nint(f[96]);
    158     if (n!=1)
    159     {
    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;
    162     }
     160    // f[96] // Number of reuse of event [max=20]
    163161
    164     fX =  f[117];   //fX = f[97];
    165     fY = -f[97];    //fY = f[117];
     162    // FIXME: Check fNumReuse<20
     163    fX =  f[117 + fNumReuse];
     164    fY = -f[ 97 + fNumReuse];
     165
     166    fWeightedNumPhotons = 0;
    166167
    167168    return !fInFormat->Eof();
    … …  
    171172// this member function is for reading the event end block
    172173
    173 Bool_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
     174Int_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat)
    174175{
    175 
    176176    if (!fInFormat->SeekNextBlock("EVTE", 1209))
    177177        return kFALSE;
    … …  
    192192    }
    193193
     194    // FIXME: What's the meaning of this if in case of reusing the event this number
     195    //        exists only once?
    194196    fWeightedNumPhotons = f[1];
    195197
    196198    return !fInFormat->Eof();
    197199}
    198 

  • trunk/Mars/mcorsika/MCorsikaEvtHeader.h

    r9616 r9937  
    2020private:
    2121    UInt_t   fEvtNumber;              // Event number
     22    UInt_t   fNumReuse;               // Counter of the reuse of the same shower
    2223//    UInt_t   fParticleID;             // Particle ID (see MMcEvtBasic or CORSIKA manual)
    2324    Float_t  fTotalEnergy;            // [GeV]
    … …  
    4647
    4748    UInt_t GetEvtNumber() const { return fEvtNumber; }
     49    UInt_t GetNumReuse() const { return fNumReuse; }
    4850//    UInt_t GetParticleID() const { return fParticleID; }
    4951
    … …  
    6264    Double_t GetImpact() const;
    6365
    64     Int_t  ReadEvt(MCorsikaFormat * fInFormat);    // read in event header block
    65     Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); // read in event end block
     66    void IncNumReuse() { fNumReuse++; }
     67    void ResetNumReuse() { fNumReuse=(UInt_t)-1; }
    6668
    67     ClassDef(MCorsikaEvtHeader, 1) // Parameter Conatiner for raw EVENT HEADER
     69    Int_t ReadEvt(MCorsikaFormat *informat);    // read in event header block
     70    Int_t ReadEvtEnd(MCorsikaFormat *informat); // read in event end block
     71
     72    ClassDef(MCorsikaEvtHeader, 3) // Parameter Conatiner for raw EVENT HEADER
    6873};
    6974

  • trunk/Mars/mcorsika/MCorsikaRead.cc

    r9616 r9937  
    319319            }
    320320
    321             fNumTotalEvents += fRunHeader->GetNumEvents();
     321            fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse();
    322322            continue;
    323323        }
    … …  
    393393// Read a single event from the stream
    394394//
    395 Bool_t MCorsikaRead::ReadEvent()
    396 {
    397     //
    398     // Read in the next EVENT HEADER (see specification),
    399     // if there is no next event anymore stop eventloop
    400     //
    401     Int_t rc = fEvtHeader->ReadEvt(fInFormat); //read event header block
    402     if (!rc)
    403         return kFALSE;
    404 
    405     rc = fEvent->ReadCorsikaEvt(fInFormat);
    406 
    407     /*
    408     // Check if we are allowed to stop reading in the middle of the data
    409     if (rc==kFALSE && !fForce)
    410     {
    411         *fLog << err;
    412         *fLog << "ERROR - End of file in the middle of the data stream!" << endl;
    413         *fLog << "        Set the force-option to force processing." << endl;
    414         return kERROR;
    415     }
    416     */
    417 
    418     return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fInFormat) : rc;
     395Int_t MCorsikaRead::ReadEvent()
     396{
     397    // Store the position to re-read a single event several times
     398    //  FIXME: Does this work with EventIO, too?
     399    if (fEvtHeader->GetNumReuse()>=fRunHeader->GetNumReuse()-1)
     400    {
     401        fInFormat->StorePos();
     402        fEvtHeader->ResetNumReuse();
     403    }
     404    else
     405        fInFormat->ResetPos();
     406
     407    fEvtHeader->IncNumReuse();
     408
     409    // Read the event header
     410    const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat);
     411    if (rc1!=kTRUE)
     412        return rc1;
     413
     414    // Check if reusage number should be reset and increase it
     415    //  Note, that the trick here is that after reset it is set to -1
     416
     417    // Read the photons corresponding to the i-th core location
     418    const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1);
     419    if (rc2!=kTRUE)
     420        return rc2;
     421
     422    // read event end
     423    return fEvtHeader->ReadEvtEnd(fInFormat);
    419424}
    420425
    … …  
    434439        {
    435440            // Read a single event from file
    436             const Bool_t rc = ReadEvent();
     441            const Int_t rc = ReadEvent();
    437442
    438443            // kFALSE means: end of file (try next one)
    … …  
    465470Int_t MCorsikaRead::PostProcess()
    466471{
     472    const UInt_t n = fNumEvents*fRunHeader->GetNumReuse();
     473
    467474    //
    468475    // Sanity check for the number of events
    469476    //
    470     if (fNumEvents==GetNumExecutions()-1 || GetNumExecutions()==0)
     477    if (n==GetNumExecutions()-1 || GetNumExecutions()==0)
    471478        return kTRUE;
    472479
    … …  
    474481    *fLog << "Warning - number of read events (" << GetNumExecutions()-1;
    475482    *fLog << ") doesn't match number in run header(s) (";
    476     *fLog << fNumEvents << ")." << endl;
     483    *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl;
    477484
    478485    return kTRUE;

  • trunk/Mars/mcorsika/MCorsikaRead.h

    r9616 r9937  
    2929    UInt_t    fNumTotalEvents; //! total number of events in all files
    3030
    31     ifstream *fIn;             //! input stream (file to read from)
    32          MCorsikaFormat * fInFormat; //! access to input corsika data
     31    ifstream       *fIn;       //! input stream (file to read from)
     32    MCorsikaFormat *fInFormat; //! access to input corsika data
    3333
    3434    MParList *fParList;        //! tasklist to call ReInit from
    … …  
    4242    Int_t  OpenNextFile(Bool_t print=kTRUE);
    4343    Bool_t CalcNumTotalEvents();
    44     Bool_t ReadEvent();
     44    Int_t ReadEvent();
    4545
    4646    Int_t PreProcess(MParList *pList);

  • trunk/Mars/mcorsika/MCorsikaRunHeader.cc

    r9616 r9937  
    4242//  + Float_t fAtmosphericCoeffC[5]
    4343//  + UInt_t  fCerenkovFlag
     44//
     45// Class Version 3:
     46// ----------------
     47//  + UInt_t  fNumReuse
    4448//
    4549////////////////////////////////////////////////////////////////////////////
    … …  
    164168    fInFormat->UnreadLastHeader();
    165169
    166     const Int_t n = TMath::Nint(g[96]);  // Number i of uses of each cherenkov event
    167     if (n!=1)
     170    fNumReuse = TMath::Nint(g[96]);  // Number i of uses of each cherenkov event
     171    if (fNumReuse!=1)
    168172    {
    169173        *fLog << err  << "ERROR   - Currently only one impact parameter per event is supported." << endl;
    … …  
    269273    *fLog << "Particle ID:    " << MMcEvt::GetParticleName(fParticleID) << endl;
    270274    if (fNumEvents>0)
    271         *fLog << "Num Events:     " << fNumEvents << endl;
     275        *fLog << "Num Events:     " << fNumEvents << " (reuse " << fNumReuse << " times)" << endl;
    272276    *fLog << "Obs Level:     ";
    273277    for (Byte_t i=0; i<fNumObsLevel; i++)

  • trunk/Mars/mcorsika/MCorsikaRunHeader.h

    r9616 r9937  
    3333
    3434    UInt_t  fRunNumber;               // Run number
     35    UInt_t  fNumReuse;                // Number of how many times the same shower is used
    3536    UInt_t  fParticleID;              // Particle ID (see MMcEvtBasic or CORSIKA manual)
    3637    UInt_t  fNumEvents;               // Number of events
    … …  
    7677    UInt_t GetParticleID() const { return fParticleID; }
    7778    UInt_t GetNumEvents() const { return fNumEvents; }
     79    UInt_t GetNumReuse() const { return fNumReuse; }
    7880
    7981    const MTime &GetRunStart() const { return fRunStart; }
    … …  
    128130    void Print(Option_t *t=NULL) const;
    129131
    130     ClassDef(MCorsikaRunHeader, 2)      // storage container for general info
     132    ClassDef(MCorsikaRunHeader, 3)      // storage container for general info
    131133};
    132134#endif
Note: See TracChangeset for help on using the changeset viewer.