source: Sophya/trunk/Cosmo/RadioBeam/treccyl.cc@ 3803

#include "sopnamsp.h"
#include "machdefs.h"
#include <math.h>
#include <iostream>
#include <typeinfo>

#include "tvector.h"
#include "srandgen.h"
#include "fioarr.h"
#include "sopemtx.h"
#include "pexceptions.h"
#include "matharr.h"

#include "sambainit.h"

// #include "tarrinit.h"

#include "timing.h"
#include "datacards.h"
#include <dvlist.h>

#include "multicyl.h"
#include "mbeamcyl.h"
#define LENGTH 1024

/* 
   Projet BAORadio / HSHS
   Programme de simulation pour reconstruction de lobe radio.
   programme principal de test

   R. Ansari - LAL      Jan 2007

*/

// Declaration des fonctions de ce fichier
static int test1cyl(string& ppfname);
static int testmulticyl(string& ppfname);
int ReadParam(const char* fileName);

//-----------------------------------------------------------
// -------------- Parametres de simulation  -----------------
//-----------------------------------------------------------
static double tClock = 2.; // should come from param file !!!!
static double cLight=0.3;       // in 1E9 m/s
//static double tClock = 1.; // should come from param file !!!!
//static double cLight=1.;      // in 1E9 m/s
//
static int MR = 256;  // Nombre de recepteur
static int NE = 64;  // Nombre d'echantillon en temps;
static double freq0 = 2.;  // frequence de base
static double da = 0.25;     // pas des antennes le long du cylindre
// ATTENTION : les parametres suivants sont relies a MR/da
static double maxangX = M_PI/3.; // angle max en X ( +/- )
static double maxangY = M_PI/60.; // angle max en Y ( +/- )
static int halfNY;
static int NX;
static int nsrcmax = 50;  // Nb total de sources - en un plan

static double snoise = 1.0;  // sigma du bruit 
static double tjit = 0.05;   // sigma du jitter en temps
static double tos = 0.02;    // sigma des offsets en temps
static double gmean = 1.;    // gain moyen
static double gsig = 0.;     // sigma des gains
static int nantgz = 0;       // nb d'antennes morts (-> gain=0)
static int prtlevel = 0;     // niveau de print 

static int nCyl;
static double xCyl[1000];
static double yCyl[1000];
//-----------------------------------------------------------


/* --------------------------------------------------------
  Le main programme de test des classes de reconstruction
  multilobe radio - R. Ansari , Sep06 -- 2007
  --------------------------------------------------------- */
   
int main(int narg, char* arg[])
{

  SophyaInit();
  InitTim();   // Initializing the CPU timer
  ReadParam("telescope.in");
  cout <<"MR="<< MR <<" NE="<<NE<<" freq0="<<freq0<<" "<<da<<" "<<maxangX <<endl;
  cout << maxangY<<" "<<nsrcmax <<" "<< snoise<<" "<< tjit<<" "<< tos<<" "<<gmean <<endl;
  cout << gsig<<" "<<nantgz <<" "<< prtlevel<<endl;
//  return 1;

  string ppfname = "treccyl.ppf";
  int act = 1;
//  int ncyl = 5; 
  if (narg < 2) {
    cout << "Usage: treccyl act ppfname \n"
       << " -act= X ou XY \n"
       << " -ppfname=  treccyl.ppf par defaut" << endl;
    return 1;
  }
  if (strcmp(arg[1],"XY") == 0) { act = 2 ;}
  if (narg > 2)  ppfname = arg[2];

  int rc = 0;
  cout << ">>>> treccyl : " << arg[1] << " PPFName=" << ppfname << endl;
  try {
    if (act == 2) rc = testmulticyl(ppfname);
    else rc = test1cyl(ppfname);
  }
  catch (PThrowable& exc) {
    cerr << " treccyl.cc catched Exception " << exc.Msg() << endl;
    rc = 77;
  }  
  catch (std::exception& sex) {
    cerr << "\n treccyl.cc std::exception :" 
         << (string)typeid(sex).name() << "\n msg= " 
         << sex.what() << endl;
  }
  catch (...) {
    cerr << " treccyl.cc catched unknown (...) exception  " << endl; 
    rc = 78; 
  } 

  cout << ">>>> treccyl ------- FIN ----------- Rc=" << rc << endl;
  return rc;
}


