| 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 [-t -g -mxr val] 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: [-t -g -mxr val] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n"
|
|---|
| 47 | << " [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDiam] \n"
|
|---|
| 48 | << " [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
|
|---|
| 49 | if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
|
|---|
| 50 | cout << "-t -g : Triangular / gaussian beam shape (def=gaussian) \n"
|
|---|
| 51 | << "-mxr val: Max beam correction factor (default=10.) \n "
|
|---|
| 52 | << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
|
|---|
| 53 | << "- convFacLSS: LSS cube conversion factor to mK (milliKelvin) \n"
|
|---|
| 54 | << "- InMapFgnd: Input 3D foreground cube (PPF file name) \n"
|
|---|
| 55 | << "- convFacFgnd: Foreground cube conversion factor to mK (milliKelvin) \n"
|
|---|
| 56 | << "- PixNoiseLevel: White noise level per pixel (mK) (default=0.) \n"
|
|---|
| 57 | << "- D_Dish/Four2DRespTableFile: Dish diameter or 2D (u,v) plane response (PPF file name) \n"
|
|---|
| 58 | << "- CorBeamDiam: Beam correction target dish diameter \n"
|
|---|
| 59 | << " These two parameters are used to correct for beam effect for a \n"
|
|---|
| 60 | << " target beam (independent of frequency) defined by D/Lambda \n"
|
|---|
| 61 | << " DoL = 100 --> beam ~ 35 arcmin (D=30m @ z~0.5 Lambda~30cm) \n"
|
|---|
| 62 | << " default : no beam correction applied \n "
|
|---|
| 63 | << " - NSigSrcThr: Point source cleaning, Nb_Sigmas on stacked 2D temperature \n"
|
|---|
| 64 | << " default (0.) : no point source cleaning, use NSigSrcThr ~ 3..5 \n"
|
|---|
| 65 | << " - P2/P1: 2nd/first degree polynomial fit on ln(Temp) = f(ln(freq)) \n "
|
|---|
| 66 | << " foreground subtraction. default is P2 \n"
|
|---|
| 67 | << "- RecMapFile: output PPF file for reconstructed foreground template \n"
|
|---|
| 68 | << " (Temperature,SpectralIndex) and extracted LSS cube \n"
|
|---|
| 69 | << endl;
|
|---|
| 70 | return 1;
|
|---|
| 71 | }
|
|---|
| 72 | else cout << " calcpk2 -h for detailed usage " << endl;
|
|---|
| 73 | return 2;
|
|---|
| 74 | }
|
|---|
| 75 | Timer tm("calcpk2");
|
|---|
| 76 | int rc = 0;
|
|---|
| 77 | try {
|
|---|
| 78 | bool fggaussian=true; // true -> gaussian beam
|
|---|
| 79 | // decodage argument optionnel
|
|---|
| 80 | bool fgoptarg=true;
|
|---|
| 81 | double maxratio=10.;
|
|---|
| 82 | while (fgoptarg) {
|
|---|
| 83 | string fbo = arg[1];
|
|---|
| 84 | if (fbo=="-t") { fggaussian=false; arg++; narg--; }
|
|---|
| 85 | else if (fbo=="-g") { fggaussian=true; arg++; narg--; }
|
|---|
| 86 | else if (fbo=="-mxr") { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
|
|---|
| 87 | else fgoptarg=false;
|
|---|
| 88 | }
|
|---|
| 89 | if (narg < 6) {
|
|---|
| 90 | cout << " calcpk2/error arguments , applobe -h for help " << endl;
|
|---|
| 91 | return 2;
|
|---|
| 92 | }
|
|---|
| 93 |
|
|---|
| 94 | string inppflss = arg[1];
|
|---|
| 95 | r_4 rfaclss = atof(arg[2]);
|
|---|
| 96 | string inppfsync = arg[3];
|
|---|
| 97 | r_4 rfacsync = atof(arg[4]);
|
|---|
| 98 | string outname = arg[5];
|
|---|
| 99 |
|
|---|
| 100 | double pixsignoise = 0.;
|
|---|
| 101 | bool fgaddnoise=false;
|
|---|
| 102 | if (narg>6) {
|
|---|
| 103 | pixsignoise=atof(arg[6]);
|
|---|
| 104 | if (pixsignoise>1.e-6) fgaddnoise=true;
|
|---|
| 105 | }
|
|---|
| 106 |
|
|---|
| 107 | bool fgcorrbeam=true;
|
|---|
| 108 |
|
|---|
| 109 | bool fgresptbl=false;
|
|---|
| 110 | double DIAMETRE=100.;
|
|---|
| 111 | string resptblname;
|
|---|
| 112 | if (narg>7) {
|
|---|
| 113 | if (isdigit(*arg[7])) {
|
|---|
| 114 | fgresptbl=false;
|
|---|
| 115 | DIAMETRE=atof(arg[7]);
|
|---|
| 116 | }
|
|---|
| 117 | else {
|
|---|
| 118 | resptblname=arg[7];
|
|---|
| 119 | fgresptbl=true;
|
|---|
| 120 | }
|
|---|
| 121 | }
|
|---|
| 122 | double tbeamDiam=0.;
|
|---|
| 123 | if (narg>8) {
|
|---|
| 124 | tbeamDiam=atof(arg[8]);
|
|---|
| 125 | if (tbeamDiam<1.) fgcorrbeam=false;
|
|---|
| 126 | }
|
|---|
| 127 | bool fgclnsrc=true;
|
|---|
| 128 | double nsigsrc=5.;
|
|---|
| 129 | if (narg>9) {
|
|---|
| 130 | nsigsrc=atof(arg[9]);
|
|---|
| 131 | if (nsigsrc<1.e-6) fgclnsrc=false;
|
|---|
| 132 | }
|
|---|
| 133 | bool fgpoly2=true; // true -> soustraction polynome degre 2
|
|---|
| 134 | if ((narg>0)&&(strcmp(arg[10],"P1")==0)) fgpoly2=false;
|
|---|
| 135 | bool fgsavemaps=false;
|
|---|
| 136 | string outmap_ppfname="extlss.ppf";
|
|---|
| 137 | if (narg>11) {
|
|---|
| 138 | outmap_ppfname=arg[11];
|
|---|
| 139 | fgsavemaps=true;
|
|---|
| 140 | }
|
|---|
| 141 |
|
|---|
| 142 | TArray<r_4> maplss, mapsync;
|
|---|
| 143 | const char * tits[2]={"LSS", "Sync/RadioSrc"};
|
|---|
| 144 | for(int ks=0; ks<2; ks++) {
|
|---|
| 145 | string& ppfname=inppflss;
|
|---|
| 146 | r_4 rfac=rfaclss;
|
|---|
| 147 | TArray<r_4>* inmap=&maplss;
|
|---|
| 148 | if (ks==1) { ppfname=inppfsync; rfac=rfacsync; inmap=&mapsync; }
|
|---|
| 149 | cout << "calcpk2[" << ks+1 << "] : reading 3D map " << tits[ks] << " from file " << ppfname
|
|---|
| 150 | << " RenormFactor=" << rfac << endl;
|
|---|
| 151 | PInPersist pin(ppfname);
|
|---|
| 152 | pin >> (*inmap);
|
|---|
| 153 | (*inmap) *= rfac;
|
|---|
| 154 | double mean, sigma;
|
|---|
| 155 | MeanSigma(*inmap, mean, sigma);
|
|---|
| 156 | cout << " ...InMap sizes " << inmap->InfoString() << endl;
|
|---|
| 157 | inmap->Show();
|
|---|
| 158 | cout << " ... Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 159 | }
|
|---|
| 160 |
|
|---|
| 161 | bool smo;
|
|---|
| 162 | if (!maplss.CompareSizes(mapsync,smo) ) {
|
|---|
| 163 | cout << " calcpk2/ERROR sizes " << endl;
|
|---|
| 164 | maplss.Show(); mapsync.Show();
|
|---|
| 165 | return 99;
|
|---|
| 166 | }
|
|---|
| 167 |
|
|---|
| 168 | TArray<r_4> skycube(mapsync);
|
|---|
| 169 | skycube += maplss;
|
|---|
| 170 |
|
|---|
| 171 | if (fgaddnoise) {
|
|---|
| 172 | cout << " calcpk2: adding noise to skycube cube ... " << endl;
|
|---|
| 173 | BeamEffect::AddNoise(skycube, pixsignoise);
|
|---|
| 174 | }
|
|---|
| 175 |
|
|---|
| 176 | double mean, sigma;
|
|---|
| 177 | MeanSigma(skycube, mean, sigma);
|
|---|
| 178 | cout << " input sky cube : Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 179 | tm.Split(" After input ");
|
|---|
| 180 |
|
|---|
| 181 | H21Conversions conv;
|
|---|
| 182 | conv.setFrequency(Freq0MHz);
|
|---|
| 183 | double lambda = conv.getLambda();
|
|---|
| 184 | Four2DResponse arep(2, DIAMETRE/lambda, DIAMETRE/lambda, lambda);
|
|---|
| 185 | Four2DResponse* arep_p=&arep;
|
|---|
| 186 | Four2DRespTable resptbl;
|
|---|
| 187 | if (fgresptbl) {
|
|---|
| 188 | cout << "calcpk2[3.a]: initializing Four2DRespTable from file" << resptblname << endl;
|
|---|
| 189 | resptbl.readFromPPF(resptblname);
|
|---|
| 190 | resptbl.renormalize(1.);
|
|---|
| 191 | arep_p=&resptbl;
|
|---|
| 192 | }
|
|---|
| 193 | else cout << " calcpk2[3.a]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
|
|---|
| 194 | << " DoL=" << DIAMETRE/lambda << " ) " << endl;
|
|---|
| 195 |
|
|---|
| 196 | double DoL = tbeamDiam/lambda;
|
|---|
| 197 | double tbeamarcmin = RadianToDegree(1.22/DoL)*60.;
|
|---|
| 198 | int typcb = (fggaussian)?1:2;
|
|---|
| 199 | // if (fgresptbl) typcb=22;
|
|---|
| 200 | Four2DResponse tbeam(typcb, DoL, DoL );
|
|---|
| 201 |
|
|---|
| 202 | ForegroundCleaner cleaner(*arep_p, tbeam, skycube, maxratio);
|
|---|
| 203 | if (fgcorrbeam) {
|
|---|
| 204 | cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam
|
|---|
| 205 | << " D/Lambda=" << DoL << " -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
|
|---|
| 206 | cleaner.BeamCorrections();
|
|---|
| 207 | }
|
|---|
| 208 | cout << " calcpk2[3.c] : calling cleaner.CleanNegatives() ... " << endl;
|
|---|
| 209 | cleaner.CleanNegatives();
|
|---|
| 210 | if (fgclnsrc) {
|
|---|
| 211 | cout << "calcpk2[3.d] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc << endl;
|
|---|
| 212 | cleaner.CleanPointSources(nsigsrc);
|
|---|
| 213 | }
|
|---|
| 214 |
|
|---|
| 215 | cout << "calcpk2[4] : calling cleaner.extractLSSCube(...) " << endl;
|
|---|
| 216 | TArray<r_4> synctemp, specidx;
|
|---|
| 217 | TArray<r_4> exlss;
|
|---|
| 218 | if (fgpoly2) exlss = cleaner.extractLSSCubeP2(synctemp, specidx);
|
|---|
| 219 | else exlss = cleaner.extractLSSCubeP1(synctemp, specidx);
|
|---|
| 220 |
|
|---|
| 221 | MeanSigma(exlss, mean, sigma);
|
|---|
| 222 | cout << " After cleaning/extractLSS: Mean=" << mean << " Sigma=" << sigma << endl;
|
|---|
| 223 | tm.Split(" After CleanForeground");
|
|---|
| 224 |
|
|---|
| 225 | cout << "calcpk2[5] : computing 3D Fourier coefficients ... " << endl;
|
|---|
| 226 | FFTWServer ffts;
|
|---|
| 227 | TArray< complex<r_4> > four3d;
|
|---|
| 228 | ffts.FFTForward(exlss, four3d);
|
|---|
| 229 | tm.Split(" After FFTForward ");
|
|---|
| 230 |
|
|---|
| 231 | cout << "calcpk2[6] : computing power spectrum ... " << endl;
|
|---|
| 232 | RandomGenerator rg;
|
|---|
| 233 | Four3DPk pkc(four3d, rg);
|
|---|
| 234 | double dkxmpc = DeuxPI/(double)exlss.SizeX()/XCellSizeMpc;
|
|---|
| 235 | double dkympc = DeuxPI/(double)exlss.SizeY()/YCellSizeMpc;
|
|---|
| 236 | double dkzmpc = DeuxPI/(double)exlss.SizeZ()/ZCellSizeMpc;
|
|---|
| 237 | pkc.SetCellSize(dkxmpc, dkympc, dkzmpc);
|
|---|
| 238 |
|
|---|
| 239 | HProf hp = pkc.ComputePk(0.,HPk_NBin);
|
|---|
| 240 |
|
|---|
| 241 | tm.Split(" Done ComputePk ");
|
|---|
| 242 | {
|
|---|
| 243 | cout << "calcpk2[7.a] : writing profile P(k) to " << outname << endl;
|
|---|
| 244 | POutPersist po(outname);
|
|---|
| 245 | po << hp;
|
|---|
| 246 | }
|
|---|
| 247 | if (fgsavemaps) {
|
|---|
| 248 | cout << "calcpk2[7.b] : writing foreground maps and extracted LSS to " << outmap_ppfname << endl;
|
|---|
| 249 | POutPersist pom(outmap_ppfname);
|
|---|
| 250 | pom << PPFNameTag("Tsync") << synctemp;
|
|---|
| 251 | pom << PPFNameTag("async") << specidx;
|
|---|
| 252 | pom << PPFNameTag("extlss") << exlss;
|
|---|
| 253 | }
|
|---|
| 254 |
|
|---|
| 255 | } // End of try bloc
|
|---|
| 256 | catch (PThrowable & exc) { // catching SOPHYA exceptions
|
|---|
| 257 | cerr << " calcpk2.cc: Catched Exception (PThrowable)" << (string)typeid(exc).name()
|
|---|
| 258 | << "\n...exc.Msg= " << exc.Msg() << endl;
|
|---|
| 259 | rc = 99;
|
|---|
| 260 | }
|
|---|
| 261 | catch (std::exception & e) { // catching standard C++ exceptions
|
|---|
| 262 | cerr << " calcpk2.cc: Catched std::exception " << " - what()= " << e.what() << endl;
|
|---|
| 263 | rc = 98;
|
|---|
| 264 | }
|
|---|
| 265 | catch (...) { // catching other exceptions
|
|---|
| 266 | cerr << " calcpk2.cc: some other exception (...) was caught ! " << endl;
|
|---|
| 267 | rc = 97;
|
|---|
| 268 | }
|
|---|
| 269 | cout << " ==== End of calcpk2.cc program Rc= " << rc << endl;
|
|---|
| 270 | return rc;
|
|---|
| 271 | }
|
|---|
| 272 |
|
|---|
| 273 |
|
|---|
![(please configure the [header_logo] section in trac.ini)](/PI/chrome/site/your_project_logo.png){kind=logo}
