Index: trunk/MagicSoft/Mars/msimcamera/MSimCamera.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimCamera.h (revision 9521) +++ trunk/MagicSoft/Mars/msimcamera/MSimCamera.h (revision 9525) @@ -28,5 +28,5 @@ 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? Index: trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (revision 9521) +++ trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc (revision 9525) @@ -52,6 +52,7 @@ // 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,4 +69,28 @@ // 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,5 +143,6 @@ 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,62 +195,59 @@ }*/ - 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,4 +260,130 @@ // 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 (maxGetWavelengthMax()) + { + *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,5 +423,5 @@ { // 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,5 +454,5 @@ // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant) FIXME: Reset? - if (fRunHeader) + if (fSimulateWavelength) { const Float_t wmin = fRunHeader->GetWavelengthMin(); @@ -353,4 +502,16 @@ } + 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; } Index: trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h =================================================================== --- trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h (revision 9521) +++ trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.h (revision 9525) @@ -33,4 +33,5 @@ TString fNameGeomCam; + TString fFileNameNSB; // MTask