//-----------------------------------------------------------------------------
//--- Fonction de test : reconstruction plan AngX-Frequence (1 cylindre)
int test1cyl(string& ppfname)
{

  // BRSourceGen sg;
//  int nsrc = 60; 
  BRSourceGen sg(nsrcmax, maxangX, 0.);
  //  sg.WritePPF(string("brsrc1.ppf"));

  cout << "=== test1cyl: BRSourceGen NbSrc= " << sg.NbSources()
       << " NbRecep=" << MR << " NSamples=" << NE << endl;

  // BRSourceGen  sg(string("brsrc1.ppf"));
  if (prtlevel > 1)  sg.Print(cout);
  
  
  MultiBeamCyl  mb(MR, NE);
  mb.SetPrintLevel(prtlevel);
  mb.SetBaseFreqDa(freq0, da);
  mb.SetNoiseSigma(snoise);
  mb.SetTimeJitter(tjit);
  mb.SetTimeOffsetSigma(tos);
  mb.SetGains(gmean, gsig, nantgz);

  mb.SetSources(sg);
  
  mb.ComputeTimeVectors();
  mb.ComputeSignalVector(0, true);
  cout << "treccy/test1cyl:  signal vectors OK " << endl;
  PrtTim("test1cyl:[1] ");
  
  cout << "--- treccy/test1cyl: Saving to PPF file " << ppfname << endl;

  POutPersist po(ppfname);
//              direct access to variables members !!!!
  po << PPFNameTag("signal") << mb.signal_;
  po << PPFNameTag("sigjitt") << mb.sigjitt_;
  po << PPFNameTag("f_sig") << mb.f_sig_;
  po << PPFNameTag("f_sigjit") << mb.f_sigjit_;

  NTuple ntsrc = sg.Convert2Table(freq0);
  po << PPFNameTag("ntsrc") << ntsrc;

  cout << "treccy/test1cyl: - sig/f_sig,ntsrc to OutPPF OK  " << endl;
  PrtTim("test1cyl[2] ");
  
  mb.ReconstructSourcePlane(true);
  {
  TMatrix<r_4> srcplane = module(mb.getRecSrcPlane() );
  po << PPFNameTag("recsrcplane") << srcplane;
  }
  PrtTim("test1cyl[3] ");

  return 0;

}


