/* ======================================================================== *\ ! ! * ! * This file is part of CheObs, the Modular Analysis and Reconstruction ! * Software. It is distributed to you in the hope that it can be a useful ! * and timesaving tool in analysing Data of imaging Cerenkov telescopes. ! * It is distributed WITHOUT ANY WARRANTY. ! * ! * Permission to use, copy, modify and distribute this software and its ! * documentation for any purpose is hereby granted without fee, ! * provided that the above copyright notice appears in all copies and ! * that both that copyright notice and this permission notice appear ! * in supporting documentation. It is provided "as is" without express ! * or implied warranty. ! * ! ! ! Author(s): Thomas Bretz, 1/2009 ! ! Copyright: CheObs Software Development, 2000-2009 ! ! \* ======================================================================== */ ////////////////////////////////////////////////////////////////////////////// // // MSimRandomPhotons // // Simulate poissonian photons. Since the distribution of the arrival time // differences of these photons is an exonential we can simulate them // using exponentially distributed time differences between two consecutive // photons. // // FIXME: We should add the wavelength distribution. // // The artificial night sky background rate is calculated as follows: // // * The photon detection efficiency vs. wavelength of the detector is obtained // from "PhotonDetectionEfficiency" of type "MParSpline" // // * The angular acceptance of the light collectors is obtained // from "ConesAngularAcceptance" of type "MParSpline" // // * The spectral acceptance of the light collectors is obtained // from "ConesTransmission" of type "MParSpline" // // * The reflectivity of the mirrors vs wavelength is obtained // from "MirrorReflectivity" of type "MParSpline" // // The rate is then calculated as // // R = R0 * Ai * f // // 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) // // with // // nm being the integral of the product of the mirror reflectivity, the cone // transmission and the photon detection efficiency. // // d the distance of the focal plane to the mirror // // Ar is the total reflective area of the reflector // // 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; // // // Input Containers: // fNameGeomCam [MGeomCam] // MPhotonEvent // MPhotonStatistics // MCorsikaEvtHeader // [MCorsikaRunHeader] // // Output Containers: // MPhotonEvent // AccidentalPhotonRate [MPedestalCam] // ////////////////////////////////////////////////////////////////////////////// #include "MSimRandomPhotons.h" #include #include "MMath.h" // RndmExp #include "MLog.h" #include "MLogManip.h" #include "MParList.h" #include "MGeomCam.h" #include "MGeom.h" #include "MPhotonEvent.h" #include "MPhotonData.h" #include "MPedestalCam.h" #include "MPedestalPix.h" #include "MCorsikaRunHeader.h" #include "MSpline3.h" #include "MParSpline.h" #include "MReflector.h" ClassImp(MSimRandomPhotons); using namespace std; // -------------------------------------------------------------------------- // // Default Constructor. // MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title) : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0), fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam"), fFileNameNSB("resmc/night-sky-la-palma.txt") { fName = name ? name : "MSimRandomPhotons"; fTitle = title ? title : "Simulate possonian photons (like NSB or dark current)"; } 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 (maxGetWavelengthMax()) { *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; } // -------------------------------------------------------------------------- // // Check for the necessary containers // Int_t MSimRandomPhotons::PreProcess(MParList *pList) { fGeom = (MGeomCam*)pList->FindObject(fNameGeomCam, "MGeomCam"); if (!fGeom) { *fLog << inf << fNameGeomCam << " [MGeomCam] not found..." << endl; fGeom = (MGeomCam*)pList->FindObject("MGeomCam"); if (!fGeom) { *fLog << err << "MGeomCam not found... aborting." << endl; return kFALSE; } } fEvt = (MPhotonEvent*)pList->FindObject("MPhotonEvent"); if (!fEvt) { *fLog << err << "MPhotonEvent not found... aborting." << endl; return kFALSE; } fStat = (MPhotonStatistics*)pList->FindObject("MPhotonStatistics"); if (!fStat) { *fLog << err << "MPhotonStatistics not found... aborting." << endl; return kFALSE; } fRates = (MPedestalCam*)pList->FindCreateObj("MPedestalCam", "AccidentalPhotonRates"); if (!fRates) return kFALSE; /* fEvtHeader = (MCorsikaEvtHeader*)pList->FindObject("MCorsikaEvtHeader"); if (!fEvtHeader) { *fLog << err << "MCorsikaEvtHeader not found... aborting." << endl; return kFALSE; }*/ fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader"); if (!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"); // 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 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 << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl; return kFALSE; } } if (fFileNameNSB.IsNull()) return kTRUE; // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader"); // Set NumPheFromDNSB // # Number of photons from the diffuse NSB (nphe / ns 0.1*0.1 deg^2 239 m^2) and // nsb_mean 0.20 // Magic pixel: 0.00885361 deg // dnsbpix = 0.2*50/15 // 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; if (!CheckWavelengthRange(flux, TString("night sky background flux from ")+fFileNameNSB)) return kFALSE; 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 of the NSB if (eff.GetSpline()->Eval(min)>1e-5) { *fLog << warn << "WARNING - Total transmission efficiency at detector at "; *fLog << min << "nm is not zero, but " << nsb.GetSpline()->Eval(min) << "... abort." << endl; } if (eff.GetSpline()->Eval(max)>1e-5) { *fLog << warn << "WARNING - Total transmission efficiency at detector at "; *fLog << max << "nm is not zero, but " << nsb.GetSpline()->Eval(max) << "... abort." << endl; } // Check if the photon flux is zero at both ends of the simulated region if (eff.GetSpline()->Eval(wmin)>1e-5) { *fLog << err << "ERROR - Total transmission efficiency at detector at minimum simulated wavelength "; *fLog << wmin << "nm is not zero... abort." << endl; *fLog << " PhotonDetectionEfficency: " << s1->GetSpline()->Eval(wmin) << endl; *fLog << " ConeTransmission: " << s2->GetSpline()->Eval(wmin) << endl; *fLog << " MirrorReflectivity: " << s3->GetSpline()->Eval(wmin) << endl; *fLog << " TotalEfficiency: " << eff.GetSpline()->Eval(wmin) << endl; return kFALSE; } if (eff.GetSpline()->Eval(wmax)>1e-5) { *fLog << err << "ERROR - Total transmission efficiency at detector at maximum simulated wavelength "; *fLog << wmax << "nm is not zero... abort." << endl; *fLog << " PhotonDetectionEfficency: " << s1->GetSpline()->Eval(wmax) << endl; *fLog << " ConeTransmission: " << s2->GetSpline()->Eval(wmax) << endl; *fLog << " MirrorReflectivity: " << s3->GetSpline()->Eval(wmax) << endl; *fLog << " TotalEfficiency: " << eff.GetSpline()->Eval(wmax) << endl; return kFALSE; } // 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*ConeSolidAngle (f1): " << f1 << " m^2 sr" << endl; *fLog << "MirrorArea*ConeSkyAngle (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 f2." << endl; // Scale the rate per mm^2 and to GHz fScale = rate * f2 / (*fGeom)[0].GetA() / 1000; // FIXME: Scale with the number of pixels if (rate*f2>1000) { *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl; return kFALSE; } return kTRUE; } Bool_t MSimRandomPhotons::ReInit(MParList *pList) { // Overwrite the default set by MGeomApply fRates->Init(*fGeom); return kTRUE; } // -------------------------------------------------------------------------- // // Check for the necessary containers // Int_t MSimRandomPhotons::Process() { // Get array from event container // const Int_t num = fEvt->GetNumPhotons(); // // Do not produce pure pedestal events! // if (num==0) // return kTRUE; // Get array from event container // FIXME: Use statistics container instead const UInt_t npix = fGeom->GetNumPixels(); // This is the possible window in which the triggered digitization // may take place. const Double_t start = fStat->GetTimeFirst(); const Double_t end = fStat->GetTimeLast(); // Loop over all pixels for (UInt_t idx=0; idxend) break; // Add a new photon // FIXME: SLOW! MPhotonData &ph = fEvt->Add(); // Set source to NightSky, time to t and tag to pixel index ph.SetPrimary(MMcEvtBasic::kNightSky); ph.SetWeight(); ph.SetTime(t); ph.SetTag(idx); // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant) FIXME: Reset? if (fSimulateWavelength) { const Float_t wmin = fRunHeader->GetWavelengthMin(); const Float_t wmax = fRunHeader->GetWavelengthMax(); ph.SetWavelength(TMath::Nint(gRandom->Uniform(wmin, wmax))); } } } // Re-sort the photons by time! fEvt->Sort(kTRUE); // Update maximum index fStat->SetMaxIndex(npix-1); // Shrink return kTRUE; } // -------------------------------------------------------------------------- // // Read the parameters from the resource file. // // FrequencyFixed: 0.040 // FrequencyNSB: 5.8 // // The fixed frequency is given in units fitting the units of the time. // Usually the time is given in nanoseconds thus, e.g., 0.040 means 40MHz. // // The FrequencyNSB is scaled by the area of the pixel in cm^2. Therefore // 0.040 would mean 40MHz/cm^2 // Int_t MSimRandomPhotons::ReadEnv(const TEnv &env, TString prefix, Bool_t print) { Bool_t rc = kFALSE; if (IsEnvDefined(env, prefix, "FrequencyFixed", print)) { rc = kTRUE; fFreqFixed = GetEnvValue(env, prefix, "FrequencyFixed", fFreqFixed); } if (IsEnvDefined(env, prefix, "FrequencyNSB", print)) { rc = kTRUE; fFreqNSB = GetEnvValue(env, prefix, "FrequencyNSB", fFreqNSB); } 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; }