Ignore:
Timestamp:
12/10/09 10:08:55 (15 years ago)
Author:
tbretz
Message:
*** empty log message ***
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/MagicSoft/Mars/msimcamera/MSimRandomPhotons.cc

    r9521 r9525  
    5252//   R = R0 * Ai * f
    5353//
    54 // R0 is the night sky background rate as given in Eckart's paper. Ai the
    55 // area of the cones acceptance window, f is given as:
     54// R0 is the night sky background rate as given in Eckart's paper (divided
     55// by the wavelength window). Ai the area of the cones acceptance window,
     56// f is given as:
    5657//
    5758//   f = nm * Min(Ar, sr*d^2)
    … …  
    6869//   sr is the effective solid angle corresponding to the integral of the
    6970//   cones angular acceptance
     71//
     72// Alternatively, the night-sky background rate can be calculated from
     73// a spectrum as given in Fig. 1 (but versus Nanometers) in
     74//
     75//   Chris R. Benn & Sara L. Ellison La Palma Night-Sky Brightness
     76//
     77// After proper conversion of the units, the rate of the pixel 0
     78// is then calculated by
     79//
     80//     rate = f * nsb
     81//
     82// With nsb
     83//
     84//   nsb = Integral(nsb spectrum * combines efficiencies)
     85//
     86// and f can be either
     87//
     88//   Eff. angular acceptance Cones (e.g. 20deg) * Cone-Area (mm^2)
     89//   f = sr * A0
     90//
     91// or
     92//
     93//   Mirror-Area * Field of view of cones (deg^2)
     94//   f = Ar * A0;
    7095//
    7196//
    … …  
    118143MSimRandomPhotons::MSimRandomPhotons(const char* name, const char *title)
    119144    : fGeom(0), fEvt(0), fStat(0), /*fEvtHeader(0),*/ fRunHeader(0),
    120     fRates(0), fSimulateWavelength(kFALSE),  fNameGeomCam("MGeomCam")
     145    fRates(0), fSimulateWavelength(kFALSE), fNameGeomCam("MGeomCam"),
     146    fFileNameNSB("resmc/night-sky-la-palma.txt")
    121147{
    122148    fName  = name  ? name  : "MSimRandomPhotons";
    … …  
    169195    }*/
    170196
    171     fRunHeader = 0;
    172     if (fSimulateWavelength)
    173     {
    174         fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
    175         if (!fRunHeader)
     197    fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
     198    if (fSimulateWavelength && !fRunHeader)
     199    {
     200        *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
     201        return kFALSE;
     202    }
     203
     204    MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
     205    if (!r)
     206    {
     207        *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
     208        return kFALSE;
     209    }
     210
     211    const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
     212    const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesTransmission",         "MParSpline");
     213    const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity",        "MParSpline");
     214    const MParSpline *s4 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
     215
     216    // Multiply all relevant efficiencies to get the total tarnsmission
     217    MParSpline *eff = (MParSpline*)s1->Clone();
     218    eff->Multiply(*s2->GetSpline());
     219    eff->Multiply(*s3->GetSpline());
     220
     221    // Effectively transmitted wavelength band
     222    const Double_t nm = eff && eff->GetSpline() ? eff->GetSpline()->Integral() : 1;
     223
     224    // Angular acceptance of the cones
     225    const Double_t sr = s4 && s4->GetSpline() ? s4->GetSpline()->IntegralSolidAngle() : 1;
     226
     227    {
     228        const Double_t d2   = fGeom->GetCameraDist()*fGeom->GetCameraDist();
     229        const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();
     230        const Double_t f1   = TMath::Min(r->GetA()/1e4, sr*d2) * conv*conv;
     231
     232        // Rate in GHz / mm^2
     233        fScale = fFreqNSB * nm * f1; // [GHz/mm^2] efficiency * m^2 *rad^2 *mm^2
     234
     235        const Double_t freq0 = fScale*(*fGeom)[0].GetA()*1000;
     236
     237        *fLog << inf << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", freq0) << "MHz" << endl;
     238
     239        // FIXME: Scale with the number of pixels
     240        if (freq0>1000)
    176241        {
    177             *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
     242            *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;
    178243            return kFALSE;
    179244        }
    180245    }
    181246
    182     MReflector *r = (MReflector*)pList->FindObject("Reflector", "MReflector");
    183     if (!r)
    184     {
    185         *fLog << err << "Reflector [MReflector] not found... aborting." << endl;
    186         return kFALSE;
    187     }
    188 
    189     *fLog << warn << "FIXME: A check is needed that the PDE is 0 at both ends!" << endl;
    190 
    191     const MParSpline *s1 = (MParSpline*)pList->FindObject("PhotonDetectionEfficiency", "MParSpline");
    192     const MParSpline *s2 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
    193     const MParSpline *s3 = (MParSpline*)pList->FindObject("MirrorReflectivity",        "MParSpline");
    194     const MParSpline *s5 = (MParSpline*)pList->FindObject("ConesTransmission",         "MParSpline");
    195 
    196     const Double_t d2  = fGeom->GetCameraDist()*fGeom->GetCameraDist();
    197 //    const Double_t pde = s1 && s1->GetSpline() ? s1->GetSpline()->Integral() : 1;
    198     const Double_t sr  = s2 && s2->GetSpline() ? s2->GetSpline()->IntegralSolidAngle() : 1;
    199 //    const Double_t mir = s3 && s3->GetSpline() ? s3->GetSpline()->Integral() : 1;
    200     const Double_t Ar  = r->GetA()/1e4;
    201 
    202     // Conversion factor to convert pixel area to steradians (because it
    203     // is a rather small area we can assume it is flat)
    204     const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad();
    205 
    206     // Multiply all relevant efficiencies
    207     MParSpline *s4 = (MParSpline*)s1->Clone();
    208     s4->Multiply(*s3->GetSpline());
    209     s4->Multiply(*s5->GetSpline());
    210 
    211     const Double_t nm = s4 && s4->GetSpline() ? s4->GetSpline()->Integral() : 1;
    212 
    213     delete s4;
    214 
    215     // /100 to convert the pixel area from mm^2 to cm^2
    216     fScale =  nm * TMath::Min(Ar, sr*d2) * conv*conv;
    217 
    218     *fLog << inf;
    219     *fLog << "Effective cone acceptance:      " << Form("%.2f", sr*d2) << "m^2" << endl;
    220     *fLog << "Reflector area:                 " << Form("%.2f", Ar) << "m^2" << endl;
    221     *fLog << "Resulting eff. collection area: " << Form("%.1f", TMath::Min(Ar, sr*d2)) << "m^2" << endl;
    222 //    *fLog << "Eff. wavelength band (PDE):     " << Form("%.1f", pde) << "nm" << endl;
    223 //    *fLog << "Eff. wavelength band (Mirror):  " << Form("%.1f", mir) << "nm" << endl;
    224     *fLog << "Eff. wavelength band (MIR+Cone+PDE): " << Form("%.1f", nm) << "nm" << endl;
    225     *fLog << "Pixel area of " << fNameGeomCam << "[0]: " << Form("%.2e", (*fGeom)[0].GetA()*conv*conv) << "sr" << endl;
    226     //*fLog << "Effective angular acceptance:     " << sr << "sr" << endl;
    227     //*fLog << "Resulting NSB frequency:           " << fFreqNSB*nm*Ar*1000 << "MHz/sr" << endl;
    228     *fLog << "Resulting Freq. in " << fNameGeomCam << "[0]: " << Form("%.2f", fFreqNSB*(*fGeom)[0].GetA()*fScale*1000) << "MHz" << endl;
     247    if (fFileNameNSB.IsNull())
     248    {
     249        delete eff;
     250        return kTRUE;
     251    }
    229252
    230253    // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
    … …  
    237260    // ampl = MMcFadcHeader->GetAmplitud()
    238261    // sqrt(pedrms*pedrms + dnsbpix*ampl*ampl/ratio)
     262
     263    // Conversion of the y-axis
     264    // ------------------------
     265    // Double_t ff = 1;                               // myJy / arcsec^2 per nm
     266    // ff *= 1e-6;                                    // Jy   / arcsec^2 per nm
     267    // ff *= 3600*3600;                               // Jy   / deg^2
     268    // ff *= 1./TMath::DegToRad()/TMath::DegToRad();  // Jy/sr = 1e-26J/s/m^2/Hz/sr
     269    // ff *= 1e-26;                                   // J/s/m^2/Hz/sr   per nm
     270
     271    const Double_t arcsec2rad = TMath::DegToRad()/3600.;
     272    const Double_t f = 1e-32 / (arcsec2rad*arcsec2rad);
     273
     274    // Read night sky background flux from file
     275    MParSpline flux;
     276    if (!flux.ReadFile(fFileNameNSB))
     277        return kFALSE;
     278
     279    const Int_t min = TMath::FloorNint(flux.GetXmin());
     280    const Int_t max = TMath::CeilNint( flux.GetXmax());
     281
     282    if (fRunHeader)
     283    {
     284        if (min>fRunHeader->GetWavelengthMin())
     285        {
     286            *fLog << warn << "WARNING - Minimum wavelength of night sky background flux (";
     287            *fLog << min << "nm) from " << fFileNameNSB;
     288            *fLog << " exceeds minimum wavelength simulated ";
     289            *fLog << fRunHeader->GetWavelengthMin() << "nm." << endl;
     290        }
     291        if (max<fRunHeader->GetWavelengthMax())
     292        {
     293            *fLog << warn << "WARNING - Maximum wavelength of night sky background flux (";
     294            *fLog << max << "nm) from " << fFileNameNSB;
     295            *fLog << " undershoots maximum wavelength simulated ";
     296            *fLog << fRunHeader->GetWavelengthMax() << "nm." << endl;
     297        }
     298    }
     299
     300    MParSpline nsb;
     301
     302    // Normalization to our units,
     303    // conversion from energy flux to photon flux
     304    nsb.SetFunction(Form("%.12e/(x*TMath::H())", f), max-min, min, max);
     305
     306    // multiply night sky background flux with normalization
     307    nsb.Multiply(*flux.GetSpline());
     308
     309    // Multiply with the total transmission
     310    nsb.Multiply(*eff->GetSpline());
     311
     312    // Check if the photon flux is zero at both ends
     313    if (nsb.GetSpline()->Eval(min)>1e-5)
     314    {
     315        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << min << "nm is not zero... abort." << endl;
     316        return kFALSE;
     317    }
     318    if (nsb.GetSpline()->Eval(max)>1e-5)
     319    {
     320        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << max << "nm is not zero... abort." << endl;
     321        return kFALSE;
     322    }
     323
     324    if (fRunHeader)
     325    {
     326        if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMin())>1e-5)
     327            *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at " << fRunHeader->GetWavelengthMin() << "nm is not zero... abort." << endl;
     328        if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMax())>1e-5)
     329            *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at " << fRunHeader->GetWavelengthMax() << "nm is not zero... abort." << endl;
     330    }
     331
     332    // Conversion from m to radians
     333    const Double_t conv = fGeom->GetConvMm2Deg()*TMath::DegToRad()*1e3;
     334
     335    // Angular acceptance of the cones
     336    //const Double_t sr = s5.GetSpline()->IntegralSolidAngle(); // sr
     337    // Absolute reflector area
     338    const Double_t Ar = r->GetA()/1e4;                      // m^2
     339    // Size of the cone's entrance window
     340    const Double_t A0 = (*fGeom)[0].GetA()*1e-6;                  // m^2
     341
     342    // Rate * m^2 * Solid Angle
     343    // -------------------------
     344
     345    // Angular acceptance Cones (e.g. 20deg) * Cone-Area
     346    const Double_t f1 = A0 * sr;                // m^2 sr
     347
     348    // Mirror-Area * Field of view of cones (e.g. 0.1deg)
     349    const Double_t f2 = Ar * A0*conv*conv;      // m^2 sr
     350
     351    // FIXME: Calculate the reflectivity of the bottom by replacing
     352    // MirrorReflectivity by bottom reflectivity and reflect
     353    // and use it to reflect the difference between f1 and f2
     354    // if any.
     355
     356    // Total NSB rate in MHz per m^2 and sr
     357    const Double_t rate = nsb.GetSpline()->Integral() * 1e-6;
     358
     359    *fLog << inf;
     360
     361    // Resulting rates as if Razmick's constant had been used
     362    // *fLog << 1.75e6/(600-300) * f1 * eff->GetSpline()->Integral() << " MHz" << endl;
     363    // *fLog << 1.75e6/(600-300) * f2 * eff->GetSpline()->Integral() << " MHz" << endl;
     364
     365    *fLog << "Conversion factor Fnu:      " << f  << endl;
     366    *fLog << "Total reflective area:      " << Form("%.2f", Ar) << " m" << UTF8::kSquare << endl;
     367    *fLog << "Acceptance area of cone 0:  " << Form("%.2f", A0*1e6) << " mm" << UTF8::kSquare << " = ";
     368    *fLog << A0*conv*conv << " sr" << endl;
     369    *fLog << "Cones angular acceptance:   " << sr << " sr" << endl;
     370    *fLog << "ConeArea*MirrorAngle (f1):  " << f1 << " m^2 sr" << endl;
     371    *fLog << "MirrorArea*ConeAngle (f2):  " << f2 << " m^2 sr" << endl;
     372    *fLog << "Effective. transmission:    " << Form("%.1f", nm) << " nm" << endl;
     373    *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f1): " << Form("%.2f", rate * f1) << " MHz" << endl;
     374    *fLog << "NSB freq. in " << fNameGeomCam << "[0] (f2): " << Form("%.2f", rate * f2) << " MHz" << endl;
     375    *fLog << "Using f1." << endl;
     376
     377    // Scale the rate per mm^2 and to GHz
     378    fScale = rate * f1 / (*fGeom)[0].GetA() / 1000;
     379
     380    // FIXME: Scale with the number of pixels
     381    if (rate*f1>1000)
     382    {
     383        *fLog << err << "ERROR - Frequency exceeds 1GHz, this might leed to too much memory consumption." << endl;
     384        return kFALSE;
     385    }
     386
     387    delete eff;
    239388
    240389    return kTRUE;
    … …  
    274423    {
    275424        // Scale the rate with the pixel size.
    276         const Double_t rate = fFreqFixed+fFreqNSB*(*fGeom)[idx].GetA()*fScale;
     425        const Double_t rate = fFreqFixed + fScale*(*fGeom)[idx].GetA();
    277426
    278427        (*fRates)[idx].SetPedestal(rate);
    … …  
    305454            // fProductionHeight, fPosX, fPosY, fCosU, fCosV (irrelevant)  FIXME: Reset?
    306455
    307             if (fRunHeader)
     456            if (fSimulateWavelength)
    308457            {
    309458                const Float_t wmin = fRunHeader->GetWavelengthMin();
    … …  
    353502    }
    354503
     504    if (IsEnvDefined(env, prefix, "FileNameNSB", print))
     505    {
     506        rc = kTRUE;
     507        fFileNameNSB = GetEnvValue(env, prefix, "FileNameNSB", fFileNameNSB);
     508    }
     509
     510    if (IsEnvDefined(env, prefix, "SimulateCherenkovSpectrum", print))
     511    {
     512        rc = kTRUE;
     513        fSimulateWavelength = GetEnvValue(env, prefix, "SimulateCherenkovSpectrum", fSimulateWavelength);
     514    }
     515
    355516    return rc;
    356517}
Note: See TracChangeset for help on using the changeset viewer.