ItemizeSimulation.cpp

Go to the documentation of this file.

//  ************************************************************************************************
//
//  BornAgain: simulate and fit reflection and scattering
//
//! @file      GUI/Model/FromCore/ItemizeSimulation.cpp
//! @brief     Implements functions that convert ISimulation from core to GUI items
//!
//! @homepage  http://www.bornagainproject.org
//! @license   GNU General Public License v3 or higher (see COPYING)
//! @copyright Forschungszentrum Jülich GmbH 2018
//! @authors   Scientific Computing Group at MLZ (see CITATION, AUTHORS)
//
//  ************************************************************************************************
 
#include "GUI/Model/FromCore/ItemizeSimulation.h"
#include "Base/Const/Units.h"
#include "Base/Util/Assert.h"
#include "Device/Beam/Beam.h"
#include "Device/Beam/FootprintGauss.h"
#include "Device/Beam/FootprintSquare.h"
#include "Device/Detector/RectangularDetector.h"
#include "Device/Detector/SphericalDetector.h"
#include "Device/Mask/DetectorMask.h"
#include "Device/Mask/Ellipse.h"
#include "Device/Mask/InfinitePlane.h"
#include "Device/Mask/Line.h"
#include "Device/Mask/Polygon.h"
#include "Device/Mask/Rectangle.h"
#include "Device/Resolution/ConvolutionDetectorResolution.h"
#include "Device/Resolution/ResolutionFunction2DGaussian.h"
#include "GUI/Model/Device/AxesItems.h"
#include "GUI/Model/Device/BackgroundItems.h"
#include "GUI/Model/Device/BeamAngleItems.h"
#include "GUI/Model/Device/BeamWavelengthItem.h"
#include "GUI/Model/Device/InstrumentCollection.h"
#include "GUI/Model/Device/InstrumentItems.h"
#include "GUI/Model/Device/MaskItems.h"
#include "GUI/Model/Device/RectangularDetectorItem.h"
#include "GUI/Model/Device/ResolutionFunctionItems.h"
#include "GUI/Model/Device/SphericalDetectorItem.h"
#include "GUI/Model/FromCore/ItemizeSample.h"
#include "GUI/Model/Sample/LayerItem.h"
#include "GUI/Model/Sample/MesoCrystalItem.h"
#include "GUI/Model/Sample/ParticleItem.h"
#include "GUI/Model/Sample/ParticleLayoutItem.h"
#include "GUI/Support/Data/SimulationOptionsItem.h"
#include "GUI/Util/Error.h"
#include "Param/Distrib/Distributions.h"
#include "Param/Distrib/RangedDistributions.h"
#include "Param/Node/NodeUtils.h"
#include "Resample/Options/SimulationOptions.h"
#include "Sim/Background/ConstantBackground.h"
#include "Sim/Background/PoissonBackground.h"
#include "Sim/Scan/AlphaScan.h"
#include "Sim/Scan/ScanResolution.h"
#include "Sim/Simulation/includeSimulations.h"
 
