31 #include "Eigen/Sparse"
45 void PhotometryFit::leastSquareDerivativesMeasurement(CcdImage
const &ccdImage, TripletList &tripletList,
46 Eigen::VectorXd &grad,
47 MeasuredStarList
const *measuredStarList)
const {
56 if (measuredStarList) assert(&(measuredStarList->front()->getCcdImage()) == &ccdImage);
58 std::size_t nparModel = (_fittingModel) ? _photometryModel->getNpar(ccdImage) : 0;
62 if (_fittingModel) _photometryModel->getMappingIndices(ccdImage, indices);
64 Eigen::VectorXd H(nparTotal);
66 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
67 const MeasuredStarList &catalog = (measuredStarList) ? *measuredStarList : ccdImage.getCatalogForFit();
69 for (
auto const &measuredStar : catalog) {
70 if (!measuredStar->isValid())
continue;
73 double residual = _photometryModel->computeResidual(ccdImage, *measuredStar);
74 double inverseSigma = 1.0 / _photometryModel->transformError(ccdImage, *measuredStar);
75 double W =
std::pow(inverseSigma, 2);
78 _photometryModel->computeParameterDerivatives(*measuredStar, ccdImage, H);
80 Eigen::Index l = indices[k];
81 tripletList.addTriplet(l, kTriplets, H[k] * inverseSigma);
82 grad[l] += H[k] * W * residual;
86 Eigen::Index index = measuredStar->getFittedStar()->getIndexInMatrix();
88 tripletList.addTriplet(index, kTriplets, -1.0 * inverseSigma);
89 grad[index] += -1.0 * W * residual;
94 tripletList.setNextFreeIndex(kTriplets);
97 void PhotometryFit::leastSquareDerivativesReference(FittedStarList
const &fittedStarList,
98 TripletList &tripletList, Eigen::VectorXd &grad)
const {
105 if (!_fittingFluxes)
return;
109 Eigen::Index kTriplets = tripletList.getNextFreeIndex();
111 for (
auto const &fittedStar : fittedStarList) {
112 auto refStar = fittedStar->getRefStar();
113 if (refStar ==
nullptr)
continue;
121 double inverseSigma = 1.0 / _photometryModel->getRefError(*refStar);
123 double residual = _photometryModel->computeRefResidual(*fittedStar, *refStar);
125 Eigen::Index index = fittedStar->getIndexInMatrix();
127 tripletList.addTriplet(index, kTriplets, 1.0 * inverseSigma);
128 grad(index) += 1.0 *
std::pow(inverseSigma, 2) * residual;
131 tripletList.setNextFreeIndex(kTriplets);
134 void PhotometryFit::accumulateStatImageList(
CcdImageList const &ccdImageList, Chi2Accumulator &accum)
const {
139 for (
auto const &ccdImage : ccdImageList) {
140 auto &catalog = ccdImage->getCatalogForFit();
142 for (
auto const &measuredStar : catalog) {
143 if (!measuredStar->isValid())
continue;
144 double sigma = _photometryModel->transformError(*ccdImage, *measuredStar);
145 double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
148 accum.addEntry(chi2Val, 1, measuredStar);
153 void PhotometryFit::accumulateStatRefStars(Chi2Accumulator &accum)
const {
159 FittedStarList &fittedStarList =
_associations->fittedStarList;
160 for (
auto const &fittedStar : fittedStarList) {
161 auto refStar = fittedStar->getRefStar();
162 if (refStar ==
nullptr)
continue;
163 double sigma = _photometryModel->getRefError(*refStar);
164 double residual = _photometryModel->computeRefResidual(*fittedStar, *refStar);
166 accum.addEntry(chi2, 1, fittedStar);
173 void PhotometryFit::getIndicesOfMeasuredStar(MeasuredStar
const &measuredStar,
177 _photometryModel->getMappingIndices(measuredStar.getCcdImage(), indices);
179 if (_fittingFluxes) {
181 Eigen::Index fsIndex = fs->getIndexInMatrix();
190 _fittingFluxes = (
_whatToFit.
find(
"Fluxes") != std::string::npos);
193 _nParModel = (_fittingModel) ? _photometryModel->assignIndices(whatToFit, 0) : 0;
196 if (_fittingFluxes) {
202 fittedStar->setIndexInMatrix(ipar);
206 _nParFluxes = ipar - _nParModel;
209 "nParameters total: " <<
_nParTot <<
" model: " << _nParModel <<
" fluxes: " << _nParFluxes);
215 "PhotometryFit::offsetParams : the provided vector length is not compatible with "
216 "the current whatToFit setting");
217 if (_fittingModel) _photometryModel->offsetParams(delta);
219 if (_fittingFluxes) {
224 Eigen::Index index = fittedStar->getIndexInMatrix();
225 _photometryModel->offsetFittedStar(*fittedStar, delta(index));
234 ofile <<
"#id" << separator <<
"xccd" << separator <<
"yccd" << separator;
235 ofile <<
"mag" << separator <<
"instMag" << separator <<
"instMagErr" << separator;
236 ofile <<
"instFlux" << separator <<
"instFluxErr" << separator;
237 ofile <<
"inputFlux" << separator <<
"inputFluxErr" << separator;
238 ofile <<
"transformedFlux" << separator <<
"transformedFluxErr" << separator;
239 ofile <<
"fittedFlux" << separator;
240 ofile <<
"mjd" << separator <<
"color" << separator;
241 ofile <<
"fsindex" << separator;
242 ofile <<
"ra" << separator <<
"dec" << separator;
243 ofile <<
"chi2" << separator <<
"nm" << separator;
244 ofile <<
"chip" << separator <<
"visit" << separator <<
std::endl;
246 ofile <<
"#id in source catalog" << separator <<
"coordinates in CCD" << separator << separator;
247 ofile <<
"fitted magnitude" << separator <<
"measured magnitude" << separator
248 <<
"measured magnitude error" << separator;
249 ofile <<
"measured instrumental flux (ADU)" << separator <<
"measured instrument flux error" << separator;
250 ofile <<
"measured flux (nJy)" << separator <<
"measured flux error" << separator;
251 ofile << separator << separator;
252 ofile <<
"fitted flux (nJy)" << separator;
253 ofile <<
"modified Julian date of the measurement" << separator <<
"currently unused" << separator;
254 ofile <<
"unique index of the fittedStar" << separator;
255 ofile <<
"on-sky position of fitted star" << separator << separator;
256 ofile <<
"contribution to Chi2 (1 dof)" << separator <<
"number of measurements of this FittedStar"
258 ofile <<
"chip id" << separator <<
"visit id" <<
std::endl;
261 for (
auto const &ccdImage : ccdImageList) {
263 for (
auto const &measuredStar : cat) {
264 if (!measuredStar->isValid())
continue;
266 double instFluxErr = _photometryModel->tweakFluxError(*measuredStar);
267 double flux = _photometryModel->transform(*ccdImage, *measuredStar);
268 double fluxErr = _photometryModel->transformError(*ccdImage, *measuredStar);
269 double jd = ccdImage->getMjd();
271 double residual = _photometryModel->computeResidual(*ccdImage, *measuredStar);
272 double chi2Val =
std::pow(residual / fluxErr, 2);
275 ofile << measuredStar->getId() << separator << measuredStar->x << separator << measuredStar->y
277 ofile << fittedStar->getMag() << separator << measuredStar->getInstMag() << separator
278 << measuredStar->getInstMagErr() << separator;
279 ofile << measuredStar->getInstFlux() << separator << instFluxErr << separator;
280 ofile << measuredStar->getFlux() << separator << measuredStar->getFluxErr() << separator;
281 ofile <<
flux << separator << fluxErr << separator << fittedStar->getFlux() << separator;
282 ofile << jd << separator << fittedStar->color << separator;
283 ofile << fittedStar->getIndexInMatrix() << separator;
284 ofile << fittedStar->x << separator << fittedStar->y << separator;
285 ofile << chi2Val << separator << fittedStar->getMeasurementCount() << separator;
286 ofile << ccdImage->getCcdId() << separator << ccdImage->getVisit() <<
std::endl;
295 ofile <<
"#ra" << separator <<
"dec " << separator;
296 ofile <<
"mag" << separator <<
"color" << separator;
297 ofile <<
"refFlux" << separator <<
"refFluxErr" << separator;
298 ofile <<
"fittedFlux" << separator <<
"fittedFluxErr" << separator;
299 ofile <<
"fsindex" << separator <<
"chi2" << separator <<
"nm" <<
std::endl;
301 ofile <<
"#coordinates of fittedStar" << separator << separator;
302 ofile <<
"magnitude" << separator <<
"currently unused" << separator;
303 ofile <<
"refStar flux (nJy)" << separator <<
"refStar fluxErr" << separator;
304 ofile <<
"fittedStar flux (nJy)" << separator <<
"fittedStar fluxErr" << separator;
305 ofile <<
"unique index of the fittedStar" << separator <<
"refStar contribution to Chi2 (1 dof)"
306 << separator <<
"number of measurements of this FittedStar" <<
std::endl;
310 for (
auto const &fittedStar : fittedStarList) {
311 const RefStar *refStar = fittedStar->getRefStar();
312 if (refStar ==
nullptr)
continue;
317 ofile << fittedStar->x << separator << fittedStar->y << separator;
318 ofile << fittedStar->getMag() << separator << fittedStar->color << separator;
320 ofile << fittedStar->getFlux() << separator << fittedStar->getFluxErr() << separator;
321 ofile << fittedStar->getIndexInMatrix() << separator << chi2 << separator
322 << fittedStar->getMeasurementCount() <<
std::endl;