C++ Peak::setGoniometerMatrix方法代码示例

  /** Append the peaks from a .peaks file into the workspace
   * @param outWS :: the workspace in which to place the information
   * @param filename :: path to the .peaks file
   */
  void LoadIsawPeaks::appendFile( PeaksWorkspace_sptr outWS, std::string filename )
  {

    // Open the file
    std::ifstream in( filename.c_str() );


    // Read the header, load the instrument
    double T0;
    std::string s = readHeader( outWS, in , T0);
    // set T0 in the run parameters
    API::Run & m_run = outWS->mutableRun();
    m_run.addProperty<double>("T0", T0, true);

    if( !in.good() || s.length() < 1 )
      throw std::runtime_error( "End of Peaks file before peaks" );

    if( s.compare( std::string( "0" ) ) != 0 )
      throw std::logic_error( "No header for Peak segments"  );

    readToEndOfLine( in ,  true );
    s = getWord( in , false );

    int run, bankNum;
    double chi , phi , omega , monCount;

    // Build the universal goniometer that will build the rotation matrix.
    Mantid::Geometry::Goniometer uniGonio;
    uniGonio.makeUniversalGoniometer();

    // TODO: Can we find the number of peaks to get better progress reporting?
    Progress prog(this, 0.0, 1.0, 100);

    while( in.good() )
    {
      // Read the header if necessary
      s = readPeakBlockHeader( s ,  in  , run , bankNum , chi , phi ,
          omega , monCount );
      // Build the Rotation matrix using phi,chi,omega
      uniGonio.setRotationAngle("phi", phi);
      uniGonio.setRotationAngle("chi", chi);
      uniGonio.setRotationAngle("omega", omega);
      //Put goniometer into peaks workspace
      outWS->mutableRun().setGoniometer(uniGonio, false);


      std::ostringstream oss;
      std::string bankString = "bank";
      if (outWS->getInstrument()->getName() == "WISH") bankString = "WISHpanel0";
      oss << bankString << bankNum;
      std::string bankName = oss.str();

      int seqNum = -1;

      try
      {
        // Read the peak
        Peak peak = readPeak(outWS, s, in, seqNum, bankName);

        // Get the calculated goniometer matrix
        Matrix<double> gonMat = uniGonio.getR();

        peak.setGoniometerMatrix(gonMat);
        peak.setRunNumber(run);
        peak.setMonitorCount( monCount );

        double tof = peak.getTOF();
        Kernel::Units::Wavelength wl;

        wl.initialize(peak.getL1(), peak.getL2(), peak.getScattering(), 0,
                  peak.getInitialEnergy(), 0.0);

        peak.setWavelength(wl.singleFromTOF( tof));
        // Add the peak to workspace
        outWS->addPeak(peak);
      }
      catch (std::runtime_error & e)
      {
        g_log.warning() << "Error reading peak SEQN " << seqNum << " : " << e.what() << std::endl;
      }

      prog.report();
    }

  }

