Changeset 4022 in Sophya


Ignore:
Timestamp:
Sep 28, 2011, 5:13:51 PM (13 years ago)
Author:
ansari
Message:

modifs pour pouvoir imposer la moyenne en temp des plans X,Y des cubes lors de l'extraction du signal HI, Reza 28/9/2011

Location:
trunk/Cosmo/RadioBeam
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed

  • trunk/Cosmo/RadioBeam/applobe.cc

    r3991 r4022  
    5858
    5959  if ((narg < 3)||(strcmp(arg[1],"-h")==0)) {
    60     cout << "Usage: applobe [-t -g -fib -mxr val] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n"
     60    cout << "Usage: applobe [-t -g -fib -kzf/-nzf -mxr val] Diameter/Four2DRespTableFile In3DPPFName Out3DPPFName \n"
    6161         << "               [TargetBeamDiam] [ResmapleFactor=0.5,0.333...] \n" << endl;
    6262    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
    6363      cout << "   -t -g : Triangular / gaussian beam shape (def=gaussian) \n"
    6464           << "   -fib :  Application of a fixed (freq.independent) lobe dish-triangle or gaussian \n"
     65           << "   -kzf -nzf :  Keep (default) or Not zero space frequency when applying lobes (BeamEffect class) \n"
    6566           << "   -mxr val: Max beam correction factor (default=10.) \n"
    6667           << "   Diameter/Four2DRespTableFile : dish diameter or 2D response PPF file name\n"
    … …  
    7980  bool fixedbeam=false;  // true -> apply freq. independent beam
    8081  double maxratio=10.;   // valeur max du rapport des lobes lors de la correction de lobe
     82  bool preservefreq0=true; // true -> keep zero frequency when appyling lobe
    8183
    8284  // decodage argument optionnel
    … …  
    8789    else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
    8890    else if (fbo=="-fib")  { fixedbeam=true; arg++; narg--; }
     91    else if (fbo=="-kzf")  { preservefreq0=true; arg++; narg--; }
     92    else if (fbo=="-nzf")  { preservefreq0=false; arg++; narg--; }
    8993    else if (fbo=="-mxr")  { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
    9094    else fgoptarg=false;
    … …  
    164168    else cout << " applobe[2.b]: Four2DResponse ( Diameter=" << DIAMETRE << " Lambda= " << lambda
    165169              << " DoL=" << DIAMETRE/lambda << " ) " << endl;
    166     BeamEffect beam(*fresp_p);
     170    BeamEffect beam(*fresp_p,preservefreq0);
    167171
    168172    if (fgcorrbeam) {

  • trunk/Cosmo/RadioBeam/calcpk2.cc

    r3989 r4022  
    4444{
    4545  if ( (narg<6)||((narg>1)&&(strcmp(arg[1],"-h")==0)) ) {
    46     cout << " Usage: [-t -g -mxr val] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n"
     46    cout << " Usage: [-t -g -mxr val -fmt T0,alpha] calcpk2 InMapLSS convFacLSS InMapFgnd convFacFgnd OutPkFile \n"
    4747         << "        [PixNoiseLevel] [D_Dish/Four2DRespTableFile CorBeamDiam] \n"
    4848         << "        [NSigSrcThr] [P2/P1] [RecMapFile] " << endl;
    4949    if ((narg>1)&&(strcmp(arg[1],"-h")==0)) {
    5050      cout << "-t -g : Triangular / gaussian beam shape (def=gaussian) \n"
     51           << "-fmt T0,alpha: fix mean slice temp. according to power law T0,alpha \n "
    5152           << "-mxr val: Max beam correction factor (default=10.) \n "
    5253           << "- InMapLSS: Input 3D LSS cube (PPF file name) \n "
    … …  
    7980    double maxratio=10.;   // valeur max du rapport des lobes lors de la correction de lobe
    8081
     82    bool fgfix_mXYtemp=false;
     83    double T0_pl=1.;
     84    double alpha_pl=0.;
    8185    // decodage argument optionnel
    8286    bool fgoptarg=true;
    … …  
    8690      else if (fbo=="-g")  { fggaussian=true; arg++; narg--; }
    8791      else if (fbo=="-mxr")  { arg++; maxratio=atof(arg[1]); arg++; narg-=2; }
     92      else if (fbo=="-fmt") 
     93        { fgfix_mXYtemp=true; arg++; sscanf(arg[1],"%lg,%lg",&T0_pl,&alpha_pl); arg++; narg-=2; }
    8894      else fgoptarg=false;
    8995    }
    … …  
    202208
    203209    ForegroundCleaner  cleaner(*arep_p, tbeam, skycube, maxratio);
     210    if (fgfix_mXYtemp) {
     211      cout << "calcpk2[3.b] : calling cleaner.FixMeanXYTemp(T0=" << T0_pl << ",alpha=" << alpha_pl << ")" << endl;
     212      cleaner.FixMeanXYTemp(T0_pl,alpha_pl);
     213    }
    204214    if (fgcorrbeam) {
    205       cout << "calcpk2[3.b] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam 
     215      cout << "calcpk2[3.c] : calling cleaner.BeamCorrections() for target beam Diameter=" << tbeamDiam 
    206216           << " D/Lambda=" << DoL  << "  -> arcmin " << tbeamarcmin << " TypDishResp=" << typcb << endl;
    207217      cleaner.BeamCorrections();
    208218    }
    209     cout << " calcpk2[3.c] : calling cleaner.CleanNegatives() ... " << endl;
     219    cout << " calcpk2[3.d] : calling cleaner.CleanNegatives() ... " << endl;
    210220    cleaner.CleanNegatives();
    211221    if (fgclnsrc) {
    212       cout << "calcpk2[3.d] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
     222      cout << "calcpk2[3.e] : calling cleaner.CleanPointSources() with threshold NSigma=" << nsigsrc  << endl;
    213223      cleaner.CleanPointSources(nsigsrc);
    214224    }

  • trunk/Cosmo/RadioBeam/fgndsub.cc

    r3989 r4022  
    1212
    1313#include "ctimer.h"
     14
     15/* --Methode-- */
     16PowerLawChecker::PowerLawChecker(TArray< TF >& skycube)
     17  : skycube_(skycube)
     18{
     19}
     20/* --Methode-- */
     21void PowerLawChecker::CheckXYMean()
     22{
     23  // double freq0 : Frequence premier index en k (MHz)
     24  // double dfreq :   // largeur en frequence de chaque plan (Mhz) 
     25  double freq0_ = Freq0MHz;
     26  double dfreq_ = FreqSizeMHz/(double)NFreq;
     27
     28  double tempfirst,templast;
     29  r_8 s1, sx, sx2, sy, sxy;
     30  s1=sx=sx2=sy=sxy=0.;
     31  double lnf0=0.;
     32  for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
     33    double lnf=log((double)k*dfreq_+freq0_);
     34    double ttot = Mean(skycube_(Range::all(), Range::all(), Range(k)));
     35    if (k==0) { tempfirst=ttot;  lnf0=lnf; }
     36    if (k==skycube_.SizeZ()-1) templast=ttot;
     37    if (ttot < 1.e-5) continue;
     38    double lntt=log(ttot);
     39    s1+=1.;  sx+=lnf;  sx2+=(lnf*lnf);
     40    sy+=lntt;    sxy+=(lnf*lntt);
     41  }
     42  double beta = (sx*sxy-sx2*sy)/(sx*sx-s1*sx2);
     43  double alpha = (s1*sxy-sx*sy)/(s1*sx2-sx*sx);
     44  double T0 = exp(beta+alpha*lnf0);
     45 
     46  cout << " PowerLawChecker::CheckMean() meanTemp(0 ...last) " << tempfirst << " ... "
     47       << templast << endl;
     48  bool fgnan = false;
     49  if (!isfinite(alpha)||(!isfinite(beta))) {
     50    cout << "ePowerLawChecker::CheckMean()  Not finite alpha, beta " << endl;
     51    fgnan = true;
     52  }
     53  cout << "PowerLawChecker::CheckMean() - T0=" << T0 << " alpha=" << alpha << "(beta=" << beta << ")" << endl;
     54  return;
     55}
    1456
    1557/* --Methode-- */
    … …  
    3476  beam.Correct2RefLobe(tbeam_, skycube_, dx_, dy_, freq0_, dfreq_, maxratio_);
    3577  cout << " ForegroundCleaner::BeamCorrections() done " << endl;
     78}
     79
     80/* --Methode-- */
     81int ForegroundCleaner::FixMeanXYTemp(double T0, double alpha)
     82{
     83  cout << "ForegroundCleaner::FixMeanXYTemp(T0=" << T0 << ",alpha=" << alpha << ")" << endl;
     84  double lnf0=log(freq0_);
     85  sa_size_t modprt=skycube_.SizeZ()/12;
     86  for(sa_size_t k=0; k<skycube_.SizeZ(); k++)  {
     87    double lnf=log((double)k*dfreq_+freq0_);
     88    double fittedtemp = T0*exp(alpha*(lnf-lnf0));
     89    TArray<TF> slice = skycube_(Range::all(), Range::all(), Range(k));
     90    TF deltatemp = (TF)(fittedtemp-(double)Mean(slice));
     91    slice += deltatemp;
     92    if (k%modprt == 0)
     93      cout << "FixMeanXYTemp[k=" << k << " MeanXYTemp=" << fittedtemp-deltatemp << " -> " << fittedtemp
     94           << " (DeltaTemp=" << deltatemp << ")" << endl;
     95  }
     96  cout << "ForegroundCleaner::FixMeanXYTemp done" << endl;
     97  return 0;
    3698}
    3799

  • trunk/Cosmo/RadioBeam/fgndsub.h

    r3986 r4022  
    2323#endif
    2424
     25//------------------------------------------------------------------------------------
     26// Cacarterisation du comportement en loi de puissance des cubes de temperature
     27class PowerLawChecker {
     28public:
     29  PowerLawChecker(TArray< TF >& skycube);
     30  void CheckXYMean();
     31
     32  TArray< TF > skycube_;
     33};
     34
     35//------------------------------------------------------------------------------------
     36// Classe implementant la soustraction des avant-plans sous forme de loi de puissance
    2537class ForegroundCleaner {
    2638public:
    2739  ForegroundCleaner(Four2DResponse& arrep, Four2DResponse& tbeam, TArray< TF >& skycube, double maxratio=10.);
    2840  void BeamCorrections();
     41  int FixMeanXYTemp(double T0, double alpha);
    2942  int CleanNegatives(TF seuil=1.e-6);
    3043  int CleanPointSources(double nsigmas=5.);
    … …  
    3952  double maxratio_;
    4053  double dx_, dy_;   // taille des pixels (radians) de skycube
    41   double freq0_, dfreq_;    // 1ere frequence et bin en frequence de skycube_;
    42  
     54  double freq0_, dfreq_;    // 1ere frequence et bin en frequence de skycube_; 
    4355};
    4456
     57
    4558#endif

  • trunk/Cosmo/RadioBeam/makefile

    r3825 r4022  
    11############################################################
    22##  makefile for interferometer response
    3 ##  R.Ansari - Avril 2010  
     3##  R.Ansari - Avril 2010 - Sep 2011
    44############################################################
    55
    66include $(SOPHYABASE)/include/sophyamake.inc
    77
    8 all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe gsm2cube
     8all : pknoise repicon calcpk calcpk2 syncube srcat2cube tjyk applobe gsm2cube ckpowl
    99
    1010clean :
    … …  
    4141applobe : Objs/applobe
    4242        echo 'makefile : applobe made'
     43
     44ckpowl : Objs/ckpowl
     45        echo 'makefile : ckpowl made'
    4346
    4447### programme pknoise
    … …  
    105108        $(CXXCOMPILE) -o Objs/gsm2cube.o gsm2cube.cc
    106109
     110### programme ckpowl
     111Objs/ckpowl : Objs/ckpowl.o $(PKGOLIST)
     112        $(CXXLINK) -o Objs/ckpowl Objs/ckpowl.o $(PKGOLIST) $(SOPHYAEXTSLBLIST)
     113
     114Objs/ckpowl.o : ckpowl.cc $(PKGHLIST)
     115        $(CXXCOMPILE) -o Objs/ckpowl.o ckpowl.cc
     116
    107117### les classes / fonctions
    108118Objs/fgndsub.o : fgndsub.cc $(PKGHLIST)

  • trunk/Cosmo/RadioBeam/sensfgnd21cm.tex

    r4014 r4022  
    128128cm emission. Such a 3D matter distribution map can be used to test the Cosmological model and to constrain the Dark Energy
    129129properties or its equation of state. A novel approach, called intensity mapping can be used to map the \HI distribution,
    130 using radio interferometers with large instanteneous field of view and waveband.}
     130using radio interferometers with large instantaneous field of view and waveband.}
    131131 % aims heading (mandatory)
    132132  { In this paper, we study the sensitivity of different radio interferometer configurations, or multi-beam
    … …  
    352352\end{equation}
    353353Introducing the \HI mass fraction relative to the total baryon mass $\gHI$, the
    354 neutral hydrogen number density relative fluctuations can be written as, and the corresponding
    355 21 cm emission temperature can be written as:
     354neutral hydrogen number density and the corresponding 21 cm emission temperature
     355can be written as a function of \HI relative density fluctuations:
    356356\begin{eqnarray}
    357357\etaHI (\vec{\Theta}, z(\lambda) ) & = & \gHIz \times \Omega_B  \frac{\rho_{crit}}{m_{H}}  \times
Note: See TracChangeset for help on using the changeset viewer.