Changeset 9525 for trunk/MagicSoft/Mars/msimcamera

Timestamp:

12/10/09 10:08:55 (15 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/MagicSoft/Mars/msimcamera

Files:

• MSimCamera.h (modified) (1 diff)
• MSimRandomPhotons.cc (modified) (8 diffs)
• MSimRandomPhotons.h (modified) (1 diff)

trunk/MagicSoft/Mars/msimcamera/MSimCamera.h

@@ -28 +28 @@
     MMcEvt            *fMcEvt;           //! For information stored in MMcEvt
 
-    const MSpline3    *fSpline;          // Pusle Shape
+    const MSpline3    *fSpline;          // Pulse Shape
 
     Bool_t fBaselineGain;  // Should the gain be applied to baseline and electronic noise?

trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc

@@ -52 +52 @@
 //   R = R0 * Ai * f
 //
-// R0 is the night sky background rate as given in Eckart's paper. Ai the
-// area of the cones acceptance window, f is given as:
+// R0 is the night sky background rate as given in Eckart's paper (divided
+// by the wavelength window). Ai the area of the cones acceptance window,
+// f is given as:
 //
 //   f = nm * Min(Ar, sr*d^2)
@@ -68 +69 @@
 //   sr is the effective solid angle corresponding to the integral of the
 //   cones angular acceptance
+//
+// Alternatively, the night-sky background rate can be calculated from
+// a spectrum as given in Fig. 1 (but versus Nanometers) in
+//
+//   Chris R. Benn & Sara L. Ellison La Palma Night-Sky Brightness
+//
+// After proper conversion of the units, the rate of the pixel 0
+// is then calculated by
+//
+//     rate = f * nsb
+//
+// With nsb
+//
+//   nsb = Integral(nsb spectrum * combines efficiencies)
+//
+// and f can be either
+//
+//   Eff. angular acceptance Cones (e.g. 20deg) * Cone-Area (mm^2)
+//   f = sr * A0
+//
+// or
+//
+//   Mirror-Area * Field of view of cones (deg^2)
+//   f = Ar * A0;
 //
 //
@@ -118 +143 @@
 MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title)
     : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0),
-    fRates(0), fSimulateWavelength(kFALSE),  fNameGeomCam("MGeomCam")
+    fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam"),
+    fFileNameNSB("resmc/night-sky-la-palma.txt")
 {
     fName  = name  ? name  : "MSimRandomPhotons";
@@ -169 +195 @@
     }*/
 
