Changeset 19339

Timestamp:

10/30/18 17:31:06 (6 years ago)

Author:

tbretz

Message:

Added the 8-float format as FillEvtCorsikaThin, calculate absolute wavelength (right now, I don't know what to do with the information from the negative WL, print the origin only if it is defined, added some docu from the corsika docu, removed MMCS specific stuff by a pre-compiler directive -- do we still have MMCS data?

Location:

trunk/Mars/msim

Files:

• MPhotonData.cc (modified) (7 diffs)
• MPhotonData.h (modified) (1 diff)

trunk/Mars/msim/MPhotonData.cc

@@ -190 +190 @@
 // 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 +220 @@
 // 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 +306 @@
 }
 
+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 +385 @@
 // 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 +399 @@
     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 +418 @@
 // 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 +479 @@
 {
     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;

trunk/Mars/msim/MPhotonData.h

@@ -134 +134 @@
     //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