| 1 | /* ------------------------ Projet BAORadio --------------------
|
|---|
| 2 | Programme de calcul du spectre de puissance (3D) a partir d'un
|
|---|
| 3 | cube de delta T/T LSS, d'un cube delta T/T LSS synchrotron
|
|---|
| 4 | ou radio-sources, apres ajustement / soustraction d'une loi de
|
|---|
| 5 | puissance en frequence (l'axe Z du tableau doit etre en frequence)
|
|---|
| 6 |
|
|---|
| 7 | R. Ansari , C. Magneville - Juin 2010
|
|---|
| 8 |
|
|---|
| 9 | Usage: calcpk2 InMapLSS convFacLSS InMapSync convFacSync InMapRadioSource convFacRsc OutPkFile
|
|---|
| 10 | [PixNoiseLevel] [Diameter/Four2DRespTableFile] [TargetBeamArcmin] [NSigSrcThr]
|
|---|
| 11 | --------------------------------------------------------------- */
|
|---|
| 12 |
|
|---|
| 13 | #include "machdefs.h"
|
|---|
| 14 | #include "sopnamsp.h"
|
|---|
| 15 | #include <iostream>
|
|---|
| 16 | #include <string>
|
|---|
| 17 | #include <math.h>
|
|---|
| 18 |
|
|---|
| 19 | #include <typeinfo>
|
|---|
| 20 |
|
|---|
| 21 | #include "specpk.h"
|
|---|
| 22 | #include "histats.h"
|
|---|
| 23 | #include "vector3d.h"
|
|---|
| 24 |
|
|---|
| 25 | #include "qhist.h"
|
|---|
| 26 | #include "lobe.h"
|
|---|
| 27 | #include "cubedef.h"
|
|---|
| 28 | #include "fgndsub.h"
|
|---|
| 29 | #include "radutil.h"
|
|---|
| 30 |
|
|---|
| 31 | #include "histinit.h"
|
|---|
| 32 | #include "fftwserver.h"
|
|---|
| 33 | #include "randr48.h"
|
|---|
| 34 |
|
|---|
| 35 | #include "ctimer.h"
|
|---|
| 36 |
|
|---|
| 37 | typedef ThSDR48RandGen RandomGenerator ;
|
|---|
| 38 |
|
|---|
| 39 | //-------------------------------------------------------------------------
|
|---|
| 40 | // ------------------ MAIN PROGRAM ------------------------------
|
|---|
| 41 | //-------------------------------------------------------------------------
|
|---|
| 42 | /* --Fonction-- */
|
|---|
| 43 | int main(int narg, const char* arg[])
|
|---|
| 44 | {
|
|---|
| 45 | if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
|
|---|
| 46 | cout << " Usage: calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n"
|
|---|
| 47 | << " [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDoL] \n"
|
|---|
| 48 | << " [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
|
|---|
| 49 | if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
|
|---|
| 50 | cout << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
|
|---|
| 51 | << "- convFacLSS: LSS cube conversion factor to mK (milliKelvin) \n"
|
|---|
| 52 | << "- InMapFgnd: Input 3D foreground cube (PPF file name) \n"
|
|---|
| 53 | << "- convFacFgnd: Foreground cube conversion factor to mK (milliKelvin) \n"
|
|---|
| 54 | << "- PixNoiseLevel: White noise level per pixel (mK) (default=0.) \n"
|
|---|
| 55 | << "- D_Dish/Four2DRespTableFile: Dish diameter or 2D (u,v) plane response (PPF file name) \n"
|
|---|
| 56 | << "- CorBeamDoL: Beam correction target Diameter/Lambda \n"
|
|---|
| 57 | << " These two parameters are used to correct for beam effect for a \n"
|
|---|
| 58 | << " target beam (independent of frequency) defined by D/Lambda \n"
|
|---|
| 59 | << " DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5) \n"
|
|---|
| 60 | << " default : no beam correction applied \n "
|
|---|
| 61 | << " - NSigSrcThr: Point source cleaning, Nb_Sigmas on stacked 2D temperature \n"
|
|---|
| 62 | << " default : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
|
|---|
| 63 | << " - P2/P1: 2nd/first degree polynomial fit on ln(Temp) = f(ln(freq)) \n "
|
|---|
| 64 | << " foreground subtraction. default is P2 \n"
|
|---|
| 65 | << "- RecMapFile: output PPF file for reconstructed foreground template \n"
|
|---|
| 66 | << " (Temperature,SpectralIndex) and extracted LSS cube \n"
|
|---|
| 67 | << endl;
|
|---|
| 68 | return 1;
|
|---|
| 69 | }
|
|---|
| 70 | else cout << " calcpk2 -h for detailed usage " << endl;
|
|---|
| 71 | return 2;
|
|---|
| 72 | }
|
|---|
| 73 | Timer tm("calcpk2");
|
|---|
| 74 | int rc = 0;
|
|---|
| 75 | try {
|
|---|
| 76 | string inppflss = arg[1];
|
|---|
| 77 | r_4 rfaclss = atof(arg[2]);
|
|---|
| 78 | string inppfsync = arg[3];
|
|---|
| 79 | r_4 rfacsync = atof(arg[4]);
|
|---|
| 80 | string outname = arg[5];
|
|---|
| 81 |
|
|---|
| 82 | double pixsignoise = 0.;
|
|---|
| 83 | bool fgaddnoise=false;
|
|---|
| 84 | if (narg>6) {
|
|---|
| 85 | pixsignoise=atof(arg[6]);
|
|---|
| 86 | if (pixsignoise>1.e-6) fgaddnoise=true;
|
|---|
| 87 | }
|
|---|
| 88 |
|
|---|
| 89 | bool fgcorrbeam=true;
|
|---|
| 90 |
|
|---|
| 91 | bool fgresptbl=false;
|
|---|
| 92 | double DIAMETRE=100.;
|
|---|
| 93 | string resptblname;
|
|---|
| 94 | if (narg>7) {
|
|---|
| 95 | if (isdigit(*arg[7])) {
|
|---|
| 96 | fgresptbl=false;
|
|---|
| 97 | DIAMETRE=atof(arg[7]);
|
|---|
| 98 | }
|
|---|
| 99 | else {
|
|---|
| 100 | resptblname=arg[7];
|
|---|
| 101 | fgresptbl=true;
|
|---|
| 102 | }
|
|---|
| 103 | }
|
|---|
| 104 | double tbeamDoL=0.;
|
|---|
| 105 | if (narg>8) {
|
|---|
| 106 | tbeamDoL=atof(arg[8]);
|
|---|
| 107 | if (tbeamDoL<1.) fgcorrbeam=false;
|
|---|
| 108 | }
|
|---|
| 109 | bool fgclnsrc=true;
|
|---|
| 110 | double nsigsrc=5.;
|
|---|
| 111 | if (narg>9) {
|
|---|
| 112 | nsigsrc=atof(arg[9]);
|
|---|
| 113 | if (nsigsrc<1.e-6) fgclnsrc=false;
|
|---|
| 114 | }
|
|---|
| 115 | bool fgpoly2=true; // true -> soustraction polynome degre 2
|
|---|
| 116 | if ((narg>0)&&(strcmp(arg[10],"P1")==0)) fgpoly2=false;
|
|---|
| 117 | bool fgsavemaps=false;
|
|---|
| 118 | string outmap_ppfname="extlss.ppf";
|
|---|
| 119 | if (narg>11) {
|
|---|
| 120 | outmap_ppfname=arg[11];
|
|---|
| 121 | fgsavemaps=true;
|
|---|
| 122 | }
|
|---|
| 123 |
|
|---|
| 124 | TArray<r_4> maplss, mapsync;
|
|---|
| 125 | const char * tits[2]={"LSS", "Sync/RadioSrc"};
|
|---|
| 126 | for(int ks=0; ks<2; ks++) {
|
|---|
| 127 | string& ppfname=inppflss;
|
|---|
| 128 | r_4 rfac=rfaclss;
|
|---|
| 129 | TArray<r_4>* inmap=&maplss;
|
|---|
| 130 | if (ks==1) { ppfname=inppfsync; rfac=rfacsync; inmap=&mapsync; }
|
|---|
| 131 | cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname
|
|---|
| 132 | << " RenormFactor=" << rfac << endl;
|
|---|
| 133 | PInPersist pin(ppfname);
|
|---|
| 134 | pin >> (*inmap);
|
|---|
| 135 | (*inmap) *= rfac;
|
|---|
| 136 | double mean, sigma;
|
|---|
| 137 | MeanSigma(*inmap, mean, sigma);
|
|---|
| 138 | cout << " ...InMap sizes " << inmap->InfoString() << endl;
|
|---|
| 139 | inmap->Show();
|
|---|
| 140 | cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 141 | }
|
|---|
| 142 |
|
|---|
| 143 | bool smo;
|
|---|
| 144 | if (!maplss.CompareSizes(mapsync,smo) ) {
|
|---|
| 145 | cout << " calcpk2/ERROR sizes " << endl;
|
|---|
| 146 | maplss.Show(); mapsync.Show();
|
|---|
| 147 | return 99;
|
|---|
| 148 | }
|
|---|
| 149 |
|
|---|
| 150 | TArray<r_4> skycube(mapsync);
|
|---|
| 151 | skycube += maplss;
|
|---|
| 152 |
|
|---|
| 153 | if (fgaddnoise) {
|
|---|
| 154 | cout << " calcpk2: adding noise to skycube cube ... " << endl;
|
|---|
| 155 | BeamEffect::AddNoise(skycube, pixsignoise);
|
|---|
| 156 | }
|
|---|
| 157 |
|
|---|
| 158 | double mean, sigma;
|
|---|
| 159 | MeanSigma(skycube, mean, sigma);
|
|---|
| 160 | cout << " input sky cube : Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 161 | tm.Split(" After input ");
|
|---|
| 162 |
|
|---|
| 163 | H21Conversions conv;
|
|---|
| 164 | conv.setFrequency(Freq0MHz);
|
|---|
| 165 | double lambda = conv.getLambda();
|
|---|
| 166 | Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
|
|---|
| 167 | Four2DResponse* arep_p=&arep;
|
|---|
| 168 | Four2DRespTable resptbl;
|
|---|
| 169 | if (fgresptbl) {
|
|---|
| 170 | cout << "calcpk2[3.a]: initializing Four2DRespTable from file" << resptblname << endl;
|
|---|
| 171 | resptbl.readFromPPF(resptblname);
|
|---|
| 172 | arep_p=&resptbl;
|
|---|
| 173 | }
|
|---|
| 174 | else cout << " calcpk2[3.a]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
|
|---|
| 175 | << " DoL=" << DIAMETRE/lambda << " ) " << endl;
|
|---|
| 176 |
|
|---|
| 177 | double DoL = tbeamDoL;
|
|---|
| 178 | double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
|
|---|
| 179 | int typcb = 2;
|
|---|
| 180 | // if (fgresptbl) typcb=22;
|
|---|
| 181 | Four2DResponse tbeam(typcb, DoL, DoL );
|
|---|
| 182 |
|
|---|
| 183 | ForegroundCleaner cleaner(*arep_p, tbeam, skycube);
|
|---|
| 184 | if (fgcorrbeam) {
|
|---|
| 185 | cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam D/Lambda=" << DoL
|
|---|
| 186 | << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
|
|---|
| 187 | cleaner.BeamCorrections();
|
|---|
| 188 | }
|
|---|
| 189 | cout << " calcpk2[3.c] : calling cleaner.CleanNegatives() ... " << endl;
|
|---|
| 190 | cleaner.CleanNegatives();
|
|---|
| 191 | if (fgclnsrc) {
|
|---|
| 192 | cout << "calcpk2[3.d] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc << endl;
|
|---|
| 193 | cleaner.CleanPointSources(nsigsrc);
|
|---|
| 194 | }
|
|---|
| 195 |
|
|---|
| 196 | cout << "calcpk2[4] : calling cleaner.extractLSSCube(...) " << endl;
|
|---|
| 197 | TArray<r_4> synctemp, specidx;
|
|---|
| 198 | TArray<r_4> exlss;
|
|---|
| 199 | if (fgpoly2) exlss = cleaner.extractLSSCubeP2(synctemp, specidx);
|
|---|
| 200 | else exlss = cleaner.extractLSSCubeP1(synctemp, specidx);
|
|---|
| 201 |
|
|---|
| 202 | MeanSigma(exlss, mean, sigma);
|
|---|
| 203 | cout << " After cleaning/extractLSS: Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 204 | tm.Split(" After CleanForeground");
|
|---|
| 205 |
|
|---|
| 206 | cout << "calcpk2[5] : computing 3D Fourier coefficients ... " << endl;
|
|---|
| 207 | FFTWServer ffts;
|
|---|
| 208 | TArray< complex<r_4> > four3d;
|
|---|
| 209 | ffts.FFTForward(exlss, four3d);
|
|---|
| 210 | tm.Split(" After FFTForward ");
|
|---|
| 211 |
|
|---|
| 212 | cout << "calcpk2[6] : computing power spectrum ... " << endl;
|
|---|
| 213 | RandomGenerator rg;
|
|---|
| 214 | Four3DPk pkc(four3d, rg);
|
|---|
| 215 | double dkxmpc = DeuxPI/(double)exlss.SizeX()/XCellSizeMpc;
|
|---|
| 216 | double dkympc = DeuxPI/(double)exlss.SizeY()/YCellSizeMpc;
|
|---|
| 217 | double dkzmpc = DeuxPI/(double)exlss.SizeZ()/ZCellSizeMpc;
|
|---|
| 218 | pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
|
|---|
| 219 |
|
|---|
| 220 | HProf hp = pkc.ComputePk(0.,256);
|
|---|
| 221 |
|
|---|
| 222 | tm.Split(" Done ComputePk ");
|
|---|
| 223 | {
|
|---|
| 224 | cout << "calcpk2[7.a] : writing profile P(k) to " << outname << endl;
|
|---|
| 225 | POutPersist po(outname);
|
|---|
| 226 | po << hp;
|
|---|
| 227 | }
|
|---|
| 228 | if (fgsavemaps) {
|
|---|
| 229 | cout << "calcpk2[7.b] : writing foreground maps and extracted LSS to " << outmap_ppfname << endl;
|
|---|
| 230 | POutPersist pom(outmap_ppfname);
|
|---|
| 231 | pom << PPFNameTag("Tsync") << synctemp;
|
|---|
| 232 | pom << PPFNameTag("async") << specidx;
|
|---|
| 233 | pom << PPFNameTag("extlss") << exlss;
|
|---|
| 234 | }
|
|---|
| 235 |
|
|---|
| 236 | } // End of try bloc
|
|---|
| 237 | catch (PThrowable & exc) { // catching SOPHYA exceptions
|
|---|
| 238 | cerr << " calcpk2.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
|
|---|
| 239 | << "\n...exc.Msg= " << exc.Msg() << endl;
|
|---|
| 240 | rc = 99;
|
|---|
| 241 | }
|
|---|
| 242 | catch (std::exception & e) { // catching standard C++ exceptions
|
|---|
| 243 | cerr << " calcpk2.cc: Catched std::exception " << " - what()= " << e.what() << endl;
|
|---|
| 244 | rc = 98;
|
|---|
| 245 | }
|
|---|
| 246 | catch (...) { // catching other exceptions
|
|---|
| 247 | cerr << " calcpk2.cc: some other exception (...) was caught ! " << endl;
|
|---|
| 248 | rc = 97;
|
|---|
| 249 | }
|
|---|
| 250 | cout << " ==== End of calcpk2.cc program Rc= " << rc << endl;
|
|---|
| 251 | return rc;
|
|---|
| 252 | }
|
|---|
| 253 |
|
|---|
| 254 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
