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