Commits

Kumar Golap authored d7ab8c0b5e5 Merge 
Merge branch 'master' into feature/CAS-9004
49add4f3d7f
master

code/alma/apps/asdm2MS/asdm2MS.cc

Modified
52 52 #include "ASDM2MSFiller.h"
53 53 ​
54 54 #include "measures/Measures/Stokes.h"
55 55 #include "measures/Measures/MFrequency.h"
56 56 using namespace casacore;
57 57 #include <tables/Tables/Table.h>
58 58 #include <tables/Tables/PlainTable.h>
59 59 #include <tables/Tables/TableCopy.h>
60 60 #include <tables/Tables/TableInfo.h>
61 61 #include <casa/Arrays/MatrixMath.h>
62 -​
63 -​
62 +#include <casa/BasicMath/Math.h>
64 63 ​
65 64 #include "CBasebandName.h"
66 65 #include "CCalibrationDevice.h"
67 66 using namespace CalibrationDeviceMod;
68 67 #include "CFrequencyReferenceCode.h"
69 68 #include "CPolarizationType.h"
70 69 #include "CProcessorSubType.h"
71 70 #include "CProcessorType.h"
72 71 #include "CScanIntent.h"
73 72 #include "CSubscanIntent.h"
1182 1181 dprime[0] = d[0] / b[0];
1183 1182 for (unsigned int i = 1; i < n; i++)
1184 1183 dprime[i] = (d[i] - dprime[i-1] * a[i]) / (b[i] - cprime[i-1] * a[i]);
1185 1184 LOG("dprime computed");
1186 1185 ​
1187 1186 x.clear(); x.resize(n);
1188 1187 x[n-1] = dprime[n-1];
1189 1188 for (int i = n-2; i >=0; i--) {
1190 1189 x[i] = dprime[i] - cprime[i] * x[i+1];
1191 1190 }
1192 - cout << endl;
1193 1191 LOGEXIT("solveTridiagonalSystem");
1194 1192 }
1195 1193 ​
1196 1194 void linearInterpCoeff(uint32_t npoints,
1197 1195 const vector<double>& time_v,
1198 1196 const vector<double>& k_v,
1199 1197 vector<vector<double> >& coeff_vv) {
1200 1198 LOGENTER("linearInterpCoeff");
1201 1199 coeff_vv.clear();
1202 1200 coeff_vv.resize(npoints-1);
2559 2557 string bdfName; // The path to the BDF which contains the data,
2560 2558 int32_t index; // The (0-based) index of the row in the Main ASDM table,
2561 2559 bool uncorrected; // true if this row contains uncorrected data,
2562 2560 bool corrected; // true if this row contains corrected data,
2563 2561 } MainRowCUStruct;
2564 2562 ​
2565 2563 void fillMainLazily(const string& dsName,
2566 2564 ASDM* ds_p,
2567 2565 std::map<int, std::set<int> >& selected_eb_scan_m,
2568 2566 std::map<unsigned int , double>& effectiveBwPerDD_m,
2569 - Enum<CorrelationMode> e_query_cm) {
2567 + Enum<CorrelationMode> e_query_cm,
2568 + bool checkDupInts) {
2570 2569 ​
2571 2570 LOGENTER("fillMainLazily");
2572 2571 ​
2573 2572 ostringstream oss;
2574 2573 ​
2575 2574 MainTable& mainT = ds_p->getMain();
2576 2575 ConfigDescriptionTable& cfgDscT = ds_p->getConfigDescription();
2577 2576 ProcessorTable& procT = ds_p->getProcessor();
2578 2577 ​
2579 2578 MainRowCUStruct mRCU_s;
2690 2689 info(infostream.str());
2691 2690 ​
2692 2691 // Now traverse the BDFs :
2693 2692 // * to write the indexes for asdmstman
2694 2693 // * to populate all the columns other than the DATA's one in the non lazy way.
2695 2694 //
2696 2695 ​
2697 2696 uInt lastMSNUrows = 0;
2698 2697 uInt lastMSNCrows = 0;
2699 2698 ​
2699 + // used in checking for duplicate integrations in the WVR (Radiometer) case
2700 + // This holds the most recent last integration time for each configDescriptionId - but only for Radiometer data.
2701 + map<Tag, double> lastTimeMap;
2702 +​
2700 2703 try {
2701 - for (vector<MainRowCUStruct>::iterator iter=mRCU_s_v.begin(); iter!=mRCU_s_v.end(); iter++) {
2704 + unsigned int mainRowIndex;
2705 + vector<MainRowCUStruct>::iterator iter;
2706 + for (iter=mRCU_s_v.begin(), mainRowIndex=0; iter!=mRCU_s_v.end(); iter++, mainRowIndex++) {
2702 2707 MainRow* mR_p = iter->mR_p;
2703 2708 ​
2704 2709 SDMDataObjectStreamReader sdosr;
2705 2710 sdosr.open(iter->bdfName);
2706 2711 LOG("Processing " + iter->bdfName);
2707 2712 unsigned int numberOfAntennas = sdosr.numAntenna();
2708 2713 unsigned int numberOfBaselines = numberOfAntennas * (numberOfAntennas - 1) / 2 ;
2709 2714 ProcessorType processorType = sdosr.processorType();
2710 2715 CorrelationMode correlationMode = sdosr.correlationMode();
2711 2716 const SDMDataObject::DataStruct& dataStruct = sdosr.dataStruct();
2929 2934 ​
2930 2935 //
2931 2936 // Everything is contained in *one* SDMDataSubset.
2932 2937 //
2933 2938 const SDMDataSubset& sdmDataSubset = sdosr.getSubset();
2934 2939 ​
2935 2940 int64_t deltaTime = sdmDataSubset.interval() / sdosr.numTime();
2936 2941 int64_t startTime = (int64_t)sdmDataSubset.time() - (int64_t)sdmDataSubset.interval()/2LL + deltaTime/2LL;
2937 2942 double interval = deltaTime / 1000000000.0;
2938 2943
2944 + // should the first integration be skipped? Any actual skipping happens later.
2945 + bool skipFirstIntegration = checkDupInts && lastTimeMap[mR_p->getConfigDescriptionId()] == ArrayTime(startTime).getMJD();
2946 + if (debug && skipFirstIntegration) {
2947 + cout << "Duplicate time seen in Row : " << mainRowIndex
2948 + << " cdId : " << mR_p->getConfigDescriptionId()
2949 + << " " << mR_p->getDataUID().getEntityId().toString()
2950 + << " numTime : " << sdosr.numTime()
2951 + << " num MS rows : " << dataDescriptionIds.size()*sdosr.numTime()*antennaIds.size()
2952 + << endl;
2953 + }
2954 + lastTimeMap[mR_p->getConfigDescriptionId()] = ArrayTime(startTime+(sdosr.numTime()-1)*deltaTime).getMJD();
2955 +
2939 2956 for (unsigned int iDD = 0; iDD < dataDescriptionIds.size(); iDD++) {
2940 2957 //
2941 2958 // Prepare a pair<int, int> to transport the shape of some cells
2942 2959 //
2943 2960 pair<int,int> nChanNPol = make_pair<int, int>(numberOfChannels_v[iDD],
2944 2961 numberOfSDPolarizations_v[iDD]);
2945 2962 ​
2946 2963 //
2947 2964 // Compute weight and sigma which depend on the data description id and on the interval
2948 2965 //
2949 2966 double weight = 1.0 * effectiveBwPerDD_m[dataDescriptionIdx2Idx[dataDescriptionIds[iDD].getTagValue()]] * interval;
2950 2967 weight = (weight == 0.0) ? 1.0 : weight;
2951 2968 double sigma = 1./sqrt(weight);
2952 2969 ​
2953 2970 for (unsigned int itime = 0; itime < sdosr.numTime(); itime++) {
2971 + if (skipFirstIntegration && itime==0) continue;
2954 2972 for (unsigned int iA = 0; iA < antennaIds.size(); iA++) {
2955 2973 antenna1_vv[iDD].push_back(antennaIds[iA].getTagValue());
2956 2974 antenna2_vv[iDD].push_back(antennaIds[iA].getTagValue());
2957 2975 dataDescId_vv[iDD].push_back(dataDescriptionIdx2Idx[dataDescriptionIds[iDD].getTagValue()]);
2958 2976 exposure_vv[iDD].push_back(interval);
2959 2977 interval_vv[iDD].push_back(interval);
2960 2978 time_vv[iDD].push_back(ArrayTime(startTime + itime * deltaTime).getMJD() * 86400.0);
2961 2979 feed1_vv[iDD].push_back(feedIds[iA]);
2962 2980 feed2_vv[iDD].push_back(feedIds[iA]);
2963 2981 flagRow_vv[iDD].push_back(false);
2991 3009 numberOfChannels_v[iDD],
2992 3010 numberOfSDPolarizations_v[iDD],
2993 3011 stepSDBl, //numberOfSpectralWindows * numberOfChannels * numberOfPolarizations,
2994 3012 iDD, // this will be used as an index in the seq of windows in the BDFs
2995 3013 autoScaleFactors,
2996 3014 sdmDataSubset.autoDataPosition() + itime * numberOfAntennas * stepSDBl * sizeof(AUTODATATYPE),
2997 3015 spwSDOffset_v[iDD]);
2998 3016 }
2999 3017 }
3000 3018 }
3001 -
3019 +​
3002 3020 // If we have uncorrected data (which is in practice always the case I think) and want to output those then populate the asdmindex of uncorrected data MS.
3003 3021 if (iter->uncorrected and produceUncorrected)
3004 3022 bdf2AsdmStManIndexU.dumpAutoCross();
3005 3023 ​
3006 3024 // If we have corrected data and want to output those then populate the asdmindex of corrected data MS.
3007 3025 if (iter->corrected and produceCorrected)
3008 3026 bdf2AsdmStManIndexC.dumpAutoCross();
3009 3027 ​
3010 3028 //
3011 3029 // It's now time to populate the columns of the MAIN table but the DATA's one.
3012 3030 for (unsigned int iDD = 0; iDD < dataDescriptionIds.size(); iDD++) {
3013 3031 for (map<AtmPhaseCorrection, ASDM2MSFiller*>::iterator msfIter = msFillers.begin();
3014 3032 msfIter != msFillers.end();
3015 3033 ++msfIter) {
3016 - msfIter->second->addData(true, // Yes ! these are complex data.
3017 - time_vv[iDD],
3018 - antenna1_vv[iDD],
3019 - antenna2_vv[iDD],
3020 - feed1_vv[iDD],
3021 - feed2_vv[iDD],
3022 - dataDescId_vv[iDD],
3023 - processorId,
3024 - fieldId,
3025 - interval_vv[iDD],
3026 - exposure_vv[iDD],
3027 - timeCentroid_vv[iDD],
3028 - scanNumber,
3029 - arrayId,
3030 - observationId,
3031 - stateId_vv[iDD],
3032 - nChanNPol_vv[iDD],
3033 - uvw_vv[iDD],
3034 - weight_vv[iDD],
3035 - sigma_vv[iDD]);
3034 + if (time_vv[iDD].size() > 0) {
3035 + msfIter->second->addData(true, // Yes ! these are complex data.
3036 + time_vv[iDD],
3037 + antenna1_vv[iDD],
3038 + antenna2_vv[iDD],
3039 + feed1_vv[iDD],
3040 + feed2_vv[iDD],
3041 + dataDescId_vv[iDD],
3042 + processorId,
3043 + fieldId,
3044 + interval_vv[iDD],
3045 + exposure_vv[iDD],
3046 + timeCentroid_vv[iDD],
3047 + scanNumber,
3048 + arrayId,
3049 + observationId,
3050 + stateId_vv[iDD],
3051 + nChanNPol_vv[iDD],
3052 + uvw_vv[iDD],
3053 + weight_vv[iDD],
3054 + sigma_vv[iDD]);
3055 + }
3036 3056 }
3037 3057 }
3038 3058 }
3039 3059 ​
3040 3060 else if (!sdosr.hasPackedData() && (processorType == CORRELATOR || processorType == RADIOMETER)) {
3041 3061 ​
3062 + // duplicate time skipping is not supported here
3063 +​
3042 3064 //
3043 3065 // Declare some containers required to populate the columns of the MS MAIN table in a non lazy way.
3044 3066 //
3045 3067 // We use vectors of vectors in order to be able to build separate vectors for different data description
3046 3068 // and then output these vectors in the appropriate order.
3047 3069 //
3048 3070 // The cross correlation chapter.
3049 3071 //
3050 3072 vector<vector<int> > cross_antenna1_vv(dataDescriptionIds.size()); // Column ANTENNA1 per Data Description
3051 3073 vector<vector<int> > cross_antenna2_vv(dataDescriptionIds.size()); // Column ANTENNA2 per Data Description
3292 3314 //
3293 3315 // It's now time to populate the columns of the MAIN table but the DATA's one.
3294 3316 // This is done with data descriptions varying the more slowly.
3295 3317 //
3296 3318 for (unsigned int iDD = 0; iDD < dataDescriptionIds.size(); iDD++) {
3297 3319 if (hasAutoData)
3298 3320 for (map<AtmPhaseCorrection, ASDM2MSFiller*>::iterator msfIter = msFillers.begin();
3299 3321 msfIter != msFillers.end();
3300 3322 ++msfIter)
3301 3323 if ((msfIter->first == AP_UNCORRECTED and iter->uncorrected and produceUncorrected) ||
3302 - (msfIter->first == AP_CORRECTED and iter->corrected and produceCorrected))
3324 + (msfIter->first == AP_CORRECTED and iter->corrected and produceCorrected)) {
3303 3325 msfIter->second->addData(true, // Yes ! these are complex data.
3304 3326 auto_time_vv[iDD],
3305 3327 auto_antenna1_vv[iDD],
3306 3328 auto_antenna2_vv[iDD],
3307 3329 auto_feed1_vv[iDD],
3308 3330 auto_feed2_vv[iDD],
3309 3331 auto_dataDescId_vv[iDD],
3310 3332 processorId,
3311 3333 fieldId,
3312 3334 auto_interval_vv[iDD],
3313 3335 auto_exposure_vv[iDD],
3314 3336 auto_timeCentroid_vv[iDD],
3315 3337 scanNumber,
3316 3338 arrayId,
3317 3339 observationId,
3318 3340 auto_stateId_vv[iDD],
3319 3341 auto_nChanNPol_vv[iDD],
3320 3342 auto_uvw_vv[iDD],
3321 3343 auto_weight_vv[iDD],
3322 3344 auto_sigma_vv[iDD]);
3345 + }
3323 3346 if (hasCrossData)
3324 3347 for (map<AtmPhaseCorrection, ASDM2MSFiller*>::iterator msfIter = msFillers.begin();
3325 3348 msfIter != msFillers.end();
3326 3349 ++msfIter)
3327 3350 if ((msfIter->first == AP_UNCORRECTED and iter->uncorrected and produceUncorrected) ||
3328 - (msfIter->first == AP_CORRECTED and iter->corrected and produceCorrected))
3351 + (msfIter->first == AP_CORRECTED and iter->corrected and produceCorrected)) {
3329 3352 msfIter->second->addData(true, // Yes ! these are complex data.
3330 3353 cross_time_vv[iDD],
3331 3354 cross_antenna1_vv[iDD],
3332 3355 cross_antenna2_vv[iDD],
3333 3356 cross_feed1_vv[iDD],
3334 3357 cross_feed2_vv[iDD],
3335 3358 cross_dataDescId_vv[iDD],
3336 3359 processorId,
3337 3360 fieldId,
3338 3361 cross_interval_vv[iDD],
3339 3362 cross_exposure_vv[iDD],
3340 3363 cross_timeCentroid_vv[iDD],
3341 3364 scanNumber,
3342 3365 arrayId,
3343 3366 observationId,
3344 3367 cross_stateId_vv[iDD],
3345 3368 cross_nChanNPol_vv[iDD],
3346 3369 cross_uvw_vv[iDD],
3347 3370 cross_weight_vv[iDD],
3348 - cross_sigma_vv[iDD]);
3371 + cross_sigma_vv[iDD]);
3372 + }
3349 3373 }
3350 3374 }
3351 3375 else
3352 3376 cout << "Processor not supported in lazy mode." << endl;
3353 3377
3354 3378 sdosr.close();
3355 3379 ​
3356 3380 if (iter->uncorrected and produceUncorrected) {
3357 3381 infostream.str("");
3358 3382 infostream << "ASDM Main row #" << iter->index << " produced a total of " << msFillers[AP_UNCORRECTED]->ms()->nrow() - lastMSNUrows << " MS Main rows with uncorrected data." << endl;
3399 3423 ​
3400 3424 /**
3401 3425 * This function fills the MS Main table from an ASDM Main table which refers to correlator data.
3402 3426 *
3403 3427 * given:
3404 3428 * @parameter r_p a pointer to the MainRow being processed.
3405 3429 * @parameter sdmBinData a reference to the SDMBinData containing a lot of information about the binary data being processed. Useful to know the requested ordering of data.
3406 3430 * @parameter uvwCoords a reference to the UVW calculator.
3407 3431 * @parameter complexData a bool which says if the DATA is going to be filled (true) or if it will be the FLOAT_DATA (false).
3408 3432 * @parameter mute if the value of this parameter is false then nothing is written in the MS .
3433 + * @parameter skipFirstTime if the value of this parameter is true, then all rows with time equal to the first time seen will be skipped (not filled).
3409 3434 *
3410 3435 * !!!!! One must be carefull to the fact that fillState must have been called before fillMain. During the execution of fillState , the global vector<int> msStateID
3411 3436 * is filled and will be used by fillMain.
3412 3437 */
3413 3438 void fillMain(
3414 3439 MainRow* r_p,
3415 3440 SDMBinData& sdmBinData,
3416 3441 const VMSData* vmsData_p,
3417 3442 UvwCoords& uvwCoords,
3418 3443 std::map<unsigned int, double>& effectiveBwPerDD_m,
3419 3444 bool complexData,
3420 3445 bool mute,
3421 - bool ac_xc_per_timestamp) {
3446 + bool ac_xc_per_timestamp,
3447 + bool skipFirstTime=false) {
3422 3448
3423 3449 if (debug) cout << "fillMain : entering" << endl;
3424 3450 ​
3425 3451 // ASDM & ds = r_p -> getTable() . getContainer();
3426 3452 ​
3427 3453 // Then populate the Main table.
3428 3454 ComplexDataFilter filter; // To process the case numCorr == 3
3429 3455
3430 3456 if (vmsData_p->v_antennaId1.size() == 0) {
3431 3457 infostream.str("");
3432 3458 infostream << "No MS data produced for the current row." << endl;
3433 3459 info(infostream.str());
3434 3460 return;
3435 3461 }
3436 3462 ​
3437 -
3463 + // msRowReIndex_v maps from location in the vmsData_p vectors to the output MS row
3464 + // this may be reorderd and the first time integration may be skipped (in which case msRowReIndex_v[i] is -1
3465 +​
3438 3466 vector <int> msRowReIndex_v(vmsData_p->v_antennaId1.size());
3439 - for (unsigned int i = 0; i < msRowReIndex_v.size(); i++)
3440 - msRowReIndex_v[i] = i;
3467 +​
3468 + // cout << "skipFirstTime ; " << skipFirstTime << endl;
3441 3469 ​
3442 3470 CorrelationModeMod::CorrelationMode cm = r_p->getTable().getContainer().getConfigDescription().getRowByKey(r_p->getConfigDescriptionId())->getCorrelationMode();
3443 3471 if (cm == CorrelationModeMod::CROSS_AND_AUTO and ac_xc_per_timestamp) {
3472 + // this is the case where the rows may be reordered
3473 +​
3444 3474 unsigned int numDD = r_p->getTable().getContainer().getConfigDescription().getRowByKey(r_p->getConfigDescriptionId())->getDataDescriptionId().size();
3445 3475 unsigned int numOfMSRowsPerIntegration = 0;
3446 3476 unsigned int numAntenna = r_p->getNumAntenna();
3447 3477 unsigned int numBl = r_p->getNumAntenna() * (r_p->getNumAntenna() - 1) / 2;
3448 3478 if (cm != CorrelationModeMod::CROSS_ONLY) numOfMSRowsPerIntegration += numAntenna;
3449 3479 if (cm != CorrelationModeMod::AUTO_ONLY) numOfMSRowsPerIntegration += numBl;
3450 3480 unsigned int numIntegrations = vmsData_p->v_antennaId1.size() / numOfMSRowsPerIntegration / numDD;
3451 3481 if (vmsData_p->v_antennaId1.size() % numOfMSRowsPerIntegration != 0) {
3452 3482 errstream.str("");
3453 3483 errstream << "The total number of rows to be writtren into the MS (" << vmsData_p->v_antennaId1.size()
3454 3484 <<") is not a multiple of the number of MS rows per integration (" << numOfMSRowsPerIntegration << ").";
3455 3485 errstream << "This is not normal. Aborting !";
3456 3486 throw ASDM2MSException(errstream.str());
3457 3487 }
3458 3488 ​
3459 3489 unsigned int i = 0;
3460 - for (unsigned int iDD = 0; i < numDD; i++) {
3490 + for (unsigned int iDD = 0; iDD < numDD; iDD++) {
3461 3491 unsigned int ddOffset = iDD * numIntegrations * (numAntenna + numBl);
3462 3492 for (unsigned int integration = 0; integration < numIntegrations; integration++) {
3463 3493 // First the auto correlations.
3464 - for (unsigned int iAnt = 0; iAnt < numAntenna; iAnt++)
3465 - msRowReIndex_v[i++] = ddOffset + numAntenna * integration + iAnt;
3466 -
3494 + for (unsigned int iAnt = 0; iAnt < numAntenna; iAnt++) {
3495 + if (skipFirstTime && integration == 0) {
3496 + msRowReIndex_v[i++] = -1;
3497 + } else {
3498 + msRowReIndex_v[i++] = ddOffset + numAntenna * integration + iAnt;
3499 + }
3500 + }
3467 3501 // Then the cross correlations.
3468 3502 for (unsigned iBl = 0; iBl < numBl; iBl++)
3469 - msRowReIndex_v[i++] = ddOffset + numIntegrations * numAntenna + integration * numBl + iBl;
3503 + if (skipFirstTime && integration == 0) {
3504 + msRowReIndex_v[i++] = -1;
3505 + } else {
3506 + msRowReIndex_v[i++] = ddOffset + numIntegrations * numAntenna + integration * numBl + iBl;
3507 + }
3508 + }
3509 + }
3510 + } else {
3511 + // set them in order, skipping times as appropriate
3512 + // it's not obvious that the vmsData_p rows are time sorted, but the first time there should be the time to be skipped
3513 + double timeToSkip = vmsData_p->v_time[0];
3514 + for (unsigned int i = 0; i < msRowReIndex_v.size(); i++) {
3515 + if (skipFirstTime && (vmsData_p->v_time[i] == timeToSkip)) {
3516 + msRowReIndex_v[i] = -1;
3517 + } else {
3518 + msRowReIndex_v[i] = i;
3470 3519 }
3471 3520 }
3472 3521 }
3473 3522
3474 3523 vector<vector<unsigned int> > filteredShape_vv = vmsData_p->vv_dataShape;
3475 3524 for (unsigned int ipart = 0; ipart < filteredShape_vv.size(); ipart++) {
3476 3525 if (filteredShape_vv.at(ipart).at(0) == 3) filteredShape_vv.at(ipart).at(0) = 4;
3477 3526 }
3478 -
3527 +​
3479 3528 vector<int> filteredDD;
3480 3529 for (unsigned int idd = 0; idd < vmsData_p->v_dataDescId.size(); idd++){
3481 3530 filteredDD.push_back(dataDescriptionIdx2Idx.at(vmsData_p->v_dataDescId.at(idd)));
3482 3531 }
3483 3532 ​
3484 3533 vector<float *> uncorrectedData_v;
3485 3534 vector<float *> correctedData_v;
3486 3535 ​
3487 - /* compute the UVW */
3488 - vector<double> uvw_v(3*vmsData_p->v_time.size());
3536 + // these sizes of these are not known immediately if there skipFirstTime is true
3537 + vector<double> uvw_v; // when set here, this is reordered and can be used as is when writing to the uncorrected MS
3538 + vector<double> weight_v; // these are the weights, in the order seen. They must be reordered to use.
3539 + vector<double> correctedWeight_v; // to be put into the MS, can only be set later when other corrected column values are set
3540 + vector<double> uncorrectedWeight_v; // to be put into the MS. May be set here when skipFirstTime is not true.
3541 + vector<double> sigma_v; // the sigmas, in the order seen. They must be reordered to use.
3542 + vector<double> correctedSigma_v; // to be put into the MS, can only be set later when other corrected column values are set
3543 + vector<double> uncorrectedSigma_v; // to be put into the MS. May be set here when skipFirstTime is not true.
3544 +​
3545 + /* compute the UVW - do this for all times, time skipping happens later */
3489 3546 vector<casacore::Vector<casacore::Double> > vv_uvw(vmsData_p->v_time.size());
3490 3547 #if DDPRIORITY
3491 3548 uvwCoords.uvw_bl(r_p, sdmBinData.timeSequence(), e_query_cm,
3492 3549 sdmbin::SDMBinData::dataOrder(),
3493 3550 vv_uvw);
3494 3551 #else
3495 3552 uvwCoords.uvw_bl(r_p, vmsData_p->v_timeCentroid, e_query_cm,
3496 3553 sdmbin::SDMBinData::dataOrder(),
3497 3554 vv_uvw);
3498 3555 #endif
3499 3556 ​
3500 - /*
3501 - ** Let's apply the reindexing on the UVW coordinates !!!
3502 - */
3503 - int k = 0;
3504 - for (unsigned int iUvw = 0; iUvw < vv_uvw.size(); iUvw++) {
3505 - uvw_v[k++] = vv_uvw[msRowReIndex_v[iUvw]](0);
3506 - uvw_v[k++] = vv_uvw[msRowReIndex_v[iUvw]](1);
3507 - uvw_v[k++] = vv_uvw[msRowReIndex_v[iUvw]](2);
3508 - }
3509 -​
3510 - /*
3511 - ** Let's apply the reindexing on weight and sigma.
3512 - */
3513 - vector<double> weight_v(vmsData_p->v_time.size());
3514 - vector<double> correctedWeight_v(vmsData_p->v_time.size());
3515 -​
3516 - vector<double> sigma_v(vmsData_p->v_time.size());
3517 - vector<double> correctedSigma_v(vmsData_p->v_time.size());
3557 + // set sigma_v and weight_v first, in order.
3558 + weight_v.resize(vmsData_p->v_time.size());
3559 + sigma_v.resize(vmsData_p->v_time.size());
3518 3560 for (unsigned int i = 0; i < weight_v.size(); i++) {
3519 - weight_v[i] = vmsData_p->v_exposure.at(msRowReIndex_v[i]) * effectiveBwPerDD_m[filteredDD[msRowReIndex_v[i]]];
3520 - if (vmsData_p->v_antennaId1[msRowReIndex_v[i]] != vmsData_p->v_antennaId2[msRowReIndex_v[i]])
3561 + weight_v[i] = vmsData_p->v_exposure[i] * effectiveBwPerDD_m[filteredDD[i]];
3562 + if (vmsData_p->v_antennaId1[i] != vmsData_p->v_antennaId2[i])
3521 3563 weight_v[i] *= 2.0;
3522 - correctedWeight_v[i] = weight_v[i];
3523 -
3564 +​
3524 3565 if (weight_v[i] == 0.0) weight_v[i] = 1.0;
3566 +​
3525 3567 sigma_v[i] = 1.0 / sqrt(weight_v[i]);
3526 3568 ​
3527 - correctedSigma_v[i] = sigma_v[i];
3528 3569 }
3570 +
3571 + if (!skipFirstTime) {
3572 + // several values can be fully filled out in this case.
3573 + // it's more efficient to just do that, most of the time this block will be executed
3574 +​
3575 + /*
3576 + ** Let's apply the reindexing on the UVW coordinates !!!
3577 + */
3578 + uvw_v.resize(3*vmsData_p->v_time.size());
3579 + int k = 0;
3580 + for (unsigned int iUvw = 0; iUvw < vv_uvw.size(); iUvw++) {
3581 + uvw_v[k++] = vv_uvw[msRowReIndex_v[iUvw]](0);
3582 + uvw_v[k++] = vv_uvw[msRowReIndex_v[iUvw]](1);
3583 + uvw_v[k++] = vv_uvw[msRowReIndex_v[iUvw]](2);
3584 + }
3529 3585 ​
3530 - // Here we make the assumption that the State is the same for all the antennas and let's use the first State found in the vector stateId contained in the ASDM Main Row
3531 - // int asdmStateIdx = r_p->getStateId().at(0).getTagValue();
3532 - vector<int> msStateId_v(vmsData_p->v_m_data.size(), stateIdx2Idx[r_p]);
3586 + /*
3587 + ** Let's apply the reindexing on weight and sigma.
3588 + */
3589 + uncorrectedWeight_v.resize(weight_v.size());
3590 + uncorrectedSigma_v.resize(weight_v.size());
3591 + for (unsigned int i = 0; i < weight_v.size(); i++) {
3592 + uncorrectedWeight_v[i] = weight_v.at(msRowReIndex_v[i]);
3593 + uncorrectedSigma_v[i] = sigma_v.at(msRowReIndex_v[i]);
3594 + }
3595 + }
3533 3596 ​
3534 3597 ComplexDataFilter cdf;
3535 3598 map<AtmPhaseCorrectionMod::AtmPhaseCorrection, float*>::const_iterator iter;
3536 3599 ​
3600 + vector<double> uncorrectedTime_v;
3601 + vector<int> uncorrectedAntennaId1_v;
3602 + vector<int> uncorrectedAntennaId2_v;
3603 + vector<int> uncorrectedFeedId1_v;
3604 + vector<int> uncorrectedFeedId2_v;
3605 + vector<int> uncorrectedFieldId_v;
3606 + vector<int> uncorrectedFilteredDD_v;
3607 + vector<vector< unsigned int> > uncorrectedFilteredShape_vv;
3608 + vector<double> uncorrectedInterval_v;
3609 + vector<double> uncorrectedExposure_v;
3610 + vector<double> uncorrectedTimeCentroid_v;
3611 + vector<unsigned int> uncorrectedFlag_v;
3612 +​
3537 3613 vector<double> correctedTime_v;
3538 3614 vector<int> correctedAntennaId1_v;
3539 3615 vector<int> correctedAntennaId2_v;
3540 3616 vector<int> correctedFeedId1_v;
3541 3617 vector<int> correctedFeedId2_v;
3542 3618 vector<int> correctedFieldId_v;
3543 3619 vector<int> correctedFilteredDD_v;
3544 3620 vector<vector< unsigned int> > correctedFilteredShape_vv;
3545 3621 vector<double> correctedInterval_v;
3546 3622 vector<double> correctedExposure_v;
3547 3623 vector<double> correctedTimeCentroid_v;
3548 - vector<int> correctedMsStateId_v(msStateId_v);
3549 3624 vector<double> correctedUvw_v;
3550 3625 vector<unsigned int> correctedFlag_v;
3551 3626 ​
3552 3627 Tag configDescriptionId = r_p -> getConfigDescriptionId();
3553 3628 ConfigDescriptionTable & cfgDescT = r_p -> getTable() . getContainer() . getConfigDescription();
3554 3629 ConfigDescriptionRow * cfgDescR_p = cfgDescT.getRowByKey(configDescriptionId);
3555 3630 const vector<AtmPhaseCorrectionMod::AtmPhaseCorrection >& apc_v = cfgDescR_p->getAtmPhaseCorrection();
3556 3631 bool subscanHasCorrectedData = std::find(apc_v.begin(), apc_v.end(), AtmPhaseCorrectionMod::AP_CORRECTED)!=apc_v.end();
3557 3632 ​
3558 3633 // Do we have to fill an MS with uncorrected data + radiometric data (radiometric data are considered as uncorrected data) ?
3559 3634 // Apply here the redindexing on uncorrected data !!!!
3560 3635 //
3561 3636 for (unsigned int iData = 0; iData < vmsData_p->v_m_data.size(); iData++) {
3637 + if (skipFirstTime && msRowReIndex_v[iData] < 0) {
3638 + continue;
3639 + }
3640 +​
3562 3641 if ((iter=vmsData_p->v_m_data.at(msRowReIndex_v[iData]).find(AtmPhaseCorrectionMod::AP_UNCORRECTED)) != vmsData_p->v_m_data.at(msRowReIndex_v[iData]).end()){
3563 3642 uncorrectedData_v.push_back(cdf.to4Pol(vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(0),
3564 3643 vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(1),
3565 3644 iter->second));
3645 + if (skipFirstTime) {
3646 + // uncorrected values must also be explicitly added here since they can not be put in as entire vectors
3647 + uncorrectedTime_v.push_back(vmsData_p->v_time.at(msRowReIndex_v[iData]));
3648 + uncorrectedAntennaId1_v.push_back(vmsData_p->v_antennaId1.at(msRowReIndex_v[iData]));
3649 + uncorrectedAntennaId2_v.push_back(vmsData_p->v_antennaId2.at(msRowReIndex_v[iData]));
3650 + uncorrectedFeedId1_v.push_back(vmsData_p->v_feedId1.at(msRowReIndex_v[iData]));
3651 + uncorrectedFeedId2_v.push_back(vmsData_p->v_feedId2.at(msRowReIndex_v[iData]));
3652 + uncorrectedFilteredDD_v.push_back(filteredDD.at(msRowReIndex_v[iData]));
3653 + uncorrectedFilteredShape_vv.push_back(filteredShape_vv.at(msRowReIndex_v[iData]));
3654 + uncorrectedFieldId_v.push_back(vmsData_p->v_fieldId.at(msRowReIndex_v[iData]));
3655 + uncorrectedInterval_v.push_back(vmsData_p->v_interval.at(msRowReIndex_v[iData]));
3656 + uncorrectedExposure_v.push_back(vmsData_p->v_exposure.at(msRowReIndex_v[iData]));
3657 + uncorrectedTimeCentroid_v.push_back(vmsData_p->v_timeCentroid.at(msRowReIndex_v[iData]));
3658 + uvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](0));
3659 + uvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](1));
3660 + uvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](2));
3661 + uncorrectedFlag_v.push_back(vmsData_p->v_flag.at(msRowReIndex_v[iData]));
3662 + uncorrectedWeight_v.push_back(weight_v.at(msRowReIndex_v[iData]));
3663 + uncorrectedSigma_v.push_back(sigma_v.at(msRowReIndex_v[iData]));
3664 + }
3566 3665 }
3567 -
3666 +​
3568 3667 // Have we asked to write an MS with corrected data + radiometric data ?
3569 3668
3570 3669 // Are we with radiometric data ? Then we assume that the data are labelled AP_UNCORRECTED.
3571 3670 if (sdmBinData.processorType(r_p) == RADIOMETER) {
3572 3671 if ((iter=vmsData_p->v_m_data.at(msRowReIndex_v[iData]).find(AtmPhaseCorrectionMod::AP_UNCORRECTED)) != vmsData_p->v_m_data.at(msRowReIndex_v[iData]).end()){
3672 +
3573 3673 correctedTime_v.push_back(vmsData_p->v_time.at(msRowReIndex_v[iData]));
3574 3674 correctedAntennaId1_v.push_back(vmsData_p->v_antennaId1.at(msRowReIndex_v[iData]));
3575 3675 correctedAntennaId2_v.push_back(vmsData_p->v_antennaId2.at(msRowReIndex_v[iData]));
3576 3676 correctedFeedId1_v.push_back(vmsData_p->v_feedId1.at(msRowReIndex_v[iData]));
3577 3677 correctedFeedId2_v.push_back(vmsData_p->v_feedId2.at(msRowReIndex_v[iData]));
3578 3678 correctedFilteredDD_v.push_back(filteredDD.at(msRowReIndex_v[iData]));
3579 3679 correctedFilteredShape_vv.push_back(filteredShape_vv.at(msRowReIndex_v[iData]));
3580 3680 correctedFieldId_v.push_back(vmsData_p->v_fieldId.at(msRowReIndex_v[iData]));
3581 3681 correctedInterval_v.push_back(vmsData_p->v_interval.at(msRowReIndex_v[iData]));
3582 3682 correctedExposure_v.push_back(vmsData_p->v_exposure.at(msRowReIndex_v[iData]));
3583 3683 correctedTimeCentroid_v.push_back(vmsData_p->v_timeCentroid.at(msRowReIndex_v[iData]));
3584 3684 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](0));
3585 3685 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](1));
3586 3686 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](2));
3587 - correctedData_v.push_back(uncorrectedData_v.at(iData)); // <------------Attention here the uncorrected data have been already re ordered !
3687 + // this is probably the most recent uncorrectedData, but it might not be - else why the if blocks here
3688 + correctedData_v.push_back(cdf.to4Pol(vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(0),
3689 + vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(1),
3690 + iter->second));
3588 3691 correctedFlag_v.push_back(vmsData_p->v_flag.at(msRowReIndex_v[iData]));
3692 + correctedWeight_v.push_back(weight_v.at(msRowReIndex_v[iData]));
3693 + correctedSigma_v.push_back(sigma_v.at(msRowReIndex_v[iData]));
3589 3694 }
3590 3695 }
3591 3696 else { // We assume that we are in front of CORRELATOR data, but do we have corrected data on that specific subscan ?
3592 3697 if (subscanHasCorrectedData) {
3593 3698 // Then we know that we have AP_CORRECTED data.
3594 3699 if (vmsData_p->v_antennaId1.at(msRowReIndex_v[iData]) == vmsData_p->v_antennaId2.at(msRowReIndex_v[iData]) ) {
3595 3700 /*
3596 3701 ** do not forget to prepend the autodata copied from the uncorrected data, because the lower layers of the software do not put the (uncorrected) autodata in the
3597 3702 ** corrected data.
3598 3703 */
3601 3706 correctedAntennaId2_v.push_back(vmsData_p->v_antennaId2.at(msRowReIndex_v[iData]));
3602 3707 correctedFeedId1_v.push_back(vmsData_p->v_feedId1.at(msRowReIndex_v[iData]));
3603 3708 correctedFeedId2_v.push_back(vmsData_p->v_feedId2.at(msRowReIndex_v[iData]));
3604 3709 correctedFilteredDD_v.push_back(filteredDD.at(msRowReIndex_v[iData]));
3605 3710 correctedFieldId_v.push_back(vmsData_p->v_fieldId.at(msRowReIndex_v[iData]));
3606 3711 correctedInterval_v.push_back(vmsData_p->v_interval.at(msRowReIndex_v[iData]));
3607 3712 correctedExposure_v.push_back(vmsData_p->v_exposure.at(msRowReIndex_v[iData]));
3608 3713 correctedTimeCentroid_v.push_back(vmsData_p->v_timeCentroid.at(msRowReIndex_v[iData]));
3609 3714 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](0));
3610 3715 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](1));
3611 - correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](2));
3612 - correctedData_v.push_back(uncorrectedData_v.at(iData)); // <-------- Here we re-use the autodata already present in the uncorrected data and which have been
3613 - // already reindexed !!!
3716 + correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](2));
3717 + // this is probably the most recent uncorrectedData, but it might not be - else why the if blocks here
3718 + correctedData_v.push_back(cdf.to4Pol(vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(0),
3719 + vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(1),
3720 + iter->second));
3614 3721 correctedFlag_v.push_back(vmsData_p->v_flag.at(msRowReIndex_v[iData]));
3615 3722 correctedFilteredShape_vv.push_back(filteredShape_vv.at(msRowReIndex_v[iData]));
3723 + correctedWeight_v.push_back(weight_v.at(msRowReIndex_v[iData]));
3724 + correctedSigma_v.push_back(sigma_v.at(msRowReIndex_v[iData]));
3725 +​
3616 3726
3617 3727 }
3618 3728 else {
3619 3729 /*
3620 3730 ** And now finally the correlation corrected data.
3621 3731 */
3622 3732 correctedTime_v.push_back(vmsData_p->v_time.at(msRowReIndex_v[iData]));
3623 3733 correctedAntennaId1_v.push_back(vmsData_p->v_antennaId1.at(msRowReIndex_v[iData]));
3624 3734 correctedAntennaId2_v.push_back(vmsData_p->v_antennaId2.at(msRowReIndex_v[iData]));
3625 3735 correctedFeedId1_v.push_back(vmsData_p->v_feedId1.at(msRowReIndex_v[iData]));
3632 3742 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](0));
3633 3743 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](1));
3634 3744 correctedUvw_v.push_back(vv_uvw[msRowReIndex_v[iData]](2));
3635 3745 iter=vmsData_p->v_m_data.at(msRowReIndex_v[iData]).find(AtmPhaseCorrectionMod::AP_CORRECTED);
3636 3746 float* theData = cdf.to4Pol(vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(0),
3637 3747 vmsData_p->vv_dataShape.at(msRowReIndex_v[iData]).at(1),
3638 3748 iter->second);
3639 3749 correctedData_v.push_back(theData);
3640 3750 correctedFlag_v.push_back(vmsData_p->v_flag.at(msRowReIndex_v[iData]));
3641 3751 correctedFilteredShape_vv.push_back(filteredShape_vv.at(msRowReIndex_v[iData]));
3752 + correctedWeight_v.push_back(weight_v.at(msRowReIndex_v[iData]));
3753 + correctedSigma_v.push_back(sigma_v.at(msRowReIndex_v[iData]));
3642 3754 }
3643 3755 }
3644 3756 }
3645 3757 }
3646 3758
3647 3759 if (uncorrectedData_v.size() > 0 && (msFillers.find(AP_UNCORRECTED) != msFillers.end())) {
3648 3760 if (! mute) { // Here we make the assumption that we have always uncorrected data. This realistic even if not totally rigorous.
3649 3761 ​
3650 - vector<double> uncorrectedTime_v = reorder<double>( vmsData_p->v_time, msRowReIndex_v);
3651 - vector<int> uncorrectedAntennaId1_v = reorder<int> (vmsData_p->v_antennaId1, msRowReIndex_v ) ;
3652 - vector<int> uncorrectedAntennaId2_v = reorder<int>(vmsData_p->v_antennaId2, msRowReIndex_v );
3653 - vector<int> uncorrectedFeedId1_v = reorder<int>(vmsData_p->v_feedId1, msRowReIndex_v );
3654 - vector<int> uncorrectedFeedId2_v = reorder<int>(vmsData_p->v_feedId2, msRowReIndex_v );
3655 - vector<int> uncorrectedFieldId_v = reorder<int>(vmsData_p->v_fieldId, msRowReIndex_v );
3656 - vector<int> uncorrectedFilteredDD_v = reorder<int>(filteredDD, msRowReIndex_v );
3657 - vector<vector< unsigned int> > uncorrectedFilteredShape_vv = reorder<vector< unsigned int> >(filteredShape_vv, msRowReIndex_v );
3658 - vector<double> uncorrectedInterval_v = reorder<double>(vmsData_p->v_interval, msRowReIndex_v );
3659 - vector<double> uncorrectedExposure_v = reorder<double>(vmsData_p->v_exposure, msRowReIndex_v );
3660 - vector<double> uncorrectedTimeCentroid_v = reorder<double>(vmsData_p->v_timeCentroid, msRowReIndex_v );
3661 - vector<unsigned int> uncorrectedFlag_v = reorder<unsigned int>(vmsData_p->v_flag, msRowReIndex_v );
3762 + if (!skipFirstTime) {
3763 + // need to set these now, but they can be set as vectors that need to be reordered
3764 + uncorrectedTime_v = reorder<double>( vmsData_p->v_time, msRowReIndex_v);
3765 + uncorrectedAntennaId1_v = reorder<int> (vmsData_p->v_antennaId1, msRowReIndex_v ) ;
3766 + uncorrectedAntennaId2_v = reorder<int>(vmsData_p->v_antennaId2, msRowReIndex_v );
3767 + uncorrectedFeedId1_v = reorder<int>(vmsData_p->v_feedId1, msRowReIndex_v );
3768 + uncorrectedFeedId2_v = reorder<int>(vmsData_p->v_feedId2, msRowReIndex_v );
3769 + uncorrectedFieldId_v = reorder<int>(vmsData_p->v_fieldId, msRowReIndex_v );
3770 + uncorrectedFilteredDD_v = reorder<int>(filteredDD, msRowReIndex_v );
3771 + uncorrectedFilteredShape_vv = reorder<vector< unsigned int> >(filteredShape_vv, msRowReIndex_v );
3772 + uncorrectedInterval_v = reorder<double>(vmsData_p->v_interval, msRowReIndex_v );
3773 + uncorrectedExposure_v = reorder<double>(vmsData_p->v_exposure, msRowReIndex_v );
3774 + uncorrectedTimeCentroid_v = reorder<double>(vmsData_p->v_timeCentroid, msRowReIndex_v );
3775 + uncorrectedFlag_v = reorder<unsigned int>(vmsData_p->v_flag, msRowReIndex_v );
3776 + // uvw_v, uncorrectedSigma_v, and uncorrectedWeight_v have all been previously set
3777 + }
3778 +​
3779 + // Here we make the assumption that the State is the same for all the antennas and let's use the first State found in the vector stateId contained in the ASDM Main Row
3780 + // state must have already been filled so that the stateIdx2IDx map is available to extract the state ID for this main row pointer (r_p).
3781 + vector<int> msStateId_v(uncorrectedTime_v.size(), stateIdx2Idx[r_p]);
3662 3782 ​
3663 3783 msFillers[AP_UNCORRECTED]->addData(complexData,
3664 3784 uncorrectedTime_v , // this is already time midpoint
3665 3785 uncorrectedAntennaId1_v,
3666 3786 uncorrectedAntennaId2_v,
3667 3787 uncorrectedFeedId1_v,
3668 3788 uncorrectedFeedId2_v,
3669 3789 uncorrectedFilteredDD_v,
3670 3790 (int) vmsData_p->processorId,
3671 - (vector<int> &) vmsData_p->v_fieldId,
3672 - (vector<double> &) vmsData_p->v_interval,
3673 - (vector<double> &) vmsData_p->v_exposure,
3674 - (vector<double> &)vmsData_p->v_timeCentroid,
3791 + uncorrectedFieldId_v,
3792 + uncorrectedInterval_v,
3793 + uncorrectedExposure_v,
3794 + uncorrectedTimeCentroid_v,
3675 3795 (int) r_p->getScanNumber(),
3676 3796 0, // Array Id
3677 3797 (int) r_p->getExecBlockId().getTagValue(), // Observation Id
3678 3798 msStateId_v,
3679 - uvw_v,
3799 + uvw_v,
3680 3800 uncorrectedFilteredShape_vv, // vmsData_p->vv_dataShape after filtering the case numCorr == 3
3681 3801 uncorrectedData_v,
3682 3802 uncorrectedFlag_v,
3683 - weight_v,
3684 - sigma_v);
3803 + uncorrectedWeight_v,
3804 + uncorrectedSigma_v);
3685 3805 }
3686 3806 }
3687 3807 ​
3808 +​
3688 3809 if (correctedData_v.size() > 0 && (msFillers.find(AP_CORRECTED) != msFillers.end())) {
3689 3810 if (! mute) {
3811 + // Here we make the assumption that the State is the same for all the antennas and let's use the first State found in the vector stateId contained in the ASDM Main Row
3812 + // state must have already been filled so that the stateIdx2IDx map is available to extract the state ID for this main row pointer (r_p).
3813 + vector<int> correctedMsStateId_v(correctedTime_v.size(), stateIdx2Idx[r_p]);
3690 3814 msFillers[AP_CORRECTED]->addData(complexData,
3691 3815 correctedTime_v, // this is already time midpoint
3692 3816 correctedAntennaId1_v,
3693 3817 correctedAntennaId2_v,
3694 3818 correctedFeedId1_v,
3695 3819 correctedFeedId2_v,
3696 3820 correctedFilteredDD_v,
3697 3821 vmsData_p->processorId,
3698 3822 correctedFieldId_v,
3699 3823 correctedInterval_v,
3707 3831 correctedFilteredShape_vv, // vmsData_p->vv_dataShape after filtering the case numCorr == 3
3708 3832 correctedData_v,
3709 3833 correctedFlag_v,
3710 3834 correctedWeight_v,
3711 3835 correctedSigma_v);
3712 3836 }
3713 3837 }
3714 3838 if (debug) cout << "fillMain : exiting" << endl;
3715 3839 }
3716 3840 ​
3717 -​
3718 -void testFunc(string& tstr) {
3719 - cerr<<tstr<<endl;
3720 -}
3721 -​
3722 -/**
3723 - * compute the UVW (put in a method to keep sepearate from fillMain for
3724 - parallel case) returns a Matrix, mat_uvw for data passing thread-safe
3725 - way
3726 -*/
3727 -void calcUVW(MainRow* r_p,
3728 - SDMBinData& sdmBinData,
3729 - const VMSData* vmsData_p,
3730 - UvwCoords& uvwCoords,
3731 - casacore::Matrix< casacore::Double >& mat_uvw) {
3732 -​
3733 - vector< casacore::Vector<casacore::Double> > vv_uvw;
3734 - mat_uvw.resize(3,vmsData_p->v_time.size());
3735 -​
3736 -#if DDPRIORITY
3737 - uvwCoords.uvw_bl(r_p, sdmBinData.timeSequence(), e_query_cm,
3738 - sdmbin::SDMBinData::dataOrder(),
3739 - vv_uvw);
3740 -#else
3741 - uvwCoords.uvw_bl(r_p, vmsData_p->v_timeCentroid, e_query_cm,
3742 - sdmbin::SDMBinData::dataOrder(),
3743 - vv_uvw);
3744 -#endif
3745 -​
3746 - //put in a Matrix
3747 - for (unsigned int iUvw = 0; iUvw < vv_uvw.size(); iUvw++) {
3748 - mat_uvw(iUvw, 0) = vv_uvw[iUvw](0);
3749 - mat_uvw(iUvw, 1) = vv_uvw[iUvw](1);
3750 - mat_uvw(iUvw, 2) = vv_uvw[iUvw](2);
3751 - }
3752 -}
3753 -​
3754 3841 /**
3842 + * This function has not been used in production. It therefore has not kept up with
3843 + * ongoing development elsewhere here. It should be seen as an example for possible
3844 + * future multithreading work and should not be use as is.
3845 + *
3846 + **
3755 3847 * This function fills the MS Main table from an ASDM Main table which refers to correlator data.
3756 3848 * designed for multithreading........
3757 3849 *
3758 3850 * given:
3759 3851 * @parameter rowNum an integer expected to contain the number of the row being processed.
3760 3852 * @parameter r_p a pointer to the MainRow being processed.
3761 3853 * @parameter sdmBinData a reference to the SDMBinData containing a lot of information about the binary data being processed. Useful to know the requested ordering of data.
3762 3854 * @parameter uvwCoords a reference to the UVW calculator.
3763 3855 * @parameter complexData a bool which says if the DATA is going to be filled (true) or if it will be the FLOAT_DATA (false).
3764 3856 * @parameter mute if the value of this parameter is false then nothing is written in the MS .
3982 4074 correctedMsStateId,
3983 4075 correctedUvw,
3984 4076 filteredShape, // vmsData_p->vv_dataShape after filtering the case numCorr == 3
3985 4077 correctedData,
3986 4078 correctedFlag);
3987 4079 }
3988 4080 }
3989 4081 if (debug) cout << "fillMain_mt : exiting" << endl;
3990 4082 }
3991 4083 #endif
3992 -​
3993 -void fillSysPower_aux (const vector<SysPowerRow *>& sysPowers, map<AtmPhaseCorrection, ASDM2MSFiller*>& msFillers_m) {
3994 - LOGENTER("fillSysPower_aux");
4084 +
4085 + void fillSysPower_aux (const vector<SysPowerRow *>& sysPowers, map<AtmPhaseCorrection, ASDM2MSFiller*>& msFillers_m) {
4086 + LOGENTER("fillSysPower_aux");
3995 4087 ​
3996 4088 vector<int> antennaId;
3997 4089 vector<int> spectralWindowId;
3998 4090 vector<int> feedId;
3999 4091 vector<double> time;
4000 4092 vector<double> interval;
4001 4093 vector<int> numReceptor;
4002 4094 vector<float> switchedPowerDifference;
4003 4095 vector<float> switchedPowerSum;
4004 4096 vector<float> requantizerGain;
4295 4387 ​
4296 4388 bool mute = false;
4297 4389 ​
4298 4390 bool doparallel = false;
4299 4391 ​
4300 4392 bool ac_xc_per_timestamp = false; // for the time being the option is 'preserve the old order'
4301 4393 ​
4302 4394 bool interpolate_ephemeris = false ;
4303 4395 bool tabulate_ephemeris_polynomials = false;
4304 4396 double polyephem_tabtimestep = 0.0;
4397 +​
4398 + bool checkDupInts = true;
4305 4399
4306 4400 // Process command line options and parameters.
4307 4401 po::variables_map vm;
4308 4402 ​
4309 4403 try {
4310 4404 ​
4311 4405 ​
4312 4406 // Declare the supported options.
4313 4407 ​
4314 4408 po::options_description generic("Converts an ASDM dataset into a CASA measurement set.\n"
4351 4445 // Hidden options, will be allowed both on command line and
4352 4446 // in config file, but will not be shown to the user.
4353 4447 po::options_description hidden("Hidden options");
4354 4448 hidden.add_options()
4355 4449 ("asdm-directory", po::value< string >(), "asdm directory")
4356 4450 ;
4357 4451 hidden.add_options()
4358 4452 ("ms-directory-prefix", po::value< string >(), "ms directory prefix")
4359 4453 ;
4360 4454 ​
4455 + // This is only used by the test programs. Not indended for general use.
4456 + hidden.add_options()
4457 + ("checkdupints", po::value<bool>()->default_value("true"),"a value of false turns off checks for duplicate integration times in RADIOMETER data")
4458 + ;
4459 +​
4361 4460 po::options_description cmdline_options;
4362 4461 cmdline_options.add(generic).add(hidden);
4363 4462
4364 4463 po::positional_options_description p;
4365 4464 p.add("asdm-directory", 1);
4366 4465 p.add("ms-directory-prefix", 1);
4367 4466
4368 4467 // Parse the command line and retrieve the options and parameters.
4369 4468 po::store(po::command_line_parser(argc, argv).options(cmdline_options).positional(p).run(), vm);
4370 4469 ​
4590 4689 ​
4591 4690 if (lazy && !lazyModeOK) {
4592 4691 lazy = false;
4593 4692 infostream.str("");
4594 4693 infostream << "The lazy filler can not be used with any input data selection (icm, isrt, its)." << endl;
4595 4694 infostream << "The non-lazy version of the filler will be used." << endl;
4596 4695 // an option is being ignored, warn the user even if not verbose
4597 4696 warning(infostream.str());
4598 4697 }
4599 4698 ​
4699 + // check for duplicate integrations in RADIOMETER data?
4700 + checkDupInts = vm["checkdupints"].as< bool >();
4701 + if (debug) {
4702 + cout << "checkDupInts : " << checkDupInts << endl;
4703 + }
4704 +​
4600 4705 // Do we consider another order than ac_xc_per_timestamp ?
4601 4706 ac_xc_per_timestamp = (vm.count("ac-xc-per-timestamp") == 0) ? false :
4602 4707 boost::algorithm::to_lower_copy(vm["ac-xc-per-timestamp"].as<string>()) == "yes";
4603 4708 ​
4604 4709 // Do we want to tabulate polynomial present in the ephemeris table ?
4605 4710 tabulate_ephemeris_polynomials = (vm.count("polyephem-tabtimestep") != 0);
4606 4711
4607 4712 if (tabulate_ephemeris_polynomials) {
4608 4713 // If we tabluate then ignore all the other options about ephemeris
4609 4714 // infostream.str();
5663 5768 info(infostream.str());
5664 5769 int nPointing = v.size();
5665 5770 ​
5666 5771 if (nPointing > 0) {
5667 5772 ​
5668 5773 // Check some assertions.
5669 5774 //
5670 5775 // All rows of ASDM-Pointing must have their attribute usePolynomials equal to false
5671 5776 // and their numTerm attribute equal to 1. Use the opportunity of this check
5672 5777 // to compute the number of rows to be created in the MS-Pointing by summing
5673 - // all the numSample attributes values.
5778 + // all the numSample attributes values. Watch for duplicate times due and avoid.
5674 5779 //
5675 5780 int numMSPointingRows = 0;
5781 +​
5782 + // initialize the lastTime to 0.0 for all antennaIds
5783 + map<Tag, double> lastTime;
5784 + const AntennaTable& antennaT = ds->getAntenna();
5785 + const vector<AntennaRow *>& vAntRow = antennaT.get();
5786 + for (unsigned int i=0; i < vAntRow.size(); i++) {
5787 + lastTime[vAntRow[i]->getAntennaId()] = 0.0;
5788 + }
5789 +​
5790 + // set this to true for any rows where the first element should be skipped because the time is a duplicate
5791 + vector<bool> vSkipFirst(v.size(),false);
5792 +​
5676 5793 for (unsigned int i = 0; i < v.size(); i++) {
5677 5794 if (v[i]->getUsePolynomials()) {
5678 5795 errstream.str("");
5679 5796 errstream << "Found usePolynomials equal to true at row #" << i <<". Can't go further.";
5680 5797 error(errstream.str());
5681 5798 }
5682 5799 ​
5683 - numMSPointingRows += v[i]->getNumSample();
5800 + // look for duplicate rows - time at the start of row is near the time at the end of the previous row for that antennaId
5801 + int numSample = v[i]->getNumSample();
5802 + Tag antId = v[i]->getAntennaId();
5803 + double tfirst = 0.0;
5804 + double tlast = 0.0;
5805 + if (v[i]->isSampledTimeIntervalExists()) {
5806 + tfirst = ((double) (v[i]->getSampledTimeInterval().at(0).getStart().get())) / ArrayTime::unitsInASecond;
5807 + tlast = ((double) (v[i]->getSampledTimeInterval().at(numSample-1).getStart().get())) / ArrayTime::unitsInASecond;
5808 + } else {
5809 + double tstart = ((double) v[i]->getTimeInterval().getStart().get()) / ArrayTime::unitsInASecond;
5810 + double tint = ((double) v[i]->getTimeInterval().getDuration().get()) / ArrayTime::unitsInASecond / numSample;
5811 + tfirst = tstart + tint/2.0;
5812 + tlast = tfirst + tint*(numSample-1);
5813 + }
5814 + if (casacore::near(tfirst,lastTime[antId])) {
5815 + infostream.str("");
5816 + infostream << "First time in Pointing row " << i << " for antenna " << antId << " is near the previous last time for that antenna, skipping that first time in the output MS POINTING table" << endl;
5817 + info(infostream.str());
5818 + numSample--;
5819 + lastTime[antId] = tlast;
5820 + vSkipFirst[i] = true;
5821 + }
5822 + lastTime[antId] = tlast;
5823 + numMSPointingRows += numSample;
5684 5824 }
5685 5825 ​
5686 5826 //
5687 5827 // Ok now we have verified the assertions and we know the number of rows
5688 5828 // to be created into the MS-Pointing, we can proceed.
5689 5829 ​
5690 5830 PointingRow* r = 0;
5691 5831 ​
5692 5832 vector<int> antenna_id_(numMSPointingRows, 0);
5693 5833 vector<double> time_(numMSPointingRows, 0.0);
5762 5902 for (unsigned int ii = 0; ii < 3; ii++) {
5763 5903 matrix3x3.push_back(cartesian1); // cartesian1 is used here just as a way to get a 1D vector of size 3.
5764 5904 }
5765 5905 double PSI = M_PI_2;
5766 5906 double THETA;
5767 5907 double PHI;
5768 5908
5769 5909 vector<ArrayTimeInterval> timeInterval ;
5770 5910 if (r->isSampledTimeIntervalExists()) timeInterval = r->getSampledTimeInterval();
5771 5911 ​
5912 + // and now insert the values into the MS POINTING table
5913 + // watch for the skipped initial sample
5914 + int numSampleUsed = numSample;
5915 + if (vSkipFirst.at(i)) numSampleUsed--;
5916 +​
5772 5917 // Use 'fill' from algorithm for the cases where values remain constant.
5773 5918 // ANTENNA_ID
5774 - fill(antenna_id_.begin()+iMSPointingRow, antenna_id_.begin()+iMSPointingRow+numSample, antennaId);
5919 + fill(antenna_id_.begin()+iMSPointingRow, antenna_id_.begin()+iMSPointingRow+numSampleUsed, antennaId);
5775 5920 ​
5776 5921 // TRACKING
5777 - fill(tracking_.begin()+iMSPointingRow, tracking_.begin()+iMSPointingRow+numSample, pointingTracking);
5922 + fill(tracking_.begin()+iMSPointingRow, tracking_.begin()+iMSPointingRow+numSampleUsed, pointingTracking);
5778 5923 ​
5779 5924 // OVER_THE_TOP
5780 5925 if (overTheTopExists4All)
5781 5926 // it's present everywhere
5782 - fill(v_overTheTop_.begin()+iMSPointingRow, v_overTheTop_.begin()+iMSPointingRow+numSample,
5927 + fill(v_overTheTop_.begin()+iMSPointingRow, v_overTheTop_.begin()+iMSPointingRow+numSampleUsed,
5783 5928 r->getOverTheTop());
5784 5929 else if (r->isOverTheTopExists()) {
5785 5930 // it's present only in some rows.
5786 5931 s_overTheTop saux ;
5787 - saux.start = iMSPointingRow; saux.len = numSample; saux.value = r->getOverTheTop();
5932 + saux.start = iMSPointingRow; saux.len = numSampleUsed; saux.value = r->getOverTheTop();
5788 5933 v_s_overTheTop_.push_back(saux);
5789 5934 }
5790 5935
5791 5936 // Use an explicit loop for the other values.
5792 5937 for (int j = 0 ; j < numSample; j++) { // ... must be expanded in numSample MS-Pointing rows.
5938 + if (j == 0 && vSkipFirst.at(i)) continue; // the first element is to be skipped
5793 5939 ​
5794 5940 // TIME and INTERVAL
5795 5941 if (r->isSampledTimeIntervalExists()) { //if sampledTimeInterval is present use its values.
5796 5942 // Here the size of timeInterval will have to be checked against numSample !!
5797 5943 interval_[iMSPointingRow] = ((double) timeInterval.at(j).getDuration().get()) / ArrayTime::unitsInASecond ;
5798 5944 time_[iMSPointingRow] = ((double) timeInterval.at(j).getStart().get()) / ArrayTime::unitsInASecond
5799 5945 + interval_[iMSPointingRow]/2;
5800 5946 }
5801 5947 else { // otherwise compute TIMEs and INTERVALs from the first values.
5802 5948 interval_[iMSPointingRow] = interval;
6555 6701 }
6556 6702 catch ( std::exception & e) {
6557 6703 errstream.str("");
6558 6704 errstream << e.what();
6559 6705 error(errstream.str());
6560 6706 }
6561 6707 ​
6562 6708 // And then finally process the state and the main table.
6563 6709 //
6564 6710 if (lazy) {
6565 - fillMainLazily(dsName, ds, selected_eb_scan_m,effectiveBwPerDD_m,e_query_cm);
6711 + fillMainLazily(dsName, ds, selected_eb_scan_m,effectiveBwPerDD_m,e_query_cm,checkDupInts);
6566 6712 }
6567 6713 else {
6568 6714 ​
6569 6715 const MainTable& mainT = ds->getMain();
6570 6716 const StateTable& stateT = ds->getState();
6571 6717
6572 6718 ​
6573 6719 ​
6574 6720 vector<MainRow*> v;
6575 6721 vector<int32_t> mainRowIndex;
6597 6743 ​
6598 6744 unsigned int nMain = v.size();
6599 6745 const VMSData *vmsDataPtr = 0;
6600 6746 ​
6601 6747 // Initialize an UVW coordinates engine.
6602 6748 UvwCoords uvwCoords(ds);
6603 6749
6604 6750 ostringstream oss;
6605 6751 EnumSet<AtmPhaseCorrection> es_query_ap_uncorrected;
6606 6752 es_query_ap_uncorrected.fromString("AP_UNCORRECTED");
6607 -
6608 - // For each selected main row.
6753 +
6754 + // used in checking for duplicate integrations in the WVR (Radiometer) case
6755 + // This holds the most recent last integration time for each configDescriptionId - but only for Radiometer data.
6756 + map<Tag, double> lastTimeMap;
6757 + // for debugging - to report on what uid the lastTime being compared came from
6758 + // map<Tag, string> lastTimeUIDMap;
6759 +​
6609 6760 for (unsigned int i = 0; i < nMain; i++) {
6610 6761 try {
6611 6762 // What's the processor for this Main row ?
6612 6763 Tag cdId = v[i]->getConfigDescriptionId();
6613 6764 ConfigDescriptionTable& cT = ds->getConfigDescription();
6614 6765 ConfigDescriptionRow* cR = cT.getRowByKey(cdId);
6615 6766 Tag pId = cR->getProcessorId();
6616 6767 ​
6617 6768 ​
6618 6769 ProcessorType processorType = ds->getProcessor().getRowByKey(pId)->getProcessorType();
6649 6800 }
6650 6801 ​
6651 6802 if (processorType == RADIOMETER) {
6652 6803 if (!sdmBinData.acceptMainRow(v[i])) {
6653 6804 infostream.str("");
6654 6805 infostream <<"No data retrieved in the Main row #" << mainRowIndex[i] << " (" << sdmBinData.reasonToReject(v[i]) <<")" << endl;
6655 6806 info(infostream.str());
6656 6807 continue;
6657 6808 }
6658 6809 vmsDataPtr = sdmBinData.getDataCols();
6659 -
6810 + bool skipFirstIntegration = checkDupInts && lastTimeMap[cdId] == vmsDataPtr->v_time[0];
6811 + unsigned int skipValues = 0;
6812 + // useful just for debugging
6813 + if (skipFirstIntegration) {
6814 + // work out the number of elements in this first row
6815 + // possibly this is available reliably elsewhere, but thats' not obvious
6816 + // this should be rare and only a few elements should tested, so relatively quick
6817 + bool firstIntegration = True;
6818 + while (firstIntegration && skipValues < vmsDataPtr->v_time.size()) {
6819 + firstIntegration = lastTimeMap[cdId] == vmsDataPtr->v_time[skipValues];
6820 + if (firstIntegration) skipValues++;
6821 + }
6822 + if (debug) {
6823 + cout << "Duplicate time seen in Row : " << i
6824 + << " cdId : " << cdId
6825 + << " " << v[i]->getDataUID().getEntityId().toString()
6826 + // << " duplicates time at end of " << lastTimeUIDMap[cdId]
6827 + << " time size " << vmsDataPtr->v_time.size()
6828 + << " antennaId1 size " << vmsDataPtr->v_antennaId1.size()
6829 + << " skipValues : " << skipValues
6830 + << endl;
6831 + }
6832 + }
6833 +​
6834 + lastTimeMap[cdId] = vmsDataPtr->v_time[vmsDataPtr->v_time.size()-1];
6835 + // lastTimeUIDMap[cdId] = v[i]->getDataUID().getEntityId().toString();
6836 +
6660 6837 fillMain(
6661 6838 v[i],
6662 6839 sdmBinData,
6663 6840 vmsDataPtr,
6664 6841 uvwCoords,
6665 6842 effectiveBwPerDD_m,
6666 6843 complexData,
6667 6844 mute,
6668 - ac_xc_per_timestamp);
6845 + ac_xc_per_timestamp,
6846 + skipValues);
6669 6847
6670 6848 infostream.str("");
6671 - infostream << "ASDM Main row #" << mainRowIndex[i] << " produced a total of " << vmsDataPtr->v_antennaId1.size() << " MS Main rows." << endl;
6849 + infostream << "ASDM Main row #" << mainRowIndex[i] << " produced a total of " << (vmsDataPtr->v_antennaId1.size()-skipValues) << " MS Main rows." << endl;
6672 6850 info(infostream.str());
6673 6851 }
6674 6852 else { // Assume we are in front of a Correlator.
6675 6853 // Open its associate BDF.
6676 6854 ​
6677 6855 bool rowOK = sdmBinData.openMainRow(v[i]);
6678 6856 if (!rowOK) {
6679 6857 infostream.str("");
6680 6858 infostream << "No data retrieved in the Main row #" << mainRowIndex[i] << " (" << sdmBinData.reasonToReject(v[i]) << ")" << endl;
6681 6859 info(infostream.str());

Everything looks good. We'll let you know here if there's anything you should know about.

Add shortcut