namespace {
 
void setMaskContainer(MaskItems* destMaskItems, const IDetector& detector, double scale)
{
    const auto* detectorMask = detector.detectorMask();
    for (size_t i_mask = 0; i_mask < detectorMask->numberOfMasks(); ++i_mask) {
        const MaskPattern* pat = detectorMask->patternAt(i_mask);
        IShape2D* shape = pat->shape;
        bool mask_value = pat->doMask;
 
        if (const auto* ellipse = dynamic_cast<const Ellipse*>(shape)) {
            auto* ellipseItem = new EllipseItem();
            ellipseItem->setXCenter(scale * ellipse->getCenterX());
            ellipseItem->setYCenter(scale * ellipse->getCenterY());
            ellipseItem->setXRadius(scale * ellipse->radiusX());
            ellipseItem->setYRadius(scale * ellipse->radiusY());
            ellipseItem->setAngle(scale * ellipse->getTheta());
            ellipseItem->setMaskValue(mask_value);
            destMaskItems->insertMask(0, ellipseItem);
        }
 
        else if (const auto* rectangle = dynamic_cast<const Rectangle*>(shape)) {
            auto* rectangleItem = new RectangleItem();
            rectangleItem->setXLow(scale * rectangle->getXlow());
            rectangleItem->setYLow(scale * rectangle->getYlow());
            rectangleItem->setXUp(scale * rectangle->getXup());
            rectangleItem->setYUp(scale * rectangle->getYup());
            rectangleItem->setMaskValue(mask_value);
            destMaskItems->insertMask(0, rectangleItem);
        }
 
        else if (const auto* polygon = dynamic_cast<const Polygon*>(shape)) {
            auto* polygonItem = new PolygonItem();
            std::vector<double> xpos, ypos;
            polygon->getPoints(xpos, ypos);
            for (size_t i_point = 0; i_point < xpos.size(); ++i_point)
                polygonItem->addPoint(scale * xpos[i_point], scale * ypos[i_point]);
 
            polygonItem->setMaskValue(mask_value);
            polygonItem->setIsClosed(true);
 
            destMaskItems->insertMask(0, polygonItem);
        }
 
        else if (const auto* vline = dynamic_cast<const VerticalLine*>(shape)) {
            auto* lineItem = new VerticalLineItem();
            lineItem->setPosX(scale * vline->getXpos());
            lineItem->setMaskValue(mask_value);
            destMaskItems->insertMask(0, lineItem);
        }
 
        else if (const auto* hline = dynamic_cast<const HorizontalLine*>(shape)) {
            auto* lineItem = new HorizontalLineItem();
            lineItem->setPosY(scale * hline->getYpos());
            lineItem->setMaskValue(mask_value);
            destMaskItems->insertMask(0, lineItem);
        }
 
        else if (const auto* plane = dynamic_cast<const InfinitePlane*>(shape)) {
            Q_UNUSED(plane);
            auto* planeItem = new MaskAllItem();
            planeItem->setMaskValue(mask_value);
            destMaskItems->insertMask(-1, planeItem);
        }
 
        else
            throw Error("Cannot convert detector mask from core to GUI: Unknown shape");
    }
 
    if (detector.hasExplicitRegionOfInterest()) {
        const auto xBounds = detector.regionOfInterestBounds(0);
        const auto yBounds = detector.regionOfInterestBounds(1);
 
        auto* roiItem = new RegionOfInterestItem();
        roiItem->setXLow(scale * xBounds.first);
        roiItem->setYLow(scale * yBounds.first);
        roiItem->setXUp(scale * xBounds.second);
        roiItem->setYUp(scale * yBounds.second);
        destMaskItems->insertMask(-1, roiItem);
    }
}
 
void setDetectorMasks(DetectorItem* detector_item, const ISimulation2D& simulation)
{
    const IDetector& detector = simulation.detector();
    if ((detector.detectorMask() && detector.detectorMask()->hasMasks())
        || detector.hasExplicitRegionOfInterest()) {
        const double scale = 1.0 / detector_item->axesToCoreUnitsFactor();
        setMaskContainer(&detector_item->maskItems(), detector, scale);
    }
}
 
void setFootprintFactor(const IFootprintFactor* footprint, SpecularBeamItem* beam_item)
{
    if (!footprint)
        return;
    if (const auto* const gaussian_fp = dynamic_cast<const FootprintGauss*>(footprint))
        beam_item->setGaussianFootprint(gaussian_fp->widthRatio());
    else if (const auto* const square_fp = dynamic_cast<const FootprintSquare*>(footprint))
        beam_item->setSquareFootprint(square_fp->widthRatio());
}
 
void setDistribution(BeamDistributionItem* part_distr_item, ParameterDistribution par_distr,
                     double factor)
{
    const IDistribution1D* p_distribution = par_distr.getDistribution();
 
    DistributionItem* distItem = nullptr;
    if (const auto* distr = dynamic_cast<const DistributionGate*>(p_distribution)) {
        auto* distr_gate_item = part_distr_item->setDistributionType<DistributionGateItem>();
        distr_gate_item->setRange(factor * distr->min(), factor * distr->max());
        distItem = distr_gate_item;
    } else if (const auto* distr = dynamic_cast<const DistributionLorentz*>(p_distribution)) {
        auto* distr_lorentz_item = part_distr_item->setDistributionType<DistributionLorentzItem>();
        distr_lorentz_item->setMean(factor * distr->mean());
        distr_lorentz_item->setHwhm(factor * distr->hwhm());
        distItem = distr_lorentz_item;
    } else if (const auto* distr = dynamic_cast<const DistributionGaussian*>(p_distribution)) {
        auto* distr_gauss_item = part_distr_item->setDistributionType<DistributionGaussianItem>();
        distr_gauss_item->setMean(factor * distr->mean());
        distr_gauss_item->setStandardDeviation(factor * distr->getStdDev());
        distItem = distr_gauss_item;
    } else if (const auto* distr = dynamic_cast<const DistributionLogNormal*>(p_distribution)) {
        auto* distr_lognorm_item =
            part_distr_item->setDistributionType<DistributionLogNormalItem>();
        distr_lognorm_item->setMedian(factor * distr->getMedian());
        distr_lognorm_item->setScaleParameter(distr->getScalePar());
        distItem = distr_lognorm_item;
    } else if (const auto* distr = dynamic_cast<const DistributionCosine*>(p_distribution)) {
        auto* distr_cos_item = part_distr_item->setDistributionType<DistributionCosineItem>();
        distr_cos_item->setMean(factor * distr->mean());
        distr_cos_item->setSigma(factor * distr->sigma());
        distItem = distr_cos_item;
    } else if (const auto* distr = dynamic_cast<const DistributionTrapezoid*>(p_distribution)) {
        auto* distr_trapez_item = part_distr_item->setDistributionType<DistributionTrapezoidItem>();
        distr_trapez_item->setCenter(factor * distr->mean());
        distr_trapez_item->setLeftWidth(factor * distr->getLeftWidth());
        distr_trapez_item->setMiddleWidth(factor * distr->getMiddleWidth());
        distr_trapez_item->setRightWidth(factor * distr->getRightWidth());
        distItem = distr_trapez_item;
    } else
        throw Error("setDistribution: -> unknown distribution");
 
    distItem->setNumberOfSamples((int)par_distr.nDraws());
 
    if (distItem->hasSigmaFactor())
        distItem->setSigmaFactor(par_distr.sigmaFactor());
 
    // TODO It's wrong if domain distribution made for angles.
    distItem->setLimits(par_distr.getLimits().scaledLimits(factor));
}
 
void setItemFromSample(BeamDistributionItem* beam_distribution_item,
                       const ParameterDistribution& parameter_distribution)
{
    ASSERT(beam_distribution_item);
 
    if (parameter_distribution.getMinValue() < parameter_distribution.getMaxValue()) {
        throw Error(
            "setItemFromSample(BeamDistributionItem* "
            "beamDistributionItem,"
            "const ParameterDistribution& parameterDistribution) -> Error. ParameterDistribution "
            "with defined min,max are not yet implemented in GUI");
    }
 
    const double unit_factor = 1.0 / beam_distribution_item->scaleFactor();
    setDistribution(beam_distribution_item, parameter_distribution, unit_factor);
}
 
void addDistributionToBeamItem(ParameterDistribution::WhichParameter which,
                               BeamDistributionItem* item,
                               const ParameterDistribution& distribution)
{
    if (distribution.whichParameter() != which)
        return;
 
    setItemFromSample(item, distribution);
}
 
void addRangedDistributionToItem(BeamDistributionItem* distributionItem,
                                 const IRangedDistribution& ranged, double mean, double std_dev)
{
    if (!distributionItem)
        return;
    ParameterDistribution par_distr(ParameterDistribution::None,
                                    *ranged.distribution(mean, std_dev), ranged.nSamples(),
                                    ranged.sigmaFactor(), ranged.limits());
    setItemFromSample(distributionItem, par_distr);
}
 
void setGISASBeamItem(BeamItem* beam_item, const ScatteringSimulation& simulation)
{
    ASSERT(beam_item);
    const Beam& beam = simulation.beam();
 
    beam_item->setIntensity(beam.intensity());
    beam_item->setWavelength(beam.wavelength());
    beam_item->setInclinationAngle(Units::rad2deg(beam.direction().alpha()));
    beam_item->setAzimuthalAngle(Units::rad2deg(beam.direction().phi()));
 
    // distribution parameters
    const std::vector<ParameterDistribution>& distributions = simulation.getDistributions();
    for (size_t i = 0; i < distributions.size(); ++i) {
        // try all, one will have the right 'whichParameter'
        addDistributionToBeamItem(ParameterDistribution::BeamWavelength,
                                  beam_item->wavelengthItem(), distributions[i]);
        addDistributionToBeamItem(ParameterDistribution::BeamInclinationAngle,
                                  beam_item->inclinationAngleItem(), distributions[i]);
        addDistributionToBeamItem(ParameterDistribution::BeamAzimuthalAngle,
                                  beam_item->azimuthalAngleItem(), distributions[i]);
    }
}
 
void setOffspecBeamItem(BeamItem* beam_item, const OffspecSimulation& simulation)
{
    const Beam& beam = simulation.beam();
 
    beam_item->setIntensity(beam.intensity());
    beam_item->setWavelength(beam.wavelength());
    beam_item->setInclinationAngle(Units::rad2deg(beam.direction().alpha()));
    beam_item->setAzimuthalAngle(Units::rad2deg(beam.direction().phi()));
    // TODO implement beam divergence
}
 
void setAxisItem(BasicAxisItem* item, const IAxis& axis, double factor)
{
    if (!dynamic_cast<const FixedBinAxis*>(&axis))
        throw Error("setAxisItem() -> Error. Unexpected axis");
 
    item->setBinCount(static_cast<int>(axis.size()));
    item->setLowerBound(factor * axis.min());
    item->setUpperBound(factor * axis.max());
    item->setTitle(QString::fromStdString(axis.axisName()));
}
 
void setSphericalDetector(SphericalDetectorItem* detectorItem, const SphericalDetector& detector)
{
    // Axes
    const IAxis& phi_axis = detector.axis(0);
    const IAxis& alpha_axis = detector.axis(1);
 
    auto& phiAxis = detectorItem->phiAxis();
    phiAxis.setNbins(phi_axis.size());
    phiAxis.setMin(Units::rad2deg(phi_axis.min()));
    phiAxis.setMax(Units::rad2deg(phi_axis.max()));
 
    auto& alphaAxis = detectorItem->alphaAxis();
    alphaAxis.setNbins(alpha_axis.size());
    alphaAxis.setMin(Units::rad2deg(alpha_axis.min()));
    alphaAxis.setMax(Units::rad2deg(alpha_axis.max()));
}
 
void setRectangularDetector(RectangularDetectorItem* detectorItem,
                            const RectangularDetector& detector)
{
    // Axes
    detectorItem->setXSize(detector.xSize());
    detectorItem->setWidth(detector.width());
 
    detectorItem->setYSize(detector.ySize());
    detectorItem->setHeight(detector.height());
 
    detectorItem->setDetectorAlignment(detector.getDetectorArrangment());
    if (detector.getDetectorArrangment() == RectangularDetector::GENERIC) {
        R3 normal = detector.getNormalVector();
        detectorItem->setNormalVector(normal);
 
        R3 direction = detector.getDirectionVector();
        detectorItem->setDirectionVector(direction);
 
        detectorItem->setU0(detector.getU0());
        detectorItem->setV0(detector.getV0());
    }
 
    else if (detector.getDetectorArrangment() == RectangularDetector::PERPENDICULAR_TO_SAMPLE) {
        detectorItem->setDistance(detector.getDistance());
        detectorItem->setU0(detector.getU0());
        detectorItem->setV0(detector.getV0());
    } else if (detector.getDetectorArrangment()
               == RectangularDetector::PERPENDICULAR_TO_DIRECT_BEAM) {
        detectorItem->setDistance(detector.getDistance());
        detectorItem->setU0(detector.getU0());
        detectorItem->setV0(detector.getV0());
    } else if (detector.getDetectorArrangment()
               == RectangularDetector::PERPENDICULAR_TO_REFLECTED_BEAM) {
        detectorItem->setDistance(detector.getDistance());
        detectorItem->setU0(detector.getU0());
        detectorItem->setV0(detector.getV0());
    } else if (detector.getDetectorArrangment()
               == RectangularDetector::PERPENDICULAR_TO_REFLECTED_BEAM_DPOS) {
        detectorItem->setDistance(detector.getDistance());
        detectorItem->setU0(detector.getDirectBeamU0());
        detectorItem->setV0(detector.getDirectBeamV0());
    } else {
        throw Error("setItemFromSample(RectangularDetectorItem* "
                    "detectorItem "
                    "Error. Unknown detector arrangement");
    }
}
 
void setDetectorGeometry(Instrument2DItem* instrument_item, const IDetector& detector)
{
    if (const auto* det = dynamic_cast<const SphericalDetector*>(&detector)) {
        auto* item = instrument_item->setDetectorType<SphericalDetectorItem>();
        setSphericalDetector(item, *det);
    } else if (const auto* det = dynamic_cast<const RectangularDetector*>(&detector)) {
        auto* item = instrument_item->setDetectorType<RectangularDetectorItem>();
        setRectangularDetector(item, *det);
    } else {
        throw Error("setDetectorGeometry() -> Unknown detector type.");
    }
}
 
void setDetectorResolution(DetectorItem* detector_item, const IDetector& detector)
{
    const IDetectorResolution* p_resfunc = detector.detectorResolution();
 
    if (!p_resfunc)
        return;
 
    if (const auto* p_convfunc = dynamic_cast<const ConvolutionDetectorResolution*>(p_resfunc)) {
        if (const auto* resfunc = dynamic_cast<const ResolutionFunction2DGaussian*>(
                p_convfunc->getResolutionFunction2D())) {
            auto* item =
                detector_item->setResolutionFunctionType<ResolutionFunction2DGaussianItem>();
            const double scale = 1.0 / detector_item->axesToCoreUnitsFactor();
            item->setSigmaX(scale * resfunc->sigmaX());
            item->setSigmaY(scale * resfunc->sigmaY());
        } else {
            throw Error("setDetectorResolution() -> Error. "
                        "Unknown detector resolution function");
        }
    } else {
        throw Error("setDetectorResolution() -> Error. "
                    "Not a ConvolutionDetectorResolution function");
    }
}
 
void setPolarizerAnalyzer(Instrument2DItem* instrument_item, const ISimulation2D& simulation)
{
    instrument_item->setPolarization(simulation.beam().polVector());
 
    const IDetector& detector = simulation.detector();
    double total_transmission = detector.analyzer().totalTransmission();
    if (total_transmission <= 0.0)
        return;
 
    R3 analyzer_dir = detector.analyzer().analyzerDirection();
    double efficiency = detector.analyzer().analyzerEfficiency();
    instrument_item->setAnalyzerDirection(analyzer_dir);
    instrument_item->setAnalyzerEfficiency(efficiency);
    instrument_item->setAnalyzerTotalTransmission(total_transmission);
}
 
void updateDetector(Instrument2DItem* instrument_item, const ISimulation2D& simulation)
{
    const IDetector& detector = simulation.detector();
    setDetectorGeometry(instrument_item, detector);
 
    auto* detector_item = instrument_item->detectorItem();
 
    setDetectorResolution(detector_item, detector);
    setDetectorMasks(detector_item, simulation);
    setPolarizerAnalyzer(instrument_item, simulation);
}
 
void setBackground(InstrumentItem* instrument_item, const ISimulation& simulation)
{
    const auto* p_bg = simulation.background();
    if (const auto* p_constant_bg = dynamic_cast<const ConstantBackground*>(p_bg)) {
        auto* constant_bg_item = instrument_item->setBackgroundType<ConstantBackgroundItem>();
        double value = p_constant_bg->backgroundValue();
        constant_bg_item->setBackgroundValue(value);
    } else if (dynamic_cast<const PoissonBackground*>(p_bg))
        instrument_item->setBackgroundType<PoissonBackgroundItem>();
}
 
GISASInstrumentItem* createGISASInstrumentItem(const ScatteringSimulation& simulation)
{
    auto* result = new GISASInstrumentItem();
    setGISASBeamItem(result->beamItem(), simulation);
    updateDetector(result, simulation);
 
    return result;
}
 
OffspecInstrumentItem* createOffspecInstrumentItem(const OffspecSimulation& simulation)
{
    auto* result = new OffspecInstrumentItem();
    setOffspecBeamItem(result->beamItem(), simulation);
    updateDetector(result, simulation);
 
    const double factor = 1. / Units::deg;
    result->alphaAxis().setNbins(simulation.beamAxis()->size());
    result->alphaAxis().setMin(factor * simulation.beamAxis()->min());
    result->alphaAxis().setMax(factor * simulation.beamAxis()->max());
 
    return result;
}
 
SpecularInstrumentItem* createSpecularInstrumentItem(const SpecularSimulation& simulation)
{
    auto* result = new SpecularInstrumentItem();
    SpecularBeamItem* beam_item = result->beamItem();
 
    const ISpecularScan* scan = simulation.scan();
 
    beam_item->setIntensity(1.0); // TODO ISSUE #303
    beam_item->setWavelength(scan->wavelength());
    beam_item->setInclinationAngle(0.0); // inclination angle is hardcoded
    beam_item->setAzimuthalAngle(0.0);   // azimuthal angle is hardcoded
 
    auto* axis_item = beam_item->inclinationAxis();
    setAxisItem(axis_item, *scan->coordinateAxis(), 1. / Units::deg);
 
    setFootprintFactor(scan->footprintFactor(), beam_item);
 
    const auto* alphaScan = dynamic_cast<const AlphaScan*>(scan);
    if (!alphaScan)
        return result;
 
    const auto* resolution = alphaScan->wavelengthResolution();
    if (!resolution->empty()) {
        double mean = alphaScan->wavelength();
        double std_dev = resolution->stdDevs(mean, 1).front();
        addRangedDistributionToItem(beam_item->wavelengthItem(), *resolution->distribution(), mean,
                                    std_dev);
    }
 
    resolution = alphaScan->angleResolution();
    if (resolution && !resolution->empty()) {
        double std_dev = resolution->stdDevs(0.0, 1).front();
        addRangedDistributionToItem(beam_item->inclinationAngleItem(), *resolution->distribution(),
                                    0.0, std_dev);
    }
 
    return result;
}
 
} // namespace
 
 
InstrumentCollection* GUI::Transform::FromCore::itemizeInstruments(const ISimulation& simulation)
{
    InstrumentItem* item;
 
    if (const auto* gisasSimulation = dynamic_cast<const ScatteringSimulation*>(&simulation))
        item = createGISASInstrumentItem(*gisasSimulation);
    else if (const auto* offspecSimulation = dynamic_cast<const OffspecSimulation*>(&simulation))
        item = createOffspecInstrumentItem(*offspecSimulation);
    else if (const auto* spec_simulation = dynamic_cast<const SpecularSimulation*>(&simulation))
        item = createSpecularInstrumentItem(*spec_simulation);
    else
        ASSERT(0);
 
    setBackground(item, simulation);
 
    auto* collection = new InstrumentCollection;
    collection->emplace_back(item);
    return collection;
}
 
SimulationOptionsItem* GUI::Transform::FromCore::itemizeOptions(const ISimulation& simulation)
{
    auto* result = new SimulationOptionsItem;
 
    if (simulation.options().isIntegrate())
        result->setUseMonteCarloIntegration(
            static_cast<unsigned>(simulation.options().getMcPoints()));
    else
        result->setUseAnalytical();
 
    result->setUseAverageMaterials(simulation.options().useAvgMaterials());
    result->setIncludeSpecularPeak(simulation.options().includeSpecular());
 
    return result;
}