- Timestamp:
- Jul 1, 2016, 10:11:46 AM (10 years ago)
- File:
-
- 1 edited
-
trunk/Ohana/src/relphot/src/setMrelCatalog.c (modified) (1 diff)
Legend:
- Unmodified
- Added
- Removed
-
trunk/Ohana/src/relphot/src/setMrelCatalog.c
r39609 r39620 963 963 } 964 964 965 # if (0) 966 // outlier warp measurements are driving bad mean values. 967 int magStatsByRankingClipped (StatDataSet *dataset, StatType *stats) { 968 969 liststats_init (stats); 970 971 if (dataset->Nlist == 0) return 0; 972 973 // we have a list of measurements for this Nsec value, along with errors, weights, and their rank 974 // first, sort by the rank (increasing) 975 sort_StatDataSet (dataset); 976 977 // find the values with the same minimum rank: 978 int Nranking = 0; 979 int minRank = dataset->ranking[0]; // MIN (5, result->psfData[Nsec].ranking[0]); -- only allow rank 5 or < 980 981 int i; 982 for (i = 0; (i < dataset->Nlist) && (dataset->ranking[i] == minRank); i++, Nranking++); 983 984 // we want to measure the weighted average, but first we want to clip points which are > 3 sigma from the median 985 // since we are only allowed to clip at most 25% of the points, we need to have at least 4 points to do this clipping 986 if (Nranking > 3) { 987 // first, measure the median of the flxlist: 988 ALLOCATE (fluxes, double, Nranking); 989 for (i = 0; i < Nranking; i++) { 990 fluxes[i] = dataset->flxlist[i]; 991 } 992 dsort (fluxes, Nranking); 993 int Nmidpoint = (int)(0.5*Nranking); // if Nranking is odd (e.g., 5), Nmidpoint is central bin (e.g., 2); if Nranking is even (e.g., 6), Nmidpoint is central bin + 1 (e.g., 3); 994 double medianFlux = (Nranking % 2) ? fluxes[Nmidpoint] : 0.5*(fluxes[Nmidpoint] + fluxes[Nmidpoint - 1]); 995 } 996 997 // get stats on only the first Nranking entries 998 liststats (dataset->flxlist, dataset->errlist, dataset->wgtlist, Nranking, stats); 999 1000 return (Nranking); 1001 } 1002 # endif 1003 965 1004 int markMeasureByRanking (StatDataSet *dataset, Measure *measure, int minrank, DVOMeasureFlags flags) { 966 1005 int i;
Note:
See TracChangeset
for help on using the changeset viewer.
