Changeset 9937 for trunk/Mars/msim

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/msim

Files:

• MPhotonData.cc (modified) (12 diffs)
• MPhotonData.h (modified) (2 diffs)
• MPhotonEvent.cc (modified) (7 diffs)
• MPhotonEvent.h (modified) (2 diffs)
• MSimMMCS.cc (modified) (2 diffs)

trunk/Mars/msim/MPhotonData.cc

@@ -16 +16 @@
 !
 !
-!   Author(s): Thomas Bretz, 12/2000 <mailto:tbretz@astro.uni-wuerzburg.de>
+!   Author(s): Thomas Bretz, 12/2000 <mailto:thomas.bretz@epfl.ch>
 !   Author(s): Qi Zhe,       06/2007 <mailto:qizhe@astro.uni-wuerzburg.de>
 !
-!   Copyright: CheObs Software Development, 2000-2009
+!   Copyright: CheObs Software Development, 2000-2010
 !
 !
@@ -34 +34 @@
 //   Version 1:
 //   ----------
-//    - First implementation
+//    * First implementation
+//
+//   Version 2:
+//   ----------
+//    - fNumPhotons
 //
 /////////////////////////////////////////////////////////////////////////////
@@ -57 +61 @@
 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 +86 @@
     MPhotonData &d = static_cast<MPhotonData&>(obj);
 
-    d.fNumPhotons       = fNumPhotons;
+//    d.fNumPhotons       = fNumPhotons;
     d.fPosX             = fPosX;
     d.fPosY             = fPosY;
@@ -185 +189 @@
 
     fPrimary    = MMcEvtBasic::kUNDEFINED;
-    fNumPhotons =  1;
+//    fNumPhotons =  1;
     fTag        = -1;
     fWeight     =  1;
@@ -201 +205 @@
 // 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 +220 @@
     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 +255 @@
 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 +267 @@
     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 +277 @@
     return kTRUE;
 }
-
+/*
 // --------------------------------------------------------------------------
 //
@@ -299 +305 @@
     return rc==kTRUE ? !fin.eof() : rc;
 }
+*/
 
 // --------------------------------------------------------------------------
@@ -308 +315 @@
 {
     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;

trunk/Mars/msim/MPhotonData.h

@@ -37 +37 @@
 
     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 +122 @@
 
     // 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
 };
 

trunk/Mars/msim/MPhotonEvent.cc

@@ -449 +449 @@
 }
 
+Double_t MPhotonEvent::GetMeanX() const
+{
+    const UInt_t n = GetNumPhotons();
+
+    Double_t mean = 0;
+    for (UInt_t i=0; i<n; i++)
+        mean += operator[](i).GetPosX();
+
+    return mean / n;
+}
+
+Double_t MPhotonEvent::GetMeanY() const
+{
+    const UInt_t n = GetNumPhotons();
+
+    Double_t mean = 0;
+    for (UInt_t i=0; i<n; i++)
+        mean += operator[](i).GetPosY();
+
+    return mean / n;
+}
+
 // --------------------------------------------------------------------------
 //
 // Read the Event section from the file
 //
-Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat * fInFormat)
+Int_t MPhotonEvent::ReadCorsikaEvt(MCorsikaFormat *fInFormat, Int_t i)
 {
     Int_t n = 0;
@@ -509 +531 @@
             {
 
-            const Int_t rc = Add(n).FillCorsika(buffer);
+            const Int_t rc = Add(n).FillCorsika(buffer, i);
             switch (rc)
             {
@@ -577 +599 @@
 }
 
-Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin)
+Int_t MPhotonEvent::ReadCorsikaEvt(istream &fin, Int_t i)
 {
     Int_t n = 0;
@@ -635 +657 @@
             // Get/Add the n-th entry from the array and
             // fill it with the current 7 floats
-            const Int_t rc = Add(n).FillCorsika(ptr);
+            const Int_t rc = Add(n).FillCorsika(ptr, i);
             ptr += 7;
 
@@ -707 +729 @@
 
 // --------------------------------------------------------------------------
-//
-Int_t MPhotonEvent::ReadRflEvt(std::istream &fin)
+/*
+Int_t MPhotonEvent::ReadRflEvt(std::istream &fin, Int_t i)
 {
     Int_t n = 0;
@@ -735 +757 @@
         // Now we read a single cherenkov bunch
         //const Int_t rc = static_cast<MPhotonData*>(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 +776 @@
     SetReadyToSave();
 
-    //*fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
+    // *fLog << all << "Number of photon bunches: " << fData.GetEntriesFast() << endl;
     return kTRUE;
 }
-
+*/
 // --------------------------------------------------------------------------
 //

trunk/Mars/msim/MPhotonEvent.h

@@ -37 +37 @@
     Double_t GetTimeMedianDev() const;
 
+    Double_t GetMeanX() const;
+    Double_t GetMeanY() const;
+
     TClonesArray &GetArray() { return fData; }
     const TClonesArray &GetArray() const { return fData; }
@@ -53 +56 @@
 
     // 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

trunk/Mars/msim/MSimMMCS.cc

@@ -131 +131 @@
 
     // 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 +211 @@
 
     fMcEvt->SetEvtNumber(fEvtHeader->GetEvtNumber());
+    fMcEvt->SetEventReuse(fEvtHeader->GetNumReuse());
     fMcEvt->SetPhotElfromShower(0);