//-----------------------------------------------------------------------------
//--- Fonction de test : reconstruction cube AngX-AngY-Frequence (multi-cylindre)
int testmulticyl(string& ppfname)
{

//.............  sources
  // BRSourceGen sg;
  int nsf = 6;
  vector<double> frq;
  frq.push_back(0.1/tClock);
  frq.push_back(0.27/tClock);
  frq.push_back(0.38/tClock);
  
  
  cout << "testmulticyl: BRSourceGen sg([frq=0.1,0.27,0.38], " << nsf 
       << "," << maxangX << "," << maxangY << ")" << endl;
  BRSourceGen sg(frq, nsf, maxangX, maxangY);
  
  int is;
  double fay[6] = {-0.7,-0.5,0.,0.,0.5,0.7};
//  double fay[6] = {-0.2,0.5,-0.3,0.6,-0.1,0.7};
//  double fax[6] = {0.6,-0.2,-0.5,0.4,-0.1,0.3};
  for(is=0; is<3*nsf; is++) {
    int ism = is%nsf;
    sg.angX(is) = maxangX*(ism-2.5)/3.;         //  accessing data member
    sg.angY(is) = maxangY*fay[ism];             //      directly !!!
//    sg.angX(is) = maxangX*fax[ism];       //  accessing data member
  }
  // sg.WritePPF(string("brsrcm.ppf"));
  // BRSourceGen  sg(string("brsrcm.ppf"));
  cout << "=== testmulticyl: NbSrc= " << sg.NbSources() 
       << " NbRecep=" << MR << " NSamples=" << NE << " NCyl=" << nCyl << endl;
  if (prtlevel > 1)  sg.Print(cout);
  
//.......................... cylinders  
  MultiCylinders  mcyl ("telescope.in");
//  MultiCylinders  mcyl (MR, NE);
//  mcyl.SetPrintLevel(prtlevel);
//  mcyl.SetBaseFreqDa(freq0, da);
//  mcyl.SetNoiseSigma(snoise);
//  mcyl.SetTimeJitter(tjit);
//  mcyl.SetTimeOffsetSigma(tos);
//  mcyl.SetGains(gmean, gsig, nantgz);

//  for (int iCyl=0; iCyl<nCyl; iCyl++)
//  {
//     mcyl.AddCylinder(xCyl[iCyl],yCyl[iCyl]);
//  }


  mcyl.SetSources(sg);
  
  PrtTim("testmulticyl[1] ");

  //  mcyl.ReconstructCylinderPlaneS(true);
  mcyl.ReconstructSourceBox(halfNY, maxangY/halfNY, NX, maxangX/NX);

  cout << "--- treccy/testmulticyl: Saving to PPF file " << ppfname << endl;
  POutPersist po(ppfname);

  DVList  dvl;
  dvl("Da") = da;
  po << PPFNameTag("dvl") <<dvl;     
  
  NTuple ntsrc = sg.Convert2Table(freq0);
  po << PPFNameTag("ntsrc") << ntsrc;

    //  TMatrix<r_4> srcplane0 = module(mcyl.GetCylinder(0).getRecSrcPlane());
  TMatrix< complex<r_4> > srcplane0 = mcyl.GetCylinder(0).getRecSrcPlane();
  po << PPFNameTag("recsrcplane0") << srcplane0;
  TMatrix< complex<r_4> > srcplane1 = mcyl.GetCylinder(1).getRecSrcPlane();
  po << PPFNameTag("recsrcplane1") << srcplane1;
//  TMatrix< complex<r_4> > srcplane3 = mcyl.GetCylinder(3).getRecSrcPlane();
//  po << PPFNameTag("recsrcplane3") << srcplane3;
  PrtTim("testmulticyl[2] ");

  po << PPFNameTag("recsrcbox") << mcyl.getRecSrcBox();

//      k= N T frq   with N=2*SizeZ()
  int kfmin = (int)(2.*frq[0]*tClock*(float)mcyl.getRecSrcBox().SizeZ() - 2.);
  int kfmax = kfmin+2;
  cout << "testmulticyl/Info: slice0 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
  TMatrix<r_4> slice0 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf0") << slice0;
  kfmin = (int)(2*frq[1]*tClock*(float)mcyl.getRecSrcBox().SizeZ() - 2.);  
  kfmax = kfmin+2;
  cout << "testmulticyl/Info: slice1 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
  TMatrix<r_4> slice1 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf1") << slice1;
  kfmin = (int)(2*frq[2]*tClock*(float)mcyl.getRecSrcBox().SizeZ() - 2.);  
  kfmax = kfmin+2;
  cout << "testmulticyl/Info: slice2 kfmin=" << kfmin << " kfmax=" << kfmax << endl;
  TMatrix<r_4> slice2 = mcyl.getRecXYSlice(kfmin, kfmax);
  po << PPFNameTag("recXYf2") << slice2;

  PrtTim("testmulticyl[3] ");

  return 0;

}

//---------------------------------------------------------------------
int ReadParam(const char* fileName)
{  
  DataCards dc;
  dc.ReadFile(fileName);
//      frequences are in units of 1/T = 0.5 GHz
//      distance are in units of cT =3E8 * 2E-9=0.60 m
//  double fUnit=0.5;   // 0.5 GHz <=> T = 2 ns
//  double dUnit=0.6;   // distance unit in m.
  double fUnit=1.;      // 0.5 GHz <=> T = 2 ns
  double dUnit=1.;      // distance unit in m.
  
  NE=dc.IParam("nSample");
  freq0=dc.DParam("freq0")/fUnit;
//  tClock=dc.DParam("tClock");
  nCyl=dc.IParam("nCyl");
  for (int i=0; i<nCyl; i++){
    xCyl[i]=dc.DParam("xCyl",i)/dUnit;
    yCyl[i]=dc.DParam("yCyl",i)/dUnit;
  }
  MR=dc.IParam("nAntenna");
  da=dc.DParam("dAntenna")/dUnit;
  maxangX=dc.DParam("angMaxX");
  double cylDiam=dc.DParam("cylinderDiam")/dUnit;
// thetaMax = lambda_M/d = c/freq_min/d;  freq_min = freq0 + 1/2T
  maxangY=cLight/(freq0+1./2./tClock)/cylDiam;
//  cout << "*************** maxangY = " <<maxangY << endl;
//  maxangY=dc.DParam("angMaxY");
  snoise=dc.DParam("noiseSigma");
  tjit=dc.DParam("sigmaTimeJitt");
  tos=dc.DParam("sigmaClockJitt");
  gmean=dc.DParam("meanGain");
  gsig=dc.DParam("sigmaGain");
  nantgz=dc.IParam("nDeadAntenna");
  prtlevel=dc.IParam("printLevel");
  halfNY=dc.IParam("halfNY");
  NX=dc.IParam("NX");
  return 1;
}