Changeset 9913 for trunk/Mars


Ignore:
Timestamp:
08/31/10 10:26:32 (14 years ago)
Author:
tbretz
Message:
Improved the range checks in MSimRandomPhotons.
Location:
trunk/Mars
Files:
3 edited

Legend:

Unmodified
Added
Removed
  • trunk/Mars/Changelog

    r9903 r9913  
    1818
    1919                                                 -*-*- END OF LINE -*-*-
     20
     21 2010/08/31 Thomas Bretz
     22
     23   * msimcamera/MSimRandomPhotons.[h,cc]:
     24     - improved the range checks. All ranges are now aligned on the
     25       range of simulated wavlengths
     26
     27
    2028
    2129 2010/08/30 Thomas Bretz
  • trunk/Mars/msimcamera/MSimRandomPhotons.cc

    r9574 r9913  
    150150}
    151151
     152Bool_t MSimRandomPhotons::CheckWavelengthRange(const MParSpline &sp, const char *txt) const
     153{
     154    const Float_t min = sp.GetXmin();
     155    const Float_t max = sp.GetXmax();
     156
     157    if (min>fRunHeader->GetWavelengthMin())
     158    {
     159        *fLog << err << "ERROR - Minimum wavelength (" << min << "nm)";
     160        *fLog << " defined for " << txt;
     161        *fLog << " exceeds minimum wavelength simulated (";
     162        *fLog << fRunHeader->GetWavelengthMin() << "nm)." << endl;
     163        return kFALSE;
     164    }
     165    if (max<fRunHeader->GetWavelengthMax())
     166    {
     167        *fLog << err << "ERROR - Maximum wavelength (" << max << "nm)";
     168        *fLog << " defined for " << txt;
     169        *fLog << " undershoots maximum wavelength simulated (";
     170        *fLog << fRunHeader->GetWavelengthMax() << "nm)." << endl;
     171        return kFALSE;
     172    }
     173
     174    return kTRUE;
     175}
     176
    152177// --------------------------------------------------------------------------
    153178//
     
    196221
    197222    fRunHeader = (MCorsikaRunHeader*)pList->FindObject("MCorsikaRunHeader");
    198     if (fSimulateWavelength && !fRunHeader)
     223    if (!fRunHeader)
    199224    {
    200225        *fLog << err << "MCorsikaRunHeader not found... aborting." << endl;
     
    214239    const MParSpline *s4 = (MParSpline*)pList->FindObject("ConesAngularAcceptance",    "MParSpline");
    215240
    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;
     241    // Check if all splines are defined in the relevant range
     242    if (!CheckWavelengthRange(*s1, "PhotonDetectionEfficiency [MParSpline]"))
     243        return kFALSE;
     244    if (!CheckWavelengthRange(*s2, "ConesTransmission [MParSpline]"))
     245        return kFALSE;
     246    if (!CheckWavelengthRange(*s3, "MirrorReflectivity [MParSpline]"))
     247        return kFALSE;
     248
     249    const Float_t wmin = fRunHeader->GetWavelengthMin();
     250    const Float_t wmax = fRunHeader->GetWavelengthMax();
     251
     252    const Int_t min = TMath::FloorNint(wmin);
     253    const Int_t max = TMath::CeilNint(wmax);
     254
     255    // Multiply all relevant efficiencies to get the total transmission
     256    MParSpline eff;
     257    eff.SetFunction("1", max-min, min, max);
     258
     259    eff.Multiply(*s1->GetSpline());
     260    eff.Multiply(*s2->GetSpline());
     261    eff.Multiply(*s3->GetSpline());
     262
     263    // Effectively transmitted wavelength band in the simulated range
     264    const Double_t nm = eff.GetSpline()->Integral();
    223265
    224266    // Angular acceptance of the cones
     
    246288
    247289    if (fFileNameNSB.IsNull())
    248     {
    249         delete eff;
    250290        return kTRUE;
    251     }
    252291
    253292    // const MMcRunHeader *mcrunheader = (MMcRunHeader*)pList->FindObject("MMcRunHeader");
     
    277316        return kFALSE;
    278317
    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     }
     318    if (!CheckWavelengthRange(flux, TString("night sky background flux from ")+fFileNameNSB))
     319        return kFALSE;
    299320
    300321    MParSpline nsb;
     
    308329
    309330    // Multiply with the total transmission
    310     nsb.Multiply(*eff->GetSpline());
    311 
    312     // Check if the photon flux is zero at both ends
     331    nsb.Multiply(*eff.GetSpline());
     332
     333    // Check if the photon flux is zero at both ends of the NSB
    313334    if (nsb.GetSpline()->Eval(min)>1e-5)
    314335    {
    315         *fLog << err << "ERROR - Transmitted NSB spectrum at detector at " << min << "nm is not zero... abort." << endl;
    316         return kFALSE;
     336        *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at ";
     337        *fLog << wmin << "nm is not zero, but " << nsb.GetSpline()->Eval(wmin) << "... abort." << endl;
    317338    }
    318339    if (nsb.GetSpline()->Eval(max)>1e-5)
    319340    {
    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;
     341        *fLog << warn << "WARNING - Transmitted NSB spectrum at detector at ";
     342        *fLog << wmax << "nm is not zero, but " << nsb.GetSpline()->Eval(wmax) << "... abort." << endl;
     343    }
     344
     345    // Check if the photon flux is zero at both ends of the simulated region
     346    if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMin())>1e-5)
     347    {
     348        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at minimum simulated wavelength ";
     349        *fLog << fRunHeader->GetWavelengthMin() << "nm is not zero... abort." << endl;
     350        return kFALSE;
     351    }
     352    if (nsb.GetSpline()->Eval(fRunHeader->GetWavelengthMax())>1e-5)
     353    {
     354        *fLog << err << "ERROR - Transmitted NSB spectrum at detector at maximum simulated wavelength ";
     355        *fLog << fRunHeader->GetWavelengthMax() << "nm is not zero... abort." << endl;
     356        return kFALSE;
    330357    }
    331358
     
    336363    //const Double_t sr = s5.GetSpline()->IntegralSolidAngle(); // sr
    337364    // Absolute reflector area
    338     const Double_t Ar = r->GetA()/1e4;                      // m^2
     365    const Double_t Ar = r->GetA()/1e4;                          // m^2
    339366    // Size of the cone's entrance window
    340     const Double_t A0 = (*fGeom)[0].GetA()*1e-6;                  // m^2
     367    const Double_t A0 = (*fGeom)[0].GetA()*1e-6;                // m^2
    341368
    342369    // Rate * m^2 * Solid Angle
     
    360387
    361388    // 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;
     389    // *fLog << 1.75e6/(600-300) * f1 * eff.GetSpline()->Integral() << " MHz" << endl;
     390    // *fLog << 1.75e6/(600-300) * f2 * eff.GetSpline()->Integral() << " MHz" << endl;
    364391
    365392    *fLog << "Conversion factor Fnu:         " << f  << endl;
     
    385412    }
    386413
    387     delete eff;
    388 
    389414    return kTRUE;
    390415}
  • trunk/Mars/msimcamera/MSimRandomPhotons.h

    r9525 r9913  
    88class MGeomCam;
    99class MParList;
     10class MParSpline;
    1011class MPhotonEvent;
    1112class MPhotonStatistics;
    12 //class MCorsikaEvtHeader;
    1313class MCorsikaRunHeader;
    1414class MPedestalCam;
     
    4343    Int_t ReadEnv(const TEnv &env, TString prefix, Bool_t print=kFALSE);
    4444
     45    // MSimRandomPhotons
     46    Bool_t CheckWavelengthRange(const MParSpline &sp, const char *txt) const;
     47
    4548public:
    4649    MSimRandomPhotons(const char *name=NULL, const char *title=NULL);
Note: See TracChangeset for help on using the changeset viewer.