source: Sophya/trunk/Cosmo/SimLSS/cmvtstagn.cc

// Distribution des flux des AGN dans un pixel de taille donne
#include "sopnamsp.h"
#include "machdefs.h"
#include <iostream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <math.h>
#include <unistd.h>
#include "sophyainit.h"
#include "timing.h"
#include "ntuple.h"

#include "geneutils.h"
#include "agnjackson.h"
#include "srandgen.h"

void usage(void);
void usage(void)
{
cout<<"cmvtstagn "<<endl
 <<" -a angsol: taille pixel (en sr si >0 ou arcmin^2 si <0)"<<endl
 <<" -x nbin,xmin,xmax: def histo des log10(S Jy)/pixel"<<endl
 <<"              si nbin<0 on a de <S>+xmin , <S>-xmax"<<endl
 <<" -n npix: nombre de tirages de pixel"<<endl;
}


int main(int narg,char *arg[])
{
SophyaInit();
 double angsol = 1.;  // sr
 double lxmin=-3.,lxmax=4.; int nbinlx = -100;
 int npix = 100;

 // --- Decodage arguments
 char c;
 while((c = getopt(narg,arg,"ha:x:n:")) != -1) {
  switch (c) {
  case 'a' :
    sscanf(optarg,"%lf",&angsol);
    break;
  case 'x' :
    sscanf(optarg,"%d,%lf,%lf",&nbinlx,&lxmin,&lxmax);
    break;
  case 'n' :
    sscanf(optarg,"%d",&npix);
    break;
  case 'h' :
  default :
    usage();
    return -1;
  }
 }

 // --- Class de distribution d'AGN
 double sr2amin = AngSol(M_PI/180./60./2.,M_PI/180./60./2.);
 if(angsol<0) angsol *= -sr2amin;
 cout<<"angsol = "<<angsol<<" sr -> "<<angsol/sr2amin<<" arcmin^2"<<endl;
 AGNJackson agn;
 agn.Print();

 cout<<endl;
 vector<double> xjack, yjack; agn.OrigJack(xjack,yjack);
 Histo hdndls0 = agn.dNdlS();
 Histo hdndls(hdndls0); hdndls *= angsol;
 Histo htirls = agn.TirL10S();
 double nobjang = agn.NObjAng()*angsol;
 double fluxang = agn.FluxAng()*angsol;
 cout<<"Nombre moyen d\'objets dans notre pixel: "<<nobjang<<endl;
 cout<<"Flux moyen dans notre pixel: "<<fluxang<<" Jy"<<endl;

 // --- Distribution sur angsol
 // Histo du nombre d'objets / angsol
 cout<<endl;
 double mu = nobjang; if(mu<=0.) mu=1.;
 double xmin = mu-5.*sqrt(mu); if(xmin<0.) xmin=0.; xmin=floor(xmin);
 double xmax = mu+7.*sqrt(mu); xmax=floor(xmax+1.);
 long nb = long(xmax-xmin+1.01); while(nb>250) nb /= 2;
 // Histo du nombre d'objets/pixel
 Histo hnobj(xmin,xmax,nb); hnobj.ReCenterBin();
 cout<<"hnobj: "; hnobj.Show();
 // Histo de check du tirage
 Histo hlflux(htirls); hlflux.Zero();
 // Histo du flux/pixel
 if(nbinlx==0 || nbinlx==1) nbinlx = -100;
 if(nbinlx<0) {
   nbinlx = -nbinlx;
   lxmin = log10(fluxang) + lxmin;
   lxmax = log10(fluxang) + lxmax;
   cout<<"Recomputing: lxmin="<<lxmin<<" lxmax="<<lxmax<<" nbin="<<nbinlx<<endl;
 }
 Histo hlfang(lxmin,lxmax,nbinlx); hlfang.Zero();
 cout<<"hlfang: "; hlfang.Show();

 // Tirage aleatoire
 cout<<endl;
 if(npix>1) {
   long l = npix/10; if(l<=0) l=1;
   cout<<"Making "<<npix<<" trials"<<endl;
   unsigned long nzero=0, ntir=0;
   double sumpix=0., sumpix2=0., sumn=0., sumn2=0.;
   for(long i=0;i<npix;i++) {
     if(i%l==0) cout<<"...trial: "<<i<<endl;
     double n = PoissonRand(mu,10.);
     hnobj.Add(n);
     sumn += n; sumn2 += n*n;
     double fsum = 0.; unsigned long nf = 0;
     for(unsigned long j=0;j<n;j++) {
       double f = agn.FluxJY();
       hlflux.Add(log10(f));
       fsum += f;
       nf++;
     }
     ntir += nf;
     if(nf==0) {
       nzero++;
     } else {
       sumpix += fsum; sumpix2 += fsum*fsum;
       hlfang.Add(log10(fsum));
     }
   }
   sumn /= (double)npix;
   sumn2 = sumn2/(double)npix - sumn*sumn;
   sumpix /= (double)npix;
   sumpix2 = sumpix2/(double)npix - sumpix*sumpix;
   cout<<"Number of try: "<<ntir<<endl;
   cout<<"Number of pixels with no flux: "<<nzero<<" / "<<npix<<endl;
   cout<<"Mean object number/pix: "<<sumn<<" (check "<<(double)ntir/(double)npix
       <<") s^2="<<sumn2<<" s="<<sqrt(fabs(sumn2))<<endl;
   cout<<"Mean flux/pix: "<<sumpix<<" s^2="<<sumpix2<<" s="<<sqrt(fabs(sumpix2))<<endl;
   // Info de l'histo des flux
   cout<<"hlfang: "; hlfang.Show();
   // normalisation de hlflux
   hlflux *= hdndls.BinWidth()/hlflux.BinWidth()*hdndls.Sum()/hlflux.Sum();
 }

 // --- Ecriture sur fichier ppf
 string tag = "cmvtstagn.ppf";
 POutPersist pos(tag);

 TVector<r_8> vxjack(xjack), vyjack(yjack);
 tag = "xjack"; pos.PutObject(vxjack,tag);
 tag = "yjack"; pos.PutObject(vyjack,tag);

 tag = "dndls0"; pos.PutObject(hdndls0,tag);
 tag = "dndls"; pos.PutObject(hdndls,tag);

 tag = "tirls"; pos.PutObject(htirls,tag);
 tag = "hlflux"; pos.PutObject(hlflux,tag);

 tag = "nobj"; pos.PutObject(hnobj,tag);
 tag = "hlfang"; pos.PutObject(hlfang,tag);

 return 0;
}


