Index: trunk/Mars/Changelog =================================================================== --- trunk/Mars/Changelog (revision 9936) +++ trunk/Mars/Changelog (revision 9937) @@ -42,4 +42,37 @@ * msim/MHPhotonEvent.cc: - implemented types 6-8 + + * mcorsika/MCorsikaEvtHeader.[h,cc]: + - added fNumReuse + - added member functions to increase and reset the number of reusages + - increased class version number accordingly + + * mcorsika/MCorsikaRead.[h,cc]: + - count the number of showers (events) times its reusage + for fNumTotalEvents + - adapted ReadEvent to re-usage of showers (still needs to be tested for + EventIO) + - Fixed return type (should be Int_t) of ReadEvent + - take the number of reusages in PostProcess into account + + * mcorsika/MCorsikaRunHeader.[h,cc]: + - added total number of reusages fNumReuse + - increased class version number accordingly + + * mmc/MMcEvt.hxx: + - added setter for the reuse counter + + * msim/MPhotonData.[h,cc]: + - removed fNumPhotons (we allow 1-photon bunches only anyway and MMCS + doesn't even distribute this number) + - increased class version number accordingly + + * msim/MSimMMCS.cc: + - correctly propagate the number of events (number of showers times reusage) + - propagate counter for reusage to MMcEvt + + * msim/MPhotonEvent.[h,cc]: + - added functions to calculate mean x and mean y + - propagate re-usage counter through ReadEvent Index: trunk/Mars/NEWS =================================================================== --- trunk/Mars/NEWS (revision 9936) +++ trunk/Mars/NEWS (revision 9937) @@ -67,4 +67,6 @@ mapping of the regular patches, resmc/fact-trigger-all.txt the mapping of all patches. + + * Added reading of re-used corsika showers (only supported if MMCS is used) ;star: Index: trunk/Mars/mcorsika/MCorsikaEvtHeader.cc =================================================================== --- trunk/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9936) +++ trunk/Mars/mcorsika/MCorsikaEvtHeader.cc (revision 9937) @@ -32,4 +32,9 @@ // // +// Class Version 3: +// ---------------- +// - UInt_t fNumReuse +// +// ///////////////////////////////////////////////////////////////////////////// #include "MCorsikaEvtHeader.h" @@ -55,5 +60,5 @@ // MCorsikaEvtHeader::MCorsikaEvtHeader(const char *name, const char *title) - : fX(0), fY(0) + : fNumReuse((UInt_t)-1), fX(0), fY(0) { fName = name ? name : "MCorsikaEvtHeader"; @@ -104,5 +109,5 @@ { *fLog << all; - *fLog << "Event Number: " << dec << fEvtNumber << endl; + *fLog << "Event Number: " << dec << fEvtNumber << " (reused=" << fNumReuse << ")" << endl; // *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; *fLog << "Energy: " << fTotalEnergy << "GeV" << endl; @@ -122,7 +127,6 @@ // return FALSE if there is no header anymore, else TRUE // -Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat * fInFormat) +Int_t MCorsikaEvtHeader::ReadEvt(MCorsikaFormat *fInFormat) { - if (!fInFormat->SeekNextBlock("EVTH", 1202)) return kFALSE; @@ -151,17 +155,14 @@ fMomentumZ = -f[8]; // Does this minus make sense?! - // FIXME: Correct for direction of magnetic field! fZd = f[9]; fAz = TMath::Pi()-f[10]; - const Int_t n = TMath::Nint(f[96]); - if (n!=1) - { - *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; - *fLog << warn << "WARNING - This error is replaced by a warning." << endl; - } + // f[96] // Number of reuse of event [max=20] - fX = f[117]; //fX = f[97]; - fY = -f[97]; //fY = f[117]; + // FIXME: Check fNumReuse<20 + fX = f[117 + fNumReuse]; + fY = -f[ 97 + fNumReuse]; + + fWeightedNumPhotons = 0; return !fInFormat->Eof(); @@ -171,7 +172,6 @@ // this member function is for reading the event end block -Bool_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) +Int_t MCorsikaEvtHeader::ReadEvtEnd(MCorsikaFormat * fInFormat) { - if (!fInFormat->SeekNextBlock("EVTE", 1209)) return kFALSE; @@ -192,7 +192,8 @@ } + // FIXME: What's the meaning of this if in case of reusing the event this number + // exists only once? fWeightedNumPhotons = f[1]; return !fInFormat->Eof(); } - Index: trunk/Mars/mcorsika/MCorsikaEvtHeader.h =================================================================== --- trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9936) +++ trunk/Mars/mcorsika/MCorsikaEvtHeader.h (revision 9937) @@ -20,4 +20,5 @@ private: UInt_t fEvtNumber; // Event number + UInt_t fNumReuse; // Counter of the reuse of the same shower // UInt_t fParticleID; // Particle ID (see MMcEvtBasic or CORSIKA manual) Float_t fTotalEnergy; // [GeV] @@ -46,4 +47,5 @@ UInt_t GetEvtNumber() const { return fEvtNumber; } + UInt_t GetNumReuse() const { return fNumReuse; } // UInt_t GetParticleID() const { return fParticleID; } @@ -62,8 +64,11 @@ Double_t GetImpact() const; - Int_t ReadEvt(MCorsikaFormat * fInFormat); // read in event header block - Bool_t ReadEvtEnd(MCorsikaFormat * fInFormat); // read in event end block + void IncNumReuse() { fNumReuse++; } + void ResetNumReuse() { fNumReuse=(UInt_t)-1; } - ClassDef(MCorsikaEvtHeader, 1) // Parameter Conatiner for raw EVENT HEADER + Int_t ReadEvt(MCorsikaFormat *informat); // read in event header block + Int_t ReadEvtEnd(MCorsikaFormat *informat); // read in event end block + + ClassDef(MCorsikaEvtHeader, 3) // Parameter Conatiner for raw EVENT HEADER }; Index: trunk/Mars/mcorsika/MCorsikaRead.cc =================================================================== --- trunk/Mars/mcorsika/MCorsikaRead.cc (revision 9936) +++ trunk/Mars/mcorsika/MCorsikaRead.cc (revision 9937) @@ -319,5 +319,5 @@ } - fNumTotalEvents += fRunHeader->GetNumEvents(); + fNumTotalEvents += fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse(); continue; } @@ -393,28 +393,33 @@ // Read a single event from the stream // -Bool_t MCorsikaRead::ReadEvent() -{ - // - // Read in the next EVENT HEADER (see specification), - // if there is no next event anymore stop eventloop - // - Int_t rc = fEvtHeader->ReadEvt(fInFormat); //read event header block - if (!rc) - return kFALSE; - - rc = fEvent->ReadCorsikaEvt(fInFormat); - - /* - // Check if we are allowed to stop reading in the middle of the data - if (rc==kFALSE && !fForce) - { - *fLog << err; - *fLog << "ERROR - End of file in the middle of the data stream!" << endl; - *fLog << " Set the force-option to force processing." << endl; - return kERROR; - } - */ - - return rc==kTRUE ? fEvtHeader->ReadEvtEnd(fInFormat) : rc; +Int_t MCorsikaRead::ReadEvent() +{ + // Store the position to re-read a single event several times + // FIXME: Does this work with EventIO, too? + if (fEvtHeader->GetNumReuse()>=fRunHeader->GetNumReuse()-1) + { + fInFormat->StorePos(); + fEvtHeader->ResetNumReuse(); + } + else + fInFormat->ResetPos(); + + fEvtHeader->IncNumReuse(); + + // Read the event header + const Int_t rc1 = fEvtHeader->ReadEvt(fInFormat); + if (rc1!=kTRUE) + return rc1; + + // Check if reusage number should be reset and increase it + // Note, that the trick here is that after reset it is set to -1 + + // Read the photons corresponding to the i-th core location + const Int_t rc2 = fEvent->ReadCorsikaEvt(fInFormat, fEvtHeader->GetNumReuse()+1); + if (rc2!=kTRUE) + return rc2; + + // read event end + return fEvtHeader->ReadEvtEnd(fInFormat); } @@ -434,5 +439,5 @@ { // Read a single event from file - const Bool_t rc = ReadEvent(); + const Int_t rc = ReadEvent(); // kFALSE means: end of file (try next one) @@ -465,8 +470,10 @@ Int_t MCorsikaRead::PostProcess() { + const UInt_t n = fNumEvents*fRunHeader->GetNumReuse(); + // // Sanity check for the number of events // - if (fNumEvents==GetNumExecutions()-1 || GetNumExecutions()==0) + if (n==GetNumExecutions()-1 || GetNumExecutions()==0) return kTRUE; @@ -474,5 +481,5 @@ *fLog << "Warning - number of read events (" << GetNumExecutions()-1; *fLog << ") doesn't match number in run header(s) ("; - *fLog << fNumEvents << ")." << endl; + *fLog << fRunHeader->GetNumReuse() << "*" << fNumEvents << "=" << n << ")." << endl; return kTRUE; Index: trunk/Mars/mcorsika/MCorsikaRead.h =================================================================== --- trunk/Mars/mcorsika/MCorsikaRead.h (revision 9936) +++ trunk/Mars/mcorsika/MCorsikaRead.h (revision 9937) @@ -29,6 +29,6 @@ UInt_t fNumTotalEvents; //! total number of events in all files - ifstream *fIn; //! input stream (file to read from) - MCorsikaFormat * fInFormat; //! access to input corsika data + ifstream *fIn; //! input stream (file to read from) + MCorsikaFormat *fInFormat; //! access to input corsika data MParList *fParList; //! tasklist to call ReInit from @@ -42,5 +42,5 @@ Int_t OpenNextFile(Bool_t print=kTRUE); Bool_t CalcNumTotalEvents(); - Bool_t ReadEvent(); + Int_t ReadEvent(); Int_t PreProcess(MParList *pList); Index: trunk/Mars/mcorsika/MCorsikaRunHeader.cc =================================================================== --- trunk/Mars/mcorsika/MCorsikaRunHeader.cc (revision 9936) +++ trunk/Mars/mcorsika/MCorsikaRunHeader.cc (revision 9937) @@ -42,4 +42,8 @@ // + Float_t fAtmosphericCoeffC[5] // + UInt_t fCerenkovFlag +// +// Class Version 3: +// ---------------- +// + UInt_t fNumReuse // //////////////////////////////////////////////////////////////////////////// @@ -164,6 +168,6 @@ fInFormat->UnreadLastHeader(); - const Int_t n = TMath::Nint(g[96]); // Number i of uses of each cherenkov event - if (n!=1) + fNumReuse = TMath::Nint(g[96]); // Number i of uses of each cherenkov event + if (fNumReuse!=1) { *fLog << err << "ERROR - Currently only one impact parameter per event is supported." << endl; @@ -269,5 +273,5 @@ *fLog << "Particle ID: " << MMcEvt::GetParticleName(fParticleID) << endl; if (fNumEvents>0) - *fLog << "Num Events: " << fNumEvents << endl; + *fLog << "Num Events: " << fNumEvents << " (reuse " << fNumReuse << " times)" << endl; *fLog << "Obs Level: "; for (Byte_t i=0; i +! Author(s): Thomas Bretz, 12/2000 ! Author(s): Qi Zhe, 06/2007 ! -! Copyright: CheObs Software Development, 2000-2009 +! Copyright: CheObs Software Development, 2000-2010 ! ! @@ -34,5 +34,9 @@ // Version 1: // ---------- -// - First implementation +// * First implementation +// +// Version 2: +// ---------- +// - fNumPhotons // ///////////////////////////////////////////////////////////////////////////// @@ -57,5 +61,5 @@ MPhotonData::MPhotonData(/*const char *name, const char *title*/) : fPosX(0), fPosY(0), fCosU(0), fCosV(0), fTime(0), fWavelength(0), - fNumPhotons(1), fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED), + /*fNumPhotons(1),*/ fProductionHeight(0), fPrimary(MMcEvtBasic::kUNDEFINED), fTag(-1), fWeight(1) { @@ -82,5 +86,5 @@ MPhotonData &d = static_cast(obj); - d.fNumPhotons = fNumPhotons; +// d.fNumPhotons = fNumPhotons; d.fPosX = fPosX; d.fPosY = fPosY; @@ -185,5 +189,5 @@ fPrimary = MMcEvtBasic::kUNDEFINED; - fNumPhotons = 1; +// fNumPhotons = 1; fTag = -1; fWeight = 1; @@ -201,8 +205,14 @@ // system intpo our own. // -Int_t MPhotonData::FillCorsika(Float_t f[7]) +Int_t MPhotonData::FillCorsika(Float_t f[7], Int_t i) { const UInt_t n = TMath::Nint(f[0]); if (n==0) + // FIXME: Do we need to decode the rest anyway? + return kCONTINUE; + + // Check reuse + const Int_t reuse = (n/1000)%100; // Force this to be 1! + if (reuse!=i) return kCONTINUE; @@ -210,13 +220,4 @@ fWavelength = n%1000; fPrimary = MMcEvtBasic::ParticleId_t(n/100000); - fNumPhotons = (n/1000)%100; // Force this to be 1! - - if (fNumPhotons!=1) - { - // FIXME: Could be done in MPhotonEvent::ReadCorsikaEvent - gLog << err << "ERROR - MPhotonData::FillCorsika: fNumPhotons not 1, but " << fNumPhotons << endl; - gLog << " This is not yet supported." << endl; - return kERROR; - } // x=north, y=west @@ -254,4 +255,10 @@ Int_t MPhotonData::FillEventIO(Float_t f[8]) { + if (TMath::Nint(f[6])!=1) + { + gLog << err << "ERROR - Bunch sizes != 1 are not supported." << endl; + return kFALSE; + } + fPosX = f[1]; // xpos relative to telescope [cm] fPosY = -f[0]; // ypos relative to telescope [cm] @@ -260,7 +267,6 @@ fTime = f[4]; // a relative arival time [ns] fProductionHeight = f[5]; // altitude of emission [cm] - fNumPhotons = TMath::Nint(f[6]); // photons in this bunch +// fNumPhotons = TMath::Nint(f[6]); // photons in this bunch fWavelength = TMath::Nint(f[7]); // so far always zeor = unspec. [nm] - // Now reset all data members which are not in the stream @@ -271,5 +277,5 @@ return kTRUE; } - +/* // -------------------------------------------------------------------------- // @@ -299,4 +305,5 @@ return rc==kTRUE ? !fin.eof() : rc; } +*/ // -------------------------------------------------------------------------- @@ -308,5 +315,6 @@ { gLog << inf << endl; - gLog << "Num Photons: " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl; +// gLog << "Num Photons: " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl; + gLog << "Origin: " << MMcEvtBasic::GetParticleName(fPrimary) << endl; gLog << "Wavelength: " << fWavelength << "nm" << endl; gLog << "Pos X/Y Cos U/V: " << fPosX << "/" << fPosY << " " << fCosU << "/" << fCosV << endl; Index: trunk/Mars/msim/MPhotonData.h =================================================================== --- trunk/Mars/msim/MPhotonData.h (revision 9936) +++ trunk/Mars/msim/MPhotonData.h (revision 9937) @@ -37,15 +37,8 @@ Float_t fTime; // [ns] Time since first interaction or entrance into atmosphere - // 17M UShort_t fWavelength; // [nm] Wavelength - // 19M - UInt_t fNumPhotons; // Number of cherenkov photons ins bunch +// UInt_t fNumPhotons; // Number of cherenkov photons ins bunch Float_t fProductionHeight; // [cm] Height of bunch production MMcEvtBasic::ParticleId_t fPrimary; // Type of emitting particle - // 22M - // gzip - // 25M - // raw - // 32M Int_t fTag; //! A tag for external use @@ -129,12 +122,12 @@ // I/O - Int_t ReadCorsikaEvt(istream &fin); - Int_t ReadRflEvt(istream &fin); + //Int_t ReadCorsikaEvt(istream &fin); + //Int_t ReadRflEvt(istream &fin); - Int_t FillCorsika(Float_t f[7]); + Int_t FillCorsika(Float_t f[7], Int_t i); Int_t FillEventIO(Float_t f[7]); Int_t FillRfl(Float_t f[8]); - ClassDef(MPhotonData, 1) //Container to store a cherenkov photon bunch from a CORSUKA file + ClassDef(MPhotonData, 2) //Container to store a cherenkov photon bunch from a CORSUKA file }; Index: trunk/Mars/msim/MPhotonEvent.cc =================================================================== --- trunk/Mars/msim/MPhotonEvent.cc (revision 9936) +++ trunk/Mars/msim/MPhotonEvent.cc (revision 9937) @@ -449,9 +449,31 @@ } +Double_t MPhotonEvent::GetMeanX() const +{ + const UInt_t n = GetNumPhotons(); + + Double_t mean = 0; + for (UInt_t i=0; i(o)->ReadRflEvt(fin); - const Int_t rc = Add(n).ReadRflEvt(fin); + const Int_t rc = Add(n).ReadRflEvt(fin, i); // Evaluate result from reading event @@ -754,8 +776,8 @@ SetReadyToSave(); - //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; + // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl; return kTRUE; } - +*/ // -------------------------------------------------------------------------- // Index: trunk/Mars/msim/MPhotonEvent.h =================================================================== --- trunk/Mars/msim/MPhotonEvent.h (revision 9936) +++ trunk/Mars/msim/MPhotonEvent.h (revision 9937) @@ -37,4 +37,7 @@ Double_t GetTimeMedianDev() const; + Double_t GetMeanX() const; + Double_t GetMeanY() const; + TClonesArray &GetArray() { return fData; } const TClonesArray &GetArray() const { return fData; } @@ -53,7 +56,7 @@ // I/O - Int_t ReadCorsikaEvt(MCorsikaFormat * fInFormat); - Int_t ReadCorsikaEvt(istream &fin); - Int_t ReadRflEvt(istream &fin); + Int_t ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i); + Int_t ReadCorsikaEvt(istream &fin, Int_t i); + //Int_t ReadRflEvt(istream &fin, Int_t i); // TObject Index: trunk/Mars/msim/MSimMMCS.cc =================================================================== --- trunk/Mars/msim/MSimMMCS.cc (revision 9936) +++ trunk/Mars/msim/MSimMMCS.cc (revision 9937) @@ -131,5 +131,5 @@ // FIXME: Is there a way to write them as LAST entry in the file? - fMcRunHeader->SetNumSimulatedShowers(fRunHeader->GetNumEvents()); + fMcRunHeader->SetNumSimulatedShowers(fRunHeader->GetNumEvents()*fRunHeader->GetNumReuse()); fMcRunHeader->SetCorsikaVersion(TMath::Nint(fRunHeader->GetProgramVersion()*100)); @@ -211,4 +211,5 @@ fMcEvt->SetEvtNumber(fEvtHeader->GetEvtNumber()); + fMcEvt->SetEventReuse(fEvtHeader->GetNumReuse()); fMcEvt->SetPhotElfromShower(0);