void PeakHKLErrors::functionDeriv1D(Jacobian *out, const double *xValues,
                                    const size_t nData) {
  PeaksWorkspace_sptr Peaks =
      AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(
          PeakWorkspaceName);
  boost::shared_ptr<Geometry::Instrument> instNew = getNewInstrument(Peaks);

  const DblMatrix &UB = Peaks->sample().getOrientedLattice().getUB();
  DblMatrix UBinv(UB);
  UBinv.Invert();
  UBinv /= 2 * M_PI;

  double GonRotx = getParameter("GonRotx");
  double GonRoty = getParameter("GonRoty");
  double GonRotz = getParameter("GonRotz");
  Matrix<double> InvGonRotxMat = RotationMatrixAboutRegAxis(GonRotx, 'x');
  Matrix<double> InvGonRotyMat = RotationMatrixAboutRegAxis(GonRoty, 'y');
  Matrix<double> InvGonRotzMat = RotationMatrixAboutRegAxis(GonRotz, 'z');
  Matrix<double> GonRot = InvGonRotxMat * InvGonRotyMat * InvGonRotzMat;

  InvGonRotxMat.Invert();
  InvGonRotyMat.Invert();
  InvGonRotzMat.Invert();

  std::map<int, Kernel::Matrix<double>> RunNums2GonMatrix;
  getRun2MatMap(Peaks, OptRuns, RunNums2GonMatrix);

  g_log.debug()
      << "----------------------------Derivative------------------------\n";

  V3D samplePosition = instNew->getSample()->getPos();
  IPeak &ppeak = Peaks->getPeak(0);
  double L0 = ppeak.getL1();
  double velocity = (L0 + ppeak.getL2()) / ppeak.getTOF();

  double K =
      2 * M_PI / ppeak.getWavelength() / velocity; // 2pi/lambda = K* velocity
  V3D beamDir = instNew->getBeamDirection();

  size_t paramNums[] = {parameterIndex(std::string("SampleXOffset")),
                        parameterIndex(std::string("SampleYOffset")),
                        parameterIndex(std::string("SampleZOffset"))};

  for (size_t i = 0; i < nData; i += 3) {
    int peakNum = boost::math::iround(xValues[i]);
    IPeak &peak_old = Peaks->getPeak(peakNum);
    Peak peak = createNewPeak(peak_old, instNew, 0, peak_old.getL1());

    int runNum = peak_old.getRunNumber();
    std::string runNumStr = std::to_string(runNum);

    for (int kk = 0; kk < static_cast<int>(nParams()); kk++) {
      out->set(i, kk, 0.0);
      out->set(i + 1, kk, 0.0);
      out->set(i + 2, kk, 0.0);
    }

    double chi, phi, omega;
    size_t chiParamNum, phiParamNum, omegaParamNum;

    size_t N = OptRuns.find("/" + runNumStr);
    if (N < OptRuns.size()) {
      chi = getParameter("chi" + (runNumStr));
      phi = getParameter("phi" + (runNumStr));
      omega = getParameter("omega" + (runNumStr));

      peak.setGoniometerMatrix(GonRot * RunNums2GonMatrix[runNum]);

      chiParamNum = parameterIndex("chi" + (runNumStr));
      phiParamNum = parameterIndex("phi" + (runNumStr));
      omegaParamNum = parameterIndex("omega" + (runNumStr));
    } else {

      Geometry::Goniometer Gon(peak.getGoniometerMatrix());
      std::vector<double> phichiOmega = Gon.getEulerAngles("YZY");
      chi = phichiOmega[1];
      phi = phichiOmega[2];
      omega = phichiOmega[0];
      // peak.setGoniometerMatrix( GonRot*Gon.getR());
      chiParamNum = phiParamNum = omegaParamNum = nParams() + 10;
      peak.setGoniometerMatrix(GonRot * peak.getGoniometerMatrix());
    }
    V3D sampOffsets(getParameter("SampleXOffset"),
                    getParameter("SampleYOffset"),
                    getParameter("SampleZOffset"));
    peak.setSamplePos(peak.getSamplePos() + sampOffsets);
    // NOTE:Use getQLabFrame except for below.
    // For parameters the getGoniometerMatrix should remove GonRot, for derivs
    // wrt GonRot*, wrt chi*,phi*,etc.

    // Deriv wrt chi phi and omega
    if (phiParamNum < nParams()) {
      Matrix<double> chiMatrix = RotationMatrixAboutRegAxis(chi, 'z');
      Matrix<double> phiMatrix = RotationMatrixAboutRegAxis(phi, 'y');
      Matrix<double> omegaMatrix = RotationMatrixAboutRegAxis(omega, 'y');

      Matrix<double> dchiMatrix = DerivRotationMatrixAboutRegAxis(chi, 'z');
      Matrix<double> dphiMatrix = DerivRotationMatrixAboutRegAxis(phi, 'y');
      Matrix<double> domegaMatrix = DerivRotationMatrixAboutRegAxis(omega, 'y');

//.........这里部分代码省略.........

//.........这里部分代码省略.........

        // Distance between this box and a box we already put in.
        coord_t distSquared = 0.0;
        for (size_t d=0; d<nd; d++)
        {
          coord_t dist = otherCenter[d] - boxCenter[d];
          distSquared += (dist * dist);
        }

        // Reject this box if it is too close to another previously found box.
        if (distSquared < peakRadiusSquared)
        {
          badBox = true;
          break;
        }
      }

      // The box was not rejected for another reason.
      if (!badBox)
      {
        if (numBoxesFound++ >= MaxPeaks)
        {
          g_log.notice() << "Number of peaks found exceeded the limit of " << MaxPeaks << ". Stopping peak finding." << std::endl;
          break;
        }

        peakBoxes.push_back(box);
        g_log.information() << "Found box at ";
        for (size_t d=0; d<nd; d++)
          g_log.information() << (d>0?",":"") << boxCenter[d];
        g_log.information() << "; Density = " << density << std::endl;
        // Report progres for each box found.
        prog->report("Finding Peaks");
      }
    }

    prog->resetNumSteps(numBoxesFound, 0.95, 1.0);

    // Copy the instrument, sample, run to the peaks workspace.
    peakWS->copyExperimentInfoFrom(ei.get());

    // --- Convert the "boxes" to peaks ----
    for (typename std::vector<boxPtr>::iterator it3=peakBoxes.begin(); it3 != peakBoxes.end(); it3++)
    {
      // The center of the box = Q in the lab frame
      boxPtr box = *it3;
#ifndef MDBOX_TRACK_CENTROID
      coord_t boxCenter[nd];
      box->calculateCentroid(boxCenter);
#else
      const coord_t * boxCenter = box->getCentroid();
#endif

      V3D Q(boxCenter[0], boxCenter[1], boxCenter[2]);

      // Create a peak and add it
      // Empty starting peak.
      Peak p;
      try
      {
        if (dimType == QLAB)
        {
          // Build using the Q-lab-frame constructor
          p = Peak(inst, Q);
          // Save gonio matrix for later
          p.setGoniometerMatrix(goniometer);
        }
        else if (dimType == QSAMPLE)
        {
          // Build using the Q-sample-frame constructor
          p = Peak(inst, Q, goniometer);
        }
      }
      catch (std::exception &e)
      {
        g_log.notice() << "Error creating peak at " << Q << " because of '" << e.what() << "'. Peak will be skipped." << std::endl;
        continue;
      }

      try
      { // Look for a detector
        p.findDetector();
      }
      catch (...)
      { /* Ignore errors in ray-tracer TODO: Handle for WISH data later */ }

      // The "bin count" used will be the box density.
      p.setBinCount( box->getSignalNormalized() * densityScalingFactor);

      // Save the run number found before.
      p.setRunNumber(runNumber);

      peakWS->addPeak(p);

      // Report progres for each box found.
      prog->report("Adding Peaks");

    } // for each box found

  }

