Changeset 9913 for trunk/Mars/msimcamera/MSimRandomPhotons.cc

Timestamp:

08/31/10 10:26:32 (14 years ago)

Author:

tbretz

Message:

Improved the range checks in MSimRandomPhotons.

File:

• trunk/Mars/msimcamera/MSimRandomPhotons.cc (modified) (9 diffs)

trunk/Mars/msimcamera/MSimRandomPhotons.cc

@@ -150 +150 @@
 }
 
+Bool_t MSimRandomPhotons::CheckWavelengthRange(const MParSpline &sp, const char *txt) const
+{
+    const Float_t min = sp.GetXmin();
+    const Float_t max = sp.GetXmax();
+
+    if (min>fRunHeader->GetWavelengthMin())
+    {
+        *fLog << err << "ERROR - Minimum wavelength (" << min << "nm)";
+        *fLog << " defined for " << txt;
+        *fLog << " exceeds minimum wavelength simulated (";
+        *fLog << fRunHeader->GetWavelengthMin() << "nm)." << endl;
+        return kFALSE;
+    }
+    if (max<fRunHeader->GetWavelengthMax())
+    {
+        *fLog << err << "ERROR - Maximum wavelength (" << max << "nm)";
+        *fLog << " defined for " << txt;
+        *fLog << " undershoots maximum wavelength simulated (";
+        *fLog << fRunHeader->GetWavelengthMax() << "nm)." << endl;
+        return kFALSE;
+    }
+
+    return kTRUE;
+}
+
 // --------------------------------------------------------------------------
 //
@@ -196 +221 @@
 
     fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
-    if (fSimulateWavelength && !fRunHeader)
+    if (!fRunHeader)
     {
         *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
@@ -214 +239 @@
     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;
+    // Check if all splines are defined in the relevant range
+    if (!CheckWavelengthRange(*s1, "PhotonDetectionEfficiency [MParSpline]"))
+        return kFALSE;
+    if (!CheckWavelengthRange(*s2, "ConesTransmission [MParSpline]"))
+        return kFALSE;
+    if (!CheckWavelengthRange(*s3, "MirrorReflectivity [MParSpline]"))
+        return kFALSE;
+
+    const Float_t wmin = fRunHeader->GetWavelengthMin();
+    const Float_t wmax = fRunHeader->GetWavelengthMax();
+
+    const Int_t min = TMath::FloorNint(wmin);
+    const Int_t max = TMath::CeilNint(wmax);
+
+    // Multiply all relevant efficiencies to get the total transmission
+    MParSpline eff;
+    eff.SetFunction("1", max-min, min, max);
+
+    eff.Multiply(*s1->GetSpline());
+    eff.Multiply(*s2->GetSpline());
+    eff.Multiply(*s3->GetSpline());
+
+    // Effectively transmitted wavelength band in the simulated range
+    const Double_t nm = eff.GetSpline()->Integral();
 
     // Angular acceptance of the cones
@@ -246 +288 @@
 
     if (fFileNameNSB.IsNull())
-    {
-        delete eff;
         return kTRUE;
-    }
 
     // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
@@ -277 +316 @@
         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;
-        }
-    }
+    if (!CheckWavelengthRange(flux, TString("night sky background flux from ")+fFileNameNSB))
+        return kFALSE;
 
     MParSpline nsb;
@@ -308 +329 @@
 
     // Multiply with the total transmission
-    nsb.Multiply(*eff->GetSpline());
-
-    // Check if the photon flux is zero at both ends
+    nsb.Multiply(*eff.GetSpline());
+
+    // Check if the photon flux is zero at both ends of the NSB
     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;
+        *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at ";
+        *fLog << wmin << "nm is not zero, but " << nsb.GetSpline()->Eval(wmin) << "... abort." << endl;
     }
     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;
+        *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at ";
+        *fLog << wmax << "nm is not zero, but " << nsb.GetSpline()->Eval(wmax) << "... abort." << endl;
+    }
+
+    // Check if the photon flux is zero at both ends of the simulated region
+    if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMin())>1e-5)
+    {
+        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at minimum simulated wavelength ";
+        *fLog << fRunHeader->GetWavelengthMin() << "nm is not zero... abort." << endl;
+        return kFALSE;
+    }
+    if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMax())>1e-5)
+    {
+        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at maximum simulated wavelength ";
+        *fLog << fRunHeader->GetWavelengthMax() << "nm is not zero... abort." << endl;
+        return kFALSE;
     }
 
@@ -336 +363 @@
     //const Double_t sr = s5.GetSpline()->IntegralSolidAngle(); // sr
     // Absolute reflector area
-    const Double_t Ar = r->GetA()/1e4;                      // m^2
+    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
+    const Double_t A0 = (*fGeom)[0].GetA()*1e-6;                // m^2
 
     // Rate * m^2 * Solid Angle
@@ -360 +387 @@
 
     // 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 << 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;
@@ -385 +412 @@
     }
 
-    delete eff;
-
     return kTRUE;
 }