#include "sopnamsp.h" #include "machdefs.h" #include #include #include #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 #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="< 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 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 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 srcplane0 = module(mcyl.GetCylinder(0).getRecSrcPlane()); TMatrix< complex > srcplane0 = mcyl.GetCylinder(0).getRecSrcPlane(); po << PPFNameTag("recsrcplane0") << srcplane0; TMatrix< complex > srcplane1 = mcyl.GetCylinder(1).getRecSrcPlane(); po << PPFNameTag("recsrcplane1") << srcplane1; // TMatrix< complex > 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 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 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 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