/**
 * Calculates the h,k, and l offsets from an integer for (some of )the peaks,
 *given the parameter values.
 *
 * @param out  For each peak there are 3 consecutive elements in this array. The
 *first is for the h offset from an
 *             integer, the second is the k offset and the 3rd is the l offset
 * @param xValues  xValues give the index in the PeaksWorkspace for the peak.
 *For each peak considered there are
 *              three consecutive entries all with the same index
 * @param nData The size of the xValues and out arrays
 */
void PeakHKLErrors::function1D(double *out, const double *xValues,
                               const size_t nData) const {
  PeaksWorkspace_sptr Peaks =
      AnalysisDataService::Instance().retrieveWS<PeaksWorkspace>(
          PeakWorkspaceName);

  boost::shared_ptr<Geometry::Instrument> instNew = getNewInstrument(Peaks);

  if (!Peaks)
    throw std::invalid_argument("Peaks not stored under the name " +
                                PeakWorkspaceName);

  std::map<int, Mantid::Kernel::Matrix<double>> RunNum2GonMatrixMap;
  getRun2MatMap(Peaks, OptRuns, RunNum2GonMatrixMap);
  const DblMatrix &UBx = Peaks->sample().getOrientedLattice().getUB();

  DblMatrix UBinv(UBx);
  UBinv.Invert();
  UBinv /= (2 * M_PI);

  double GonRotx = getParameter("GonRotx");
  double GonRoty = getParameter("GonRoty");
  double GonRotz = getParameter("GonRotz");
  Matrix<double> GonRot = RotationMatrixAboutRegAxis(GonRotx, 'x') *
                          RotationMatrixAboutRegAxis(GonRoty, 'y') *
                          RotationMatrixAboutRegAxis(GonRotz, 'z');

  double ChiSqTot = 0.0;
  for (size_t i = 0; i < nData; i += 3) {
    int peakNum = boost::math::iround(xValues[i]);
    IPeak &peak_old = Peaks->getPeak(peakNum);

    int runNum = peak_old.getRunNumber();
    std::string runNumStr = std::to_string(runNum);
    Peak peak = createNewPeak(peak_old, instNew, 0, peak_old.getL1());

    size_t N = OptRuns.find("/" + runNumStr + "/");
    if (N < OptRuns.size()) {
      peak.setGoniometerMatrix(GonRot * RunNum2GonMatrixMap[runNum]);

    } else {
      peak.setGoniometerMatrix(GonRot * peak.getGoniometerMatrix());
    }
    V3D sampOffsets(getParameter("SampleXOffset"),
                    getParameter("SampleYOffset"),
                    getParameter("SampleZOffset"));
    peak.setSamplePos(peak.getSamplePos() + sampOffsets);

    V3D hkl = UBinv * peak.getQSampleFrame();

    for (int k = 0; k < 3; k++) {
      double d1 = hkl[k] - floor(hkl[k]);
      if (d1 > .5)
        d1 = d1 - 1;
      if (d1 < -.5)
        d1 = d1 + 1;

      out[i + k] = d1;
      ChiSqTot += d1 * d1;
    }
  }

  g_log.debug() << "------------------------Function---------------------------"
                   "--------------------\n";
  for (size_t p = 0; p < nParams(); p++) {
    g_log.debug() << parameterName(p) << "(" << getParameter(p) << "),";
    if ((p + 1) % 6 == 0)
      g_log.debug() << '\n';
  }
  g_log.debug() << '\n';
  g_log.debug() << "Off constraints=";
  for (size_t p = 0; p < nParams(); p++) {
    IConstraint *constr = getConstraint(p);
    if (constr)
      if ((constr->check() > 0))
        g_log.debug() << "(" << parameterName(p) << "=" << constr->check()
                      << ");";
  }
  g_log.debug() << '\n';

  g_log.debug() << "    Chi**2 = " << ChiSqTot << "     nData = " << nData
                << '\n';
}