Changeset 9937 for trunk


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
Files:
14 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/Changelog

    r9936 r9937  
    4242   * msim/MHPhotonEvent.cc:
    4343     - implemented types 6-8
     44
     45   * mcorsika/MCorsikaEvtHeader.[h,cc]:
     46     - added fNumReuse
     47     - added member functions to increase and reset the number of reusages
     48     - increased class version number accordingly
     49
     50   * mcorsika/MCorsikaRead.[h,cc]:
     51     - count the number of showers (events) times its reusage
     52       for fNumTotalEvents
     53     - adapted ReadEvent to re-usage of showers (still needs to be tested for
     54       EventIO)
     55     - Fixed return type (should be Int_t) of ReadEvent
     56     - take the number of reusages in PostProcess into account
     57
     58   * mcorsika/MCorsikaRunHeader.[h,cc]:
     59     - added total number of reusages fNumReuse
     60     - increased class version number accordingly
     61
     62   * mmc/MMcEvt.hxx:
     63     - added setter for the reuse counter
     64
     65   * msim/MPhotonData.[h,cc]:
     66     - removed fNumPhotons (we allow 1-photon bunches only anyway and MMCS
     67       doesn't even distribute this number)
     68     - increased class version number accordingly
     69
     70   * msim/MSimMMCS.cc:
     71     - correctly propagate the number of events (number of showers times reusage)
     72     - propagate counter for reusage to MMcEvt
     73
     74   * msim/MPhotonEvent.[h,cc]:
     75     - added functions to calculate mean x and mean y
     76     - propagate re-usage counter through ReadEvent
    4477
    4578
  • trunk/Mars/NEWS

    r9878 r9937  
    6767     mapping of the regular patches, resmc/fact-trigger-all.txt the mapping
    6868     of all patches.
     69
     70   * Added reading of re-used corsika showers (only supported if MMCS is used)
    6971
    7072 ;star:
  • 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
  • trunk/Mars/mmc/MMcEvt.hxx

    r9595 r9937  
    103103
    104104    void SetEvtNumber(UInt_t n) { fEvtNumber=n; }
     105    void SetEventReuse(UInt_t n) { fEventReuse=n; }
    105106    void SetPhotElfromShower(UInt_t n) { fPhotElfromShower=n; }
    106107
  • trunk/Mars/msim/MPhotonData.cc

    r9820 r9937  
    1616!
    1717!
    18 !   Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
     18!   Author(s): Thomas Bretz, 12/2000 <mailto:thomas.bretz@epfl.ch>
    1919!   Author(s): Qi Zhe,       06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
    2020!
    21 !   Copyright: CheObs Software Development, 2000-2009
     21!   Copyright: CheObs Software Development, 2000-2010
    2222!
    2323!
     
    3434//   Version 1:
    3535//   ----------
    36 //    - First implementation
     36//    * First implementation
     37//
     38//   Version 2:
     39//   ----------
     40//    - fNumPhotons
    3741//
    3842/////////////////////////////////////////////////////////////////////////////
     
    5761MPhotonData::MPhotonData(/*const char *name, const char *title*/)
    5862    : fPosX(0), fPosY(0), fCosU(0), fCosV(0), fTime(0), fWavelength(0),
    59     fNumPhotons(1), fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED),
     63    /*fNumPhotons(1),*/ fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED),
    6064    fTag(-1), fWeight(1)
    6165{
     
    8286    MPhotonData &d = static_cast<MPhotonData&>(obj);
    8387
    84     d.fNumPhotons       = fNumPhotons;
     88//    d.fNumPhotons       = fNumPhotons;
    8589    d.fPosX             = fPosX;
    8690    d.fPosY             = fPosY;
     
    185189
    186190    fPrimary    = MMcEvtBasic::kUNDEFINED;
    187     fNumPhotons =  1;
     191//    fNumPhotons =  1;
    188192    fTag        = -1;
    189193    fWeight     =  1;
     
    201205// system intpo our own.
    202206//
    203 Int_t MPhotonData::FillCorsika(Float_t f[7])
     207Int_t MPhotonData::FillCorsika(Float_t f[7], Int_t i)
    204208{
    205209    const UInt_t n = TMath::Nint(f[0]);
    206210    if (n==0)
     211        // FIXME: Do we need to decode the rest anyway?
     212        return kCONTINUE;
     213
     214    // Check reuse
     215    const Int_t reuse = (n/1000)%100; // Force this to be 1!
     216    if (reuse!=i)
    207217        return kCONTINUE;
    208218
     
    210220    fWavelength = n%1000;
    211221    fPrimary    = MMcEvtBasic::ParticleId_t(n/100000);
    212     fNumPhotons = (n/1000)%100; // Force this to be 1!
    213 
    214     if (fNumPhotons!=1)
    215     {
    216         // FIXME: Could be done in MPhotonEvent::ReadCorsikaEvent
    217         gLog << err << "ERROR - MPhotonData::FillCorsika: fNumPhotons not 1, but " << fNumPhotons << endl;
    218         gLog << "        This is not yet supported." << endl;
    219         return kERROR;
    220     }
    221222
    222223    // x=north, y=west
     
    254255Int_t MPhotonData::FillEventIO(Float_t f[8])
    255256{
     257    if (TMath::Nint(f[6])!=1)
     258    {
     259        gLog << err << "ERROR - Bunch sizes != 1 are not supported." << endl;
     260        return kFALSE;
     261    }
     262
    256263    fPosX             =  f[1];              // xpos relative to telescope [cm]
    257264    fPosY             = -f[0];              // ypos relative to telescope [cm]
     
    260267    fTime             =  f[4];              // a relative arival time [ns]
    261268    fProductionHeight =  f[5];              // altitude of emission [cm]
    262     fNumPhotons       =  TMath::Nint(f[6]); // photons in this bunch
     269//    fNumPhotons       =  TMath::Nint(f[6]); // photons in this bunch
    263270    fWavelength       =  TMath::Nint(f[7]); // so far always zeor = unspec. [nm]
    264 
    265271
    266272    // Now reset all data members which are not in the stream
     
    271277    return kTRUE;
    272278}
    273 
     279/*
    274280// --------------------------------------------------------------------------
    275281//
     
    299305    return rc==kTRUE ? !fin.eof() : rc;
    300306}
     307*/
    301308
    302309// --------------------------------------------------------------------------
     
    308315{
    309316    gLog << inf << endl;
    310     gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
     317//    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
     318    gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
    311319    gLog << "Wavelength:       " << fWavelength << "nm" << endl;
    312320    gLog << "Pos X/Y  Cos U/V: " << fPosX << "/" << fPosY << "   " << fCosU << "/" << fCosV << endl;
  • trunk/Mars/msim/MPhotonData.h

    r9616 r9937  
    3737
    3838    Float_t fTime;                       // [ns] Time since first interaction or entrance into atmosphere
    39     // 17M
    4039    UShort_t fWavelength;                // [nm] Wavelength
    41     // 19M
    42     UInt_t   fNumPhotons;                // Number of cherenkov photons ins bunch
     40//    UInt_t   fNumPhotons;                // Number of cherenkov photons ins bunch
    4341    Float_t  fProductionHeight;          // [cm] Height of bunch production
    4442    MMcEvtBasic::ParticleId_t fPrimary;  // Type of emitting particle
    45     // 22M
    46     // gzip
    47     // 25M
    48     // raw
    49     // 32M
    5043
    5144    Int_t   fTag;    //! A tag for external use
     
    129122
    130123    // I/O
    131     Int_t ReadCorsikaEvt(istream &fin);
    132     Int_t ReadRflEvt(istream &fin);
     124    //Int_t ReadCorsikaEvt(istream &fin);
     125    //Int_t ReadRflEvt(istream &fin);
    133126
    134     Int_t FillCorsika(Float_t f[7]);
     127    Int_t FillCorsika(Float_t f[7], Int_t i);
    135128    Int_t FillEventIO(Float_t f[7]);
    136129    Int_t FillRfl(Float_t f[8]);
    137130
    138     ClassDef(MPhotonData, 1) //Container to store a cherenkov photon bunch from a CORSUKA file
     131    ClassDef(MPhotonData, 2) //Container to store a cherenkov photon bunch from a CORSUKA file
    139132};
    140133
  • trunk/Mars/msim/MPhotonEvent.cc

    r9616 r9937  
    449449}
    450450
     451Double_t MPhotonEvent::GetMeanX() const
     452{
     453    const UInt_t n = GetNumPhotons();
     454
     455    Double_t mean = 0;
     456    for (UInt_t i=0; i<n; i++)
     457        mean += operator[](i).GetPosX();
     458
     459    return mean / n;
     460}
     461
     462Double_t MPhotonEvent::GetMeanY() const
     463{
     464    const UInt_t n = GetNumPhotons();
     465
     466    Double_t mean = 0;
     467    for (UInt_t i=0; i<n; i++)
     468        mean += operator[](i).GetPosY();
     469
     470    return mean / n;
     471}
     472
    451473// --------------------------------------------------------------------------
    452474//
    453475// Read the Event section from the file
    454476//
    455 Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat * fInFormat)
     477Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i)
    456478{
    457479    Int_t n = 0;
     
    509531            {
    510532
    511             const Int_t rc = Add(n).FillCorsika(buffer);
     533            const Int_t rc = Add(n).FillCorsika(buffer, i);
    512534            switch (rc)
    513535            {
     
    577599}
    578600
    579 Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin)
     601Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin, Int_t i)
    580602{
    581603    Int_t n = 0;
     
    635657            // Get/Add the n-th entry from the array and
    636658            // fill it with the current 7 floats
    637             const Int_t rc = Add(n).FillCorsika(ptr);
     659            const Int_t rc = Add(n).FillCorsika(ptr, i);
    638660            ptr += 7;
    639661
     
    707729
    708730// --------------------------------------------------------------------------
    709 //
    710 Int_t MPhotonEvent::ReadRflEvt(std::istream &fin)
     731/*
     732Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i)
    711733{
    712734    Int_t n = 0;
     
    735757        // Now we read a single cherenkov bunch
    736758        //const Int_t rc = static_cast<MPhotonData*>(o)->ReadRflEvt(fin);
    737         const Int_t rc = Add(n).ReadRflEvt(fin);
     759        const Int_t rc = Add(n).ReadRflEvt(fin, i);
    738760
    739761        // Evaluate result from reading event
     
    754776    SetReadyToSave();
    755777
    756     //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     778    // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
    757779    return kTRUE;
    758780}
    759 
     781*/
    760782// --------------------------------------------------------------------------
    761783//
  • trunk/Mars/msim/MPhotonEvent.h

    r9616 r9937  
    3737    Double_t GetTimeMedianDev() const;
    3838
     39    Double_t GetMeanX() const;
     40    Double_t GetMeanY() const;
     41
    3942    TClonesArray &GetArray() { return fData; }
    4043    const TClonesArray &GetArray() const { return fData; }
     
    5356
    5457    // I/O
    55     Int_t ReadCorsikaEvt(MCorsikaFormat * fInFormat);
    56     Int_t ReadCorsikaEvt(istream &fin);
    57     Int_t ReadRflEvt(istream &fin);
     58    Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i);
     59    Int_t ReadCorsikaEvt(istream &fin, Int_t i);
     60    //Int_t ReadRflEvt(istream &fin, Int_t i);
    5861
    5962    // TObject
  • trunk/Mars/msim/MSimMMCS.cc

    r9595 r9937  
    131131
    132132    // FIXME: Is there a way to write them as LAST entry in the file?
    133     fMcRunHeader->SetNumSimulatedShowers(fRunHeader->GetNumEvents());
     133    fMcRunHeader->SetNumSimulatedShowers(fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse());
    134134    fMcRunHeader->SetCorsikaVersion(TMath::Nint(fRunHeader->GetProgramVersion()*100));
    135135
     
    211211
    212212    fMcEvt->SetEvtNumber(fEvtHeader->GetEvtNumber());
     213    fMcEvt->SetEventReuse(fEvtHeader->GetNumReuse());
    213214    fMcEvt->SetPhotElfromShower(0);
    214215
Note: See TracChangeset for help on using the changeset viewer.