72 LOGLS_DEBUG(_log,
"Catalog " << ccdImage->getName() <<
" has " << ccdImage->getWholeCatalog().size()
80 centers.
push_back(ccdImage->getBoresightRaDec());
88 _commonTangentPoint =
Point(commonTangentPoint.getX(), commonTangentPoint.getY());
89 for (
auto &ccdImage :
ccdImageList) ccdImage->setCommonTangentPoint(_commonTangentPoint);
94 Frame tangentPlaneFrame;
97 Frame CTPFrame = ccdImage->getPixelToCommonTangentPlane()->apply(ccdImage->getImageFrame(),
false);
98 if (tangentPlaneFrame.
getArea() == 0)
99 tangentPlaneFrame = CTPFrame;
101 tangentPlaneFrame += CTPFrame;
106 TanPixelToRaDec commonTangentPlaneToRaDec(identity, _commonTangentPoint);
107 Frame raDecFrame = commonTangentPlaneToRaDec.
apply(tangentPlaneFrame,
false);
122 const bool enlargeFittedList) {
128 item->clearBeforeAssoc();
138 ccdImage->resetCatalogForFit();
144 Frame ccdImageFrameCPT = toCommonTangentPlane->apply(ccdImage->getImageFrame(),
false);
145 ccdImageFrameCPT = ccdImageFrameCPT.
rescale(1.10);
151 if (ccdImageFrameCPT.
inFrame(*fittedStar)) {
158 toCommonTangentPlane.
get(), matchCutInArcSec / 3600.);
161 LOGLS_DEBUG(_log,
"Measured-to-Fitted matches before removing ambiguities " << starMatchList->size());
162 starMatchList->removeAmbiguities(*toCommonTangentPlane);
163 LOGLS_DEBUG(_log,
"Measured-to-Fitted matches after removing ambiguities " << starMatchList->size());
167 int matchedCount = 0;
168 for (
auto const &starMatch : *starMatchList) {
169 auto bs = starMatch.s1;
170 auto ms_const = std::dynamic_pointer_cast<const MeasuredStar>(bs);
171 auto ms = std::const_pointer_cast<MeasuredStar>(ms_const);
172 auto bs2 = starMatch.s2;
173 auto fs_const = std::dynamic_pointer_cast<const FittedStar>(bs2);
174 auto fs = std::const_pointer_cast<FittedStar>(fs_const);
175 ms->setFittedStar(fs);
178 LOGLS_INFO(_log,
"Matched " << matchedCount <<
" objects in " << ccdImage->getName());
181 int unMatchedCount = 0;
182 for (
auto const &mstar : catalog) {
184 if (mstar->getFittedStar())
continue;
185 if (enlargeFittedList) {
186 auto fs = std::make_shared<FittedStar>(*mstar);
188 toCommonTangentPlane->transformPosAndErrors(*fs, *fs);
190 mstar->setFittedStar(fs);
194 LOGLS_INFO(_log,
"Unmatched objects: " << unMatchedCount);
206 std::string const &fluxField,
float refCoordinateErr,
207 bool rejectBadFluxes) {
208 if (refCat.size() == 0) {
210 " reference catalog is empty : stop here "));
218 raErrKey = refCat.getSchema()[
"coord_raErr"];
219 decErrKey = refCat.getSchema()[
"coord_decErr"];
222 auto fluxKey = refCat.getSchema().
find<
double>(fluxField).
key;
226 fluxErrKey = refCat.getSchema().
find<
double>(fluxField +
"Err").
key;
229 << fluxField <<
"Err"
230 <<
") not found in reference catalog. Not using ref flux errors.");
234 for (
size_t i = 0; i < refCat.size(); i++) {
235 auto const &record = refCat.get(i);
237 auto coord = record->get(coordKey);
238 double flux = record->get(fluxKey);
240 if (fluxErrKey.isValid()) {
241 fluxErr = record->get(fluxErrKey);
247 auto star = std::make_shared<RefStar>(ra,
dec,
flux, fluxErr);
255 star->vx =
std::pow(refCoordinateErr / 1000. / 3600. /
std::cos(coord.getLatitude()), 2);
256 star->vy =
std::pow(refCoordinateErr / 1000. / 3600., 2);
268 TanRaDecToPixel raDecToCommonTangentPlane(identity, _commonTangentPoint);
270 associateRefStars(matchCut.
asArcseconds(), &raDecToCommonTangentPlane);
277 matchCutInArcSec / 3600.);
279 LOGLS_DEBUG(_log,
"Refcat matches before removing ambiguities " << starMatchList->size());
280 starMatchList->removeAmbiguities(*
transform);
281 LOGLS_DEBUG(_log,
"Refcat matches after removing ambiguities " << starMatchList->size());
284 for (
auto const &starMatch : *starMatchList) {
288 const BaseStar &bs2 = *starMatch.s2;
296 "Associated " << starMatchList->size() <<
" reference stars among " <<
refStarList.
size());
300 selectFittedStars(minMeasurements);
301 normalizeFittedStars();
304 void Associations::selectFittedStars(
int minMeasurements) {
307 int totalMeasured = 0, validMeasured = 0;
319 if (fittedStar ==
nullptr) {
326 if (!fittedStar->getRefStar() && fittedStar->getMeasurementCount() < minMeasurements) {
327 fittedStar->getMeasurementCount()--;
328 mi = catalog.
erase(mi);
338 if ((*fi)->getMeasurementCount() == 0) {
346 LOGLS_INFO(_log,
"Total, valid number of Measured stars: " << totalMeasured <<
", " << validMeasured);
349 void Associations::normalizeFittedStars()
const {
354 fittedStar->setFlux(0.0);
355 fittedStar->getMag() = 0.0;
361 MeasuredStarList &catalog = ccdImage->getCatalogForFit();
362 for (
auto &mi : catalog) {
363 auto fittedStar = mi->getFittedStar();
364 if (fittedStar ==
nullptr)
366 pex::exceptions::RuntimeError,
367 "All measuredStars must have a fittedStar: did you call selectFittedStars()?"));
368 auto point = toCommonTangentPlane->apply(*mi);
369 fittedStar->x += point.x;
370 fittedStar->y += point.y;
371 fittedStar->getFlux() += mi->getFlux();
376 auto measurementCount = fi->getMeasurementCount();
377 fi->x /= measurementCount;
378 fi->y /= measurementCount;
379 fi->getFlux() /= measurementCount;
384 void Associations::assignMags() {
386 MeasuredStarList &catalog = ccdImage->getCatalogForFit();
387 for (
auto const &mstar : catalog) {
389 if (!fstar)
continue;
400 "DeprojectFittedStars: Fitted stars are already in sidereal coordinates, nothing done ");
411 return item->getCatalogForFit().size() > 0;
418 if ((fittedStar !=
nullptr) & (fittedStar->getRefStar() !=
nullptr))
count++;
424 void Associations::collectMCStars(
int realization) {
430 string dbimdir = ccdImage.Dir();
431 string mctruth = dbimdir +
"/mc/mctruth.list";
433 if (realization >= 0) {
435 sstrm << dbimdir <<
"/mc/mctruth_" << realization <<
".list";
436 mctruth = sstrm.
str();
439 AstrometryTransformIdentity gti;
443 DicStarList mctruthlist(mctruth);
447 for (smI = starMatchList->begin(); smI != starMatchList->end(); smI++) {
448 StarMatch &sm = *smI;
449 BaseStar *bs = sm.s1;
450 MeasuredStar *mstar =
dynamic_cast<MeasuredStar *
>(bs);
452 DicStar *dstar =
dynamic_cast<DicStar *
>(bs);
454 mcstar->GetMCInfo().iflux = dstar->getval(
"iflux");
455 mcstar->GetMCInfo().tflux = dstar->getval(
"sflux");
462 LOGLS_FATAL(_log,
"CollectMCStars Unable to match MCTruth w/ catalog!");
467 double matchCutArcSec) {
469 size_t pos_minus = color.
find(
'-');
470 bool compute_diff = (pos_minus != string::npos);
472 c1 = color.
substr(0, pos_minus);
473 if (compute_diff) c2 = color.
substr(pos_minus + 1, string::npos);
474 DicStarList cList(dicStarListName);
475 if (!cList.HasKey(c1))
476 throw(GastroException(
"Associations::SetFittedstarColors : " + dicStarListName +
477 " misses a key named \"" + c1 +
"\""));
478 if (compute_diff && !cList.HasKey(c2))
479 throw(GastroException(
"Associations::SetFittedstarColors : " + dicStarListName +
480 " misses a key named \"" + c2 +
"\""));
484 AstrometryTransformIdentity
id;
487 AstrometryTransform *id_or_proj = &proj;
493 id_or_proj, matchCutArcSec / 3600);
496 <<
" FittedStars to color catalog");
498 for (
auto i = starMatchList->begin(); i != starMatchList->end(); ++i) {
499 BaseStar *s1 = i->s1;
500 FittedStar *fs =
dynamic_cast<FittedStar *
>(s1);
501 BaseStar *s2 = i->s2;
502 const TStar *ts =
dynamic_cast<const TStar *
>(s2);
503 const DicStar *ds =
dynamic_cast<const DicStar *
>(ts->get_original());
504 fs->color = ds->getval(c1);
505 if (compute_diff) fs->color -= ds->getval(c2);