Index: /trunk/Mars/msim/MPhotonData.cc
===================================================================
--- /trunk/Mars/msim/MPhotonData.cc	(revision 19338)
+++ /trunk/Mars/msim/MPhotonData.cc	(revision 19339)
@@ -190,5 +190,5 @@
 // contained in the input stream.
 //
-Int_t MPhotonData::FillRfl(Float_t f[8])
+Int_t MPhotonData::FillRfl(const Float_t f[8])
 {
     // Check coordinate system!!!!
@@ -220,23 +220,66 @@
 // system intpo our own.
 //
-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?
+Int_t MPhotonData::FillCorsika(const Float_t f[7], Int_t i)
+{
+    // From the Corsika manual:
+    //
+    // f[0] : n Number of Cherenkov photons in bunch
+    //          (For THIN option multiplied with thinning weight)
+    // f[1] : x [cm]
+    // f[2] : y [cm]
+    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
+    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
+    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
+    // f[6] : h [ch] bunch production height (except MCERFI=3)
+    // f[7] : w weight of bunch [only with THIN option]
+
+    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
+    // f[1] = XCER
+    // f[2] = YXCER
+    // f[3] = UEMIS
+    // f[4] = VEMIS
+    // f[5] = CARTIM
+    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
+    // #if __THIN__
+    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
+    // #endif
+    //
+    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
+    //
+    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
+    //
+    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
+    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
+    //
+
+    const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
+
+    // Empty bunch (this happend at the end of events when
+    // the remaining block is filled with zeroes)
+    if (f[0]==0)
         return kCONTINUE;
 
+    if (f[0]!=1)
+        cout << "WARNING - Bunch size !=1 (" << f[0] << ")" <<endl;
+
+#ifdef __MMCS__
     // Check reuse
     if (i >=0)
-      {
-       const Int_t reuse = (n/1000)%100; // Force this to be 1!
-       if (reuse!=i)
-           return kCONTINUE;
-      }
+    {
+        const Int_t reuse = (n/1000)%100; // Force this to be 1!
+        if (reuse!=i)
+        {
+            cout << "REUSE " << reuse << " " << i << " " << n << endl;
+            return kCONTINUE;
+        }
+    }
 
     // This seems to be special to mmcs
     fWavelength = n%1000;
     fPrimary    = MMcEvtBasic::ParticleId_t(n/100000);
+#else
+    fWavelength = 0;
+    fPrimary    = MMcEvtBasic::kUNDEFINED;
+#endif
 
     // x=north, y=west
@@ -263,4 +306,74 @@
 }
 
+Int_t MPhotonData::FillCorsikaThin(const Float_t f[8], Int_t i)
+{
+    // From the Corsika manual:
+    //
+    // f[0] : n Number of Cherenkov photons in bunch
+    //          (For THIN option multiplied with thinning weight)
+    // f[1] : x [cm]
+    // f[2] : y [cm]
+    // f[3] : u direction cosine (to x axis)  [u = sin(theta)cos(phi)]
+    // f[4] : v direction cosine (to y axis)  [v = sin(theta)sin(phi)]
+    // f[5] : t [ns] time to first interaction or since start of atmosphere (see TSTART)
+    // f[6] : h [ch] bunch production height (except MCERFI=3)
+    // f[7] : w weight of bunch [only with THIN option]
+
+    // f[0] = MCERFI==1/2/3 && !THIN ? PHOTCM : PHOTCM*WTCER/MAX(1e-10, PROBTH)
+    // f[1] = XCER
+    // f[2] = YXCER
+    // f[3] = UEMIS
+    // f[4] = VEMIS
+    // f[5] = CARTIM
+    // f[6] = MCERFI<3 ? ZEMIS : CERDIST
+    // #if __THIN__
+    // f[7] = MCERFI==1 ? WTCER/MAX(1e-10, PROBTH) : (CEFFIC || CERWLEN ? WL*WLFLAG : WLFLAG);
+    // #endif
+    //
+    // WLFLAG = [CEFFIC=-1] [CERWLEN=1] [ELSE=0]
+    //
+    // WL = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*RDM(IRDM))
+    //
+    // WL-MIN = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*0) = WAVLGU
+    // WL-MAX = 1./(1./WAVLGU+(1./WAVLGL-1./WAVLGU)*1) = WAVLGL
+    //
+
+    const UInt_t n /*fWeight*/ = TMath::Nint(f[0]);
+
+    // Empty bunch (this happend at the end of events when
+    // the remaining block is filled with zeroes)
+    if (f[0]==0)
+        return kCONTINUE;
+
+    if (f[0]!=1)
+        cout << "WARNING - Bunch size !=1 (" << f[0] << ")" <<endl;
+
+
+    // x=north, y=west
+    //fPosX = f[1];  // [cm]
+    //fPosY = f[2];  // [cm]
+    //fCosU = f[3];  // cos to x
+    //fCosV = f[4];  // cos to y
+    // x=west, y=south
+    fPosX =  f[2];  // [cm]
+    fPosY = -f[1];  // [cm]
+
+    fCosU =  f[4];  // cos to x
+    fCosV = -f[3];  // cos to y
+
+    fTime =  f[5];  // [ns]
+
+    fProductionHeight = f[6]; // [cm]
+
+    fWavelength = TMath::Abs(f[7]);
+
+    // Now reset all data members which are not in the stream
+    fPrimary = MMcEvtBasic::kUNDEFINED;
+    fTag     = -1;
+    fWeight  =  1;
+
+    return kTRUE;
+}
+
 // --------------------------------------------------------------------------
 //
@@ -272,5 +385,5 @@
 // system into our own.
 //
-Int_t MPhotonData::FillEventIO(Short_t f[8])
+Int_t MPhotonData::FillEventIO(const Short_t f[8])
 {
     // From 5.5 compact_bunch:
@@ -286,5 +399,5 @@
     fTime             =  f[4]/10.;              // a relative arival time [ns]
     fProductionHeight =  pow(10, f[5]/1000.);   // altitude of emission a.s.l. [cm]
-    fWavelength       =  f[7];                  // wavelength [nm]: 0 undetermined, <0 already in p.e.
+    fWavelength       =  TMath::Abs(f[7]);      // wavelength [nm]: 0 undetermined, <0 already in p.e.
 
     // Now reset all data members which are not in the stream
@@ -305,5 +418,5 @@
 // system into our own.
 //
-Int_t MPhotonData::FillEventIO(Float_t f[8])
+Int_t MPhotonData::FillEventIO(const Float_t f[8])
 {
     // photons in this bunch
@@ -366,6 +479,7 @@
 {
     gLog << inf << endl;
-//    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
-    gLog << "Origin:           " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
+    //    gLog << "Num Photons:      " << fNumPhotons << " from " << MMcEvtBasic::GetParticleName(fPrimary) << endl;
+    if (fPrimary!=MMcEvtBasic::kUNDEFINED)
+        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 19338)
+++ /trunk/Mars/msim/MPhotonData.h	(revision 19339)
@@ -134,8 +134,9 @@
     //Int_t ReadRflEvt(istream &fin);
 
-    Int_t FillCorsika(Float_t f[7], Int_t i);
-    Int_t FillEventIO(Short_t f[8]);
-    Int_t FillEventIO(Float_t f[8]);
-    Int_t FillRfl(Float_t f[8]);
+    Int_t FillCorsika(const Float_t f[7], Int_t i);
+    Int_t FillCorsikaThin(const Float_t f[8], Int_t i);
+    Int_t FillEventIO(const Short_t f[8]);
+    Int_t FillEventIO(const Float_t f[8]);
+    Int_t FillRfl(const Float_t f[8]);
 
     ClassDef(MPhotonData, 2) //Container to store a cherenkov photon bunch from a CORSUKA file