-    fRunHeader = 0;
-    if (fSimulateWavelength)
-    {
-        fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
-        if (!fRunHeader)
+    fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
+    if (fSimulateWavelength && !fRunHeader)
+    {
+        *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
+    if (!r)
+    {
+        *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
+        return kFALSE;
+    }
+
+    const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
+    const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesTransmission",         "MParSpline");
+    const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity",        "MParSpline");
+    const MParSpline *s4 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
+
+    // Multiply all relevant efficiencies to get the total tarnsmission
+    MParSpline *eff = (MParSpline*)s1->Clone();
+    eff->Multiply(*s2->GetSpline());
+    eff->Multiply(*s3->GetSpline());
+
+    // Effectively transmitted wavelength band
+    const Double_t nm = eff && eff->GetSpline() ? eff->GetSpline()->Integral() : 1;
+
+    // Angular acceptance of the cones
+    const Double_t sr = s4 && s4->GetSpline() ? s4->GetSpline()->IntegralSolidAngle() : 1;
+
+    {
+        const Double_t d2   = fGeom->GetCameraDist()*fGeom->GetCameraDist();
+        const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();
+        const Double_t f1   = TMath::Min(r->GetA()/1e4, sr*d2) * conv*conv;
+
+        // Rate in GHz / mm^2
+        fScale = fFreqNSB * nm * f1; // [GHz/mm^2] efficiency * m^2 *rad^2 *mm^2
+
+        const Double_t freq0 = fScale*(*fGeom)[0].GetA()*1000;
+
+        *fLog << inf << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", freq0) << "MHz" << endl;
+
+        // FIXME: Scale with the number of pixels
+        if (freq0>1000)
         {
-            *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
+            *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;
             return kFALSE;
         }
     }
 
-    MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
-    if (!r)
-    {
-        *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
-        return kFALSE;
-    }
-
-    *fLog << warn << "FIXME: A check is needed that the PDE is 0 at both ends!" << endl;
-
-    const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
-    const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
-    const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity",        "MParSpline");
-    const MParSpline *s5 = (MParSpline*)pList->FindObject("ConesTransmission",         "MParSpline");
-
-    const Double_t d2  = fGeom->GetCameraDist()*fGeom->GetCameraDist();
-//    const Double_t pde = s1 && s1->GetSpline() ? s1->GetSpline()->Integral() : 1;
-    const Double_t sr  = s2 && s2->GetSpline() ? s2->GetSpline()->IntegralSolidAngle() : 1;
-//    const Double_t mir = s3 && s3->GetSpline() ? s3->GetSpline()->Integral() : 1;
-    const Double_t Ar  = r->GetA()/1e4;
-
-    // Conversion factor to convert pixel area to steradians (because it
-    // is a rather small area we can assume it is flat)
-    const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();
-
-    // Multiply all relevant efficiencies
-    MParSpline *s4 = (MParSpline*)s1->Clone();
-    s4->Multiply(*s3->GetSpline());
-    s4->Multiply(*s5->GetSpline());
-
-    const Double_t nm = s4 && s4->GetSpline() ? s4->GetSpline()->Integral() : 1;
-
-    delete s4;
-
-    // /100 to convert the pixel area from mm^2 to cm^2
-    fScale =  nm * TMath::Min(Ar, sr*d2) * conv*conv;
-
-    *fLog << inf;
-    *fLog << "Effective cone acceptance:      " << Form("%.2f", sr*d2) << "m^2" << endl;
-    *fLog << "Reflector area:                 " << Form("%.2f", Ar) << "m^2" << endl;
-    *fLog << "Resulting eff. collection area: " << Form("%.1f", TMath::Min(Ar, sr*d2)) << "m^2" << endl;
-//    *fLog << "Eff. wavelength band (PDE):     " << Form("%.1f", pde) << "nm" << endl;
-//    *fLog << "Eff. wavelength band (Mirror):  " << Form("%.1f", mir) << "nm" << endl;
-    *fLog << "Eff. wavelength band (MIR+Cone+PDE): " << Form("%.1f", nm) << "nm" << endl;
-    *fLog << "Pixel area of " << fNameGeomCam << "[0]: " << Form("%.2e", (*fGeom)[0].GetA()*conv*conv) << "sr" << endl;
-    //*fLog << "Effective angular acceptance:     " << sr << "sr" << endl;
-    //*fLog << "Resulting NSB frequency:           " << fFreqNSB*nm*Ar*1000 << "MHz/sr" << endl;
-    *fLog << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", fFreqNSB*(*fGeom)[0].GetA()*fScale*1000) << "MHz" << endl;
+    if (fFileNameNSB.IsNull())
+    {
+        delete eff;
+        return kTRUE;
+    }
 
     // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
@@ -237 +260 @@
     // ampl = MMcFadcHeader->GetAmplitud()
     // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio)
+
+    // Conversion of the y-axis
+    // ------------------------
+    // Double_t ff = 1;                               // myJy / arcsec^2 per nm
+    // ff *= 1e-6;                                    // Jy   / arcsec^2 per nm
+    // ff *= 3600*3600;                               // Jy   / deg^2
+    // ff *= 1./TMath::DegToRad()/TMath::DegToRad();  // Jy/sr = 1e-26J/s/m^2/Hz/sr
+    // ff *= 1e-26;                                   // J/s/m^2/Hz/sr   per nm
+
+    const Double_t arcsec2rad = TMath::DegToRad()/3600.;
+    const Double_t f = 1e-32 / (arcsec2rad*arcsec2rad);
+
+    // Read night sky background flux from file
+    MParSpline flux;
+    if (!flux.ReadFile(fFileNameNSB))
+        return kFALSE;
+
+    const Int_t min = TMath::FloorNint(flux.GetXmin());
+    const Int_t max = TMath::CeilNint( flux.GetXmax());
+
+    if (fRunHeader)
+    {
+        if (min>fRunHeader->GetWavelengthMin())
+        {
+            *fLog << warn << "WARNING - Minimum wavelength of night sky background flux (";
+            *fLog << min << "nm) from " << fFileNameNSB;
+            *fLog << " exceeds minimum wavelength simulated ";
+            *fLog << fRunHeader->GetWavelengthMin() << "nm." << endl;
+        }
+        if (max<fRunHeader->GetWavelengthMax())
+        {
+            *fLog << warn << "WARNING - Maximum wavelength of night sky background flux (";
+            *fLog << max << "nm) from " << fFileNameNSB;
+            *fLog << " undershoots maximum wavelength simulated ";
+            *fLog << fRunHeader->GetWavelengthMax() << "nm." << endl;
+        }
+    }
+
+    MParSpline nsb;
+
+    // Normalization to our units,
+    // conversion from energy flux to photon flux
+    nsb.SetFunction(Form("%.12e/(x*TMath::H())", f), max-min, min, max);
+
+    // multiply night sky background flux with normalization
+    nsb.Multiply(*flux.GetSpline());
+
+    // Multiply with the total transmission
+    nsb.Multiply(*eff->GetSpline());
+
+    // Check if the photon flux is zero at both ends
+    if (nsb.GetSpline()->Eval(min)>1e-5)
+    {
+        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << min << "nm is not zero... abort." << endl;
+        return kFALSE;
+    }
+    if (nsb.GetSpline()->Eval(max)>1e-5)
+    {
+        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << max << "nm is not zero... abort." << endl;
+        return kFALSE;
+    }
+
+    if (fRunHeader)
+    {
+        if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMin())>1e-5)
+            *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at " << fRunHeader->GetWavelengthMin() << "nm is not zero... abort." << endl;
+        if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMax())>1e-5)
+            *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at " << fRunHeader->GetWavelengthMax() << "nm is not zero... abort." << endl;
+    }
+
+    // Conversion from m to radians
+    const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad()*1e3;
+
+    // Angular acceptance of the cones
+    //const Double_t sr = s5.GetSpline()->IntegralSolidAngle(); // sr
+    // Absolute reflector area
+    const Double_t Ar = r->GetA()/1e4;                      // m^2
+    // Size of the cone's entrance window
+    const Double_t A0 = (*fGeom)[0].GetA()*1e-6;                  // m^2
+
+    // Rate * m^2 * Solid Angle
+    // -------------------------
+
+    // Angular acceptance Cones (e.g. 20deg) * Cone-Area
+    const Double_t f1 = A0 * sr;                // m^2 sr
+
+    // Mirror-Area * Field of view of cones (e.g. 0.1deg)
+    const Double_t f2 = Ar * A0*conv*conv;      // m^2 sr
+
+    // FIXME: Calculate the reflectivity of the bottom by replacing
+    // MirrorReflectivity by bottom reflectivity and reflect
+    // and use it to reflect the difference between f1 and f2
+    // if any.
+
+    // Total NSB rate in MHz per m^2 and sr
+    const Double_t rate = nsb.GetSpline()->Integral() * 1e-6;
+
+    *fLog << inf;
+
+    // Resulting rates as if Razmick's constant had been used
+    // *fLog << 1.75e6/(600-300) * f1 * eff->GetSpline()->Integral() << " MHz" << endl;
+    // *fLog << 1.75e6/(600-300) * f2 * eff->GetSpline()->Integral() << " MHz" << endl;
+
+    *fLog << "Conversion factor Fnu:      " << f  << endl;
+    *fLog << "Total reflective area:      " << Form("%.2f", Ar) << " m" << UTF8::kSquare << endl;
+    *fLog << "Acceptance area of cone 0:  " << Form("%.2f", A0*1e6) << " mm" << UTF8::kSquare << " = ";
+    *fLog << A0*conv*conv << " sr" << endl;
+    *fLog << "Cones angular acceptance:   " << sr << " sr" << endl;
+    *fLog << "ConeArea*MirrorAngle (f1):  " << f1 << " m^2 sr" << endl;
+    *fLog << "MirrorArea*ConeAngle (f2):  " << f2 << " m^2 sr" << endl;
+    *fLog << "Effective. transmission:    " << Form("%.1f", nm) << " nm" << endl;
+    *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f1): " << Form("%.2f", rate * f1) << " MHz" << endl;
+    *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f2): " << Form("%.2f", rate * f2) << " MHz" << endl;
+    *fLog << "Using f1." << endl;
+
+    // Scale the rate per mm^2 and to GHz
+    fScale = rate * f1 / (*fGeom)[0].GetA() / 1000;
+
+    // FIXME: Scale with the number of pixels
+    if (rate*f1>1000)
+    {
+        *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;
+        return kFALSE;
+    }
+
+    delete eff;
 
     return kTRUE;
@@ -274 +423 @@
     {
         // Scale the rate with the pixel size.
-        const Double_t rate = fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()*fScale;
+        const Double_t rate = fFreqFixed + fScale*(*fGeom)[idx].GetA();
 
         (*fRates)[idx].SetPedestal(rate);
@@ -305 +454 @@
             // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant)  FIXME: Reset?
 
-            if (fRunHeader)
+            if (fSimulateWavelength)
             {
                 const Float_t wmin = fRunHeader->GetWavelengthMin();
@@ -353 +502 @@
     }
 
+    if (IsEnvDefined(env, prefix, "FileNameNSB", print))
+    {
+        rc = kTRUE;
+        fFileNameNSB = GetEnvValue(env, prefix, "FileNameNSB", fFileNameNSB);
+    }
+
+    if (IsEnvDefined(env, prefix, "SimulateCherenkovSpectrum", print))
+    {
+        rc = kTRUE;
+        fSimulateWavelength = GetEnvValue(env, prefix, "SimulateCherenkovSpectrum", fSimulateWavelength);
+    }
+
     return rc;
 }

trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h

@@ -33 +33 @@
 
     TString fNameGeomCam;
+    TString fFileNameNSB;
 
     // MTask