Changeset 3192 in Sophya for trunk/Cosmo/RadioBeam/treccyl.cc

Timestamp:

Mar 20, 2007, 11:48:15 AM (19 years ago)

Author:

legoff

Message:

input file read with DVcards by treccyl.cc

and by multiCylinders(char* filename)

using real units (meter and GHz)
using _ for member variables (e.g. x_)

File:

• trunk/Cosmo/RadioBeam/treccyl.cc (modified) (12 diffs)

trunk/Cosmo/RadioBeam/treccyl.cc

@@ -17 +17 @@
 
 #include "timing.h"
+#include "datacards.h"
+#include <dvlist.h>
 
 #include "multicyl.h"
 #include "mbeamcyl.h"
+#define LENGTH 1024
 
 /* 
@@ -32 +35 @@
 // Declaration des fonctions de ce fichier
 static int test1cyl(string& ppfname);
-static int testmulticyl(string& ppfname, int ncyl=5);
+static int testmulticyl(string& ppfname);
+int ReadParam(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 = 1024;  // Nombre d'echantillon en temps;
+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
@@ -44 +53 @@
 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
 
@@ -53 +64 @@
 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];
 //-----------------------------------------------------------
 
@@ -64 +79 @@
 {
 
-SophyaInit();
-InitTim();   // Initializing the CPU timer
-
-string ppfname = "treccyl.ppf";
-int act = 1;
-int ncyl = 5; 
-if (narg < 2) {
-  cout << "Usage: treccyl act ppfname [PrtLev=0] \n"
-       << " -act= X ou XY5 ou XY12  (5 ou 12 cylindres) \n"
+  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],"XY5") == 0) { act = 2;  ncyl = 5; }
-else if (strcmp(arg[1],"XY12") == 0) { act = 2;  ncyl = 12; }
-if (narg > 2)  ppfname = arg[2];
-if (narg > 3)  prtlevel = atoi(arg[3]);
-
-int rc = 0;
-cout << ">>>> treccyl : " << arg[1] << " PPFName=" << ppfname << endl;
-try {
-  if (act == 2) rc = testmulticyl(ppfname, ncyl);
-  else rc = test1cyl(ppfname);
-}
+    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;
@@ -101 +119 @@
   } 
 
-cout << ">>>> treccyl ------- FIN ----------- Rc=" << rc << endl;
-return rc;
+  cout << ">>>> treccyl ------- FIN ----------- Rc=" << rc << endl;
+  return rc;
 }
 
 
+//-----------------------------------------------------------------------------
 //--- Fonction de test : reconstruction plan AngX-Frequence (1 cylindre)
 int test1cyl(string& ppfname)
@@ -111 +130 @@
 
   // BRSourceGen sg;
-  int nsrc = 60; 
+//  int nsrc = 60; 
   BRSourceGen sg(nsrcmax, maxangX, 0.);
   //  sg.WritePPF(string("brsrc1.ppf"));
@@ -140 +159 @@
 
   POutPersist po(ppfname);
-  po << PPFNameTag("signal") << mb.signal;
-  po << PPFNameTag("sigjitt") << mb.sigjitt;
-  po << PPFNameTag("f_sig") << mb.f_sig;
-  po << PPFNameTag("f_sigjit") << mb.f_sigjit;
+//              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);
@@ -163 +183 @@
 
 
+//-----------------------------------------------------------------------------
 //--- Fonction de test : reconstruction cube AngX-AngY-Frequence (multi-cylindre)
-int testmulticyl(string& ppfname, int ncyl)
+int testmulticyl(string& ppfname)
 {
 
-  if ((ncyl != 5) && (ncyl != 12)) ncyl = 5;
-
+//.............  sources
   // BRSourceGen sg;
   int nsf = 6;
   vector<double> frq;
-  frq.push_back(0.1);
-  frq.push_back(0.27);
-  frq.push_back(0.38);
+  frq.push_back(0.1/tClock);
+  frq.push_back(0.27/tClock);
+  frq.push_back(0.38/tClock);
   
   
@@ -183 +203 @@
   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.;
-    sg.angY(is) = maxangY*fay[ism];
+    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;
+       << " NbRecep=" << MR << " NSamples=" << NE << " NCyl=" << nCyl << endl;
   if (prtlevel > 1)  sg.Print(cout);
   
-  
-  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);
-
-  if (ncyl == 12) {
-    cout << " --- testmulticyl/ NCyl= " << ncyl << " posY=0 ... " << (ncyl-1)*5. 
-         << " (EqualSpacing)" << endl;
-    for(int kkc=0; kkc<12; kkc++) {
-      cout << "..." << kkc << " - mcyl.AddCylinder(posY= " << 5.*kkc << " )" << endl;
-      mcyl.AddCylinder(5.*kkc,5.*kkc);
-//      mcyl.AddCylinder(0.,5.*kkc);
-    }
-  }
-  else {
-    cout << " --- testmulticyl/ NCyl= " << ncyl << " posY=0 ... 55 (IrregularSpacing)" << endl;
-    double posyc[5] = {0.,5.,15.,35.,55.};
-    for(int kkc=0; kkc<5; kkc++) {
-      cout << "..." << kkc << " - mcyl.AddCylinder(posY= " << posyc[kkc] << " )" << endl;
-      mcyl.AddCylinder(posyc[kkc],posyc[kkc]);
-    }
-  }
+//.......................... 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);
@@ -226 +238 @@
 
   //  mcyl.ReconstructCylinderPlaneS(true);
-  mcyl.ReconstructSourceBox(10, maxangY/10.);
+  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<r_4> srcplane2 = module(mcyl.GetCylinder(3).getRecSrcPlane());
-  TMatrix< complex<r_4> > srcplane2 = mcyl.GetCylinder(2).getRecSrcPlane();
-  po << PPFNameTag("recsrcplane2") << srcplane2;
-  }
-
-  {
-    //  TMatrix<r_4> srcplane3 = module(mcyl.GetCylinder(3).getRecSrcPlane());
-  TMatrix< complex<r_4> > srcplane3 = mcyl.GetCylinder(0).getRecSrcPlane();
-  po << PPFNameTag("recsrcplane3") << srcplane3;
-  }
-
+  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] ");
 
-  int kfmin, kfmax;
   po << PPFNameTag("recsrcbox") << mcyl.getRecSrcBox();
-  kfmin = mcyl.getRecSrcBox().SizeZ()/0.5*frq[0] - 2;  kfmax = kfmin+2;
+
+//      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 = mcyl.getRecSrcBox().SizeZ()/0.5*frq[1] - 2;  kfmax = kfmin+2;
+  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 = mcyl.getRecSrcBox().SizeZ()/0.5*frq[2] - 2;  kfmax = kfmin+2;
+  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] ");
@@ -280 +283 @@
 
 }
+
+//---------------------------------------------------------------------
+int ReadParam(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;
+}