/*
openppf cmvtstagn.ppf

del ntjack
c++exec const char *vname[2] = {"x","y"}; double xnt[2]; \
NTuple ntjack(2,vname); \
for(int i=0;i<xjack.Size();i++) \
  {xnt[0]=xjack(i); xnt[1]=yjack(i); ntjack.Fill(xnt);} \
KeepObj(ntjack); cout<<"OK";

# Jackson initial (Jy^1.5/sr)
n/plot ntjack.y%x ! ! "nsta connectpoints"
vecplot xjack yjack "same red fcirclemarker7"

# Check de ce qu'on a extrapole
set yxt "y*pow(10.,-1.5*x)*M_LN10"
n/plot dndls0.log10(val)%x ! ! "nsta connectpoints"
n/plot ntjack.log10($yxt)%x ! ! "nsta red fcirclemarker7 same"

# la pdf : histo en 1/log10(Jy)/sr
disp dndls0 "nsta"
n/plot ntjack.$yxt%x ! ! "nsta red fcirclemarker7 same"

# la pdf : histo en 1/log10(Jy)/angsol
disp dndls "nsta"

# la pdf en flux: dn/dlog10(S)*S en microJy/log10(Jy)/sr
n/plot dndls.val*pow(10.,x)*1e6%x ! ! "nsta connectpoints"

# la fonction de tirage en 1/Jy/(angsol sr)
disp tirls "nsta"

# Check tirages
disp dndls "nsta"
disp hlflux "nsta same red"

# Resultat: nombre d'objets (poisson) par pixel
disp nobj

# Resultat: distribution du flux par pixel
disp hlfang

 */