Changeset 3192 in Sophya

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_)

Location:

trunk/Cosmo/RadioBeam

Files:

• brsrc.cc (modified) (1 diff)
• brsrc.h (modified) (1 diff)
• makefile (modified) (2 diffs)
• mbeamcyl.cc (modified) (6 diffs)
• mbeamcyl.h (modified) (4 diffs)
• multicyl.cc (modified) (8 diffs)
• multicyl.h (modified) (6 diffs)
• treccyl.cc (modified) (12 diffs)

trunk/Cosmo/RadioBeam/brsrc.cc

@@ -71 +71 @@
 
   int j = 0;
-  for(int k=0; k<f.size(); k++) {
+  for(int k=0; k<(int)f.size(); k++) {
     if (f[k] >= 0.5) {
       cout << "BRSourceGen::BRSourceGen() f[k=" << k << "]=" << f[k] << " >=0.5 ignored " << endl;

trunk/Cosmo/RadioBeam/brsrc.h

@@ -49 +49 @@
   NTuple Convert2Table(double freq0=0.);   
 
+//private:
   //-------- Variables / objets membres 
  Vector freq; // frequence des sources

trunk/Cosmo/RadioBeam/makefile

@@ -11 +11 @@
         echo 'makefile : treccyl made'
 
+testFFT : Objs/testFFT.o
 ######
 
@@ -30 +31 @@
         $(CXXCOMPILE) -o Objs/multicyl.o  multicyl.cc
 
+testFFT :  Objs/testFFT.o
+        $(CXXLINK) -o testFFT \
+        Objs/testFFT.o \
+        $(SOPHYAEXTSLBLIST)
+
+Objs/testFFT.o : testFFT.cc
+        $(CXXCOMPILE) -o Objs/testFFT.o testFFT.cc

trunk/Cosmo/RadioBeam/mbeamcyl.cc

@@ -5 +5 @@
 #include "srandgen.h"
 
+static double cLight=0.3;       // in 1E9 m/s
+static double tClock = 2.; // should come from param file !!!!
+//static double cLight=1;       
+//static double tClock = 1.; 
+
 //=================================================
 
 MultiBeamCyl::MultiBeamCyl(int nr, int ns, double posx, double posy)
-  : texact(ns) , tjitt(ns) , toffset(nr) ,
-    signal(ns), sigjitt(ns) , gain(nr)
-{
-  NR = nr;
-  NS = ns;
-  posY = posy;
-  posX = posx;
+  : texact_(ns) , tjitt_(ns) , toffset_(nr) ,
+    gain_(nr), signal_(ns), sigjitt_(ns) 
+{
+  NR_ = nr;
+  NS_ = ns;
+  posY_ = posy;
+  posX_ = posx;
 
   SetPrintLevel(0);
@@ -23 +28 @@
   SetTimeOffsetSigma(0.);
   SetGains(1., 0., 0);
-  adfg = false; src = NULL;
+  adfg_ = false; src_ = NULL;
   SetSources(new BRSourceGen, true);
 }
 
+//-----------------------------------------------------------------------------
 MultiBeamCyl::~MultiBeamCyl()
 {
-  if (adfg && src) delete src;
-}
-
+  if (adfg_ && src_) delete src_;
+}
+
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::SetSources(BRSourceGen* brs, bool ad)
 {
   if (brs == NULL) return;
-  if (adfg && src) delete src;
-  src = brs;  adfg=ad;
-  if (PrtLev > 1) 
+  if (adfg_ && src_) delete src_;
+  src_ = brs;  adfg_=ad;
+  if (PrtLev_ > 1) 
     cout << " MultiBeamCyl::SetSources(brs=" << brs <<" ,bool) NbSrc=" 
-         << src->NbSources() << endl;
-
-}
-
+         << src_->NbSources() << endl;
+
+}
+
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::SetGains(double g, double sigg, int nzerogain)
 {
-  if (sigg < 1.e-6)  gain = g;
-  else gain = RandomSequence(RandomSequence::Gaussian, g, sigg);
+  if (sigg < 1.e-6)  gain_ = g;
+  else gain_ = RandomSequence(RandomSequence::Gaussian, g, sigg);
   int k;
-  for (k=0; k<NR; k++) if (gain(k) < 0) gain(k) = 0.;
+  for (k=0; k<NR_; k++) if (gain_(k) < 0) gain_(k) = 0.;
   for(k=0; k<nzerogain; k++) {
-    int zg = random()%NR;
-    if ((zg >=0) && (zg < NR))  gain(zg) = 0.; 
-  }
-  if (PrtLev > 1) 
+    int zg = random()%NR_;
+    if ((zg >=0) && (zg < NR_))  gain_(zg) = 0.; 
+  }
+  if (PrtLev_ > 1) 
     cout << " MultiBeamCyl::SetGains(g=" << g <<" ,sigg=" << sigg 
          << " ,nzg=" << nzerogain << " )" << endl;
 }
 
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::ComputeTimeVectors()
 {
-  texact = RegularSequence(0., 1.);
-  toffset = RandomSequence(RandomSequence::Gaussian, 0., toffsig);
+  texact_ = RegularSequence(0., tClock);        // 0, tClock, 2*tClock, ...
+//  for (int i=0;i<1024;i++) {cout << texact(i)<< endl;};
+  toffset_ = RandomSequence(RandomSequence::Gaussian, 0., toffsig_);
   NewTJitVector();
 }
 
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::NewTJitVector(int num)
 {
-  if (timejitter > 1.e-19) {
-    tjitt = RandomSequence(RandomSequence::Gaussian, 0., timejitter);
-    tjitt += texact;  
-  }
-  else tjitt = texact;  
-  if (num >= 0) tjitt += toffset(num);
+  if (timejitter_ > 1.e-19) {
+    tjitt_ = RandomSequence(RandomSequence::Gaussian, 0., timejitter_);
+    tjitt_ += texact_;  
+  }
+  else tjitt_ = texact_;  
+  if (num >= 0) tjitt_ += toffset_(num);
   
 }
 
+//-----------------------------------------------------------------------------
 int MultiBeamCyl::ComputeSignalVector(int num, bool fgsignojit)
 {
   int nok = 0;
-  signal = 0.;
-  sigjitt = 0.;
-  for(int is=0; is<src->freq.Size(); is++) {
-    double fr = src->freq(is);
+  signal_ = 0.;
+  sigjitt_ = 0.;
+  for(int is=0; is<src_->freq.Size(); is++) {
+    double fr = src_->freq(is);
     if ((fr < 0.) || (fr > 0.5)) continue;
     nok++;
@@ -88 +100 @@
     // pas celle apres shift (freq-reduite)
     // lambda = c T = c/freq avec c = 1, dephasage = 2*pi*num*Da*sin(ang)/lambda
-    double dephasage = (posX+num*Da)*sin(src->angX(is)) + 
-                                posY*sin(src->angY(is)) ;
-    dephasage *= (2*M_PI*(fr+freq0));
+    double dephasage = (posX_+num*Da_)*sin(src_->angX(is)) 
+                                + posY_*sin(src_->angY(is)) ;
+    dephasage *= 2*M_PI*(fr+freq0_)/cLight;
     // On ajoute alors la phase propre de chaque source
-    dephasage += src->phase(is);
-    double amprep = src->amp(is)*AngResponse(src->angX(is), src->angY(is));
-    for(int k=0; k<NS; k++) {
-      sigjitt(k) += amprep*sin(2.*M_PI*fr*tjitt(k)+dephasage);
+    dephasage += src_->phase(is);
+    double amprep = src_->amp(is)*AngResponse(src_->angX(is), src_->angY(is));
+    for(int k=0; k<NS_; k++) {
+      sigjitt_(k) += amprep*sin(2.*M_PI*fr*tjitt_(k)+dephasage);
       if (fgsignojit) 
-        signal(k) += amprep*sin(2.*M_PI*fr*texact(k)+dephasage);
+        signal_(k) += amprep*sin(2.*M_PI*fr*texact_(k)+dephasage);
     }
   }
   // Application du gain du detecteur 
-  r_4 ga = gain(num);
+  r_4 ga = gain_(num);
   if (fabs(ga-1.) > 1.e-9) {
-    signal *= ga;
-    sigjitt *= ga;
+    signal_ *= ga;
+    sigjitt_ *= ga;
   }
   // Ajout de bruit (ampli ...) 
-  if (signoise > 1.e-18) {
-    for(int k=0; k<NS; k++) {
-      sigjitt(k) += GauRnd(0., signoise);
-      if (fgsignojit) signal(k) += GauRnd(0., signoise);
-    }
-  }
-
+  if (signoise_ > 1.e-18) {
+    for(int k=0; k<NS_; k++) {
+      sigjitt_(k) += GauRnd(0., signoise_);
+      if (fgsignojit) signal_(k) += GauRnd(0., signoise_);
+    }
+  }
+//............      FFT in time
   FFTPackServer ffts;
-
-  ffts.FFTForward(sigjitt, f_sigjit);
-  if (fgsignojit) ffts.FFTForward(signal, f_sig);
+  ffts.FFTForward(sigjitt_, f_sigjit_);
+  if (fgsignojit) ffts.FFTForward(signal_, f_sig_);
 
   return nok;
@@ -123 +134 @@
 
 
+//-----------------------------------------------------------------------------
 /*  --- a supprimer ?
 inline float myZmodule(complex<r_4>& z) 
@@ -130 +142 @@
 ----- */
 
+//-----------------------------------------------------------------------------
 void MultiBeamCyl::ReconstructSourcePlane(bool fgzerocentre)
 {
@@ -135 +148 @@
   int noksrc = ComputeSignalVector(0, false);
   vector<TVector< complex<r_4> > >  vvfc;
-  cout << "MultiBeamCyl::ReconstructSourcePlane() NR=" << NR 
-       << " PosY=" << posY << " NFreq=" << f_sigjit.Size()-1 << " NOkSrc=" << noksrc << endl;
-  if (PrtLev > 0) {
-    cout << " ... posY= " << posY << " MeanGain=" << Mean(gain) << " MeanAmpSrc=" 
-         << Mean(src->amp) << endl;
-    cout << " ...SigNoise=" << signoise << " TimeJitter=" << timejitter 
-         << " TOffSig=" << toffsig << " Da=" << Da << " Freq0=" << freq0 << endl;
+  cout << "MultiBeamCyl::ReconstructSourcePlane() NR=" << NR_ 
+       << " PosY=" << posY_ << " NFreq=" << f_sigjit_.Size()-1 << " NOkSrc=" << noksrc << endl;
+  if (PrtLev_ > 0) {
+    cout << " ... posY= " << posY_ << " MeanGain=" << Mean(gain_) << " MeanAmpSrc=" 
+         << Mean(src_->amp) << endl;
+    cout << " ...SigNoise=" << signoise_ << " TimeJitter=" << timejitter_ 
+         << " TOffSig=" << toffsig_ << " Da=" << Da_ << " Freq0=" << freq0_ << endl;
   }
   // On ne s'occupe pas de la composante continue
-  for(int jf=1; jf<f_sigjit.Size(); jf++) {
-    TVector<complex<r_4> > cf(NR);
+  for(int jf=1; jf<f_sigjit_.Size(); jf++) {
+    TVector<complex<r_4> > cf(NR_);
     vvfc.push_back(cf);
   }
   cout << "ReconstructSourcePlane()/Info: computing s(t) for each receptor ..." << endl;
-  int pmod = NR/10;
-  for(int ir=0; ir<NR; ir++) {
-    if (timejitter > 1.e-19) NewTJitVector(ir);
+  int pmod = NR_/10;
+  for(int ir=0; ir<NR_; ir++) {
+    if (timejitter_ > 1.e-19) NewTJitVector(ir);
     noksrc = ComputeSignalVector(ir, false);
-    for(int jf=1; jf<f_sigjit.Size(); jf++)
-      vvfc[jf-1](ir) = f_sigjit(jf);
-    if ( (PrtLev>0) && (ir%pmod == 0) )
-      cout << " OK s(t) for ir=" << ir << " / NR=" << NR << " NOkSrc=" << noksrc << endl;
+    for(int jf=1; jf<f_sigjit_.Size(); jf++)
+      vvfc[jf-1](ir) = f_sigjit_(jf);
+    if ( (PrtLev_>0) && (ir%pmod == 0) )
+      cout << " OK s(t) for ir=" << ir << " / NR=" << NR_ << " NOkSrc=" << noksrc << endl;
   }
 
   cout << "ReconstructSourcePlane()/Info: computing s(ang) for each freq by FFT" << endl;
-  cmplx_srcplane.SetSize(f_sigjit.Size()-1, NR);
-  // rec_srcplane.SetSize(f_sigjit.Size()-1, NR);
+  cmplx_srcplane_.SetSize(f_sigjit_.Size()-1, NR_);
+  // rec_srcplane.SetSize(f_sigjit_.Size()-1, NR_);
   FFTPackServer ffts;
   TVector<complex<r_4> > fcf;
   pmod = vvfc.size()/10;
-  for(int jf=0; jf<vvfc.size(); jf++) {
-    ffts.FFTForward(vvfc[jf], fcf);
-    if (fcf.Size() != NR) {
+  for(int jf=0; jf<(int)vvfc.size(); jf++) {            // loop over frequencies
+    ffts.FFTForward(vvfc[jf], fcf);     // FFT alomg cylinder
+    if (fcf.Size() != NR_) {
       cout << "ReconstructSourcePlane()/BUG jf=" << jf << " fcf.Size() != NR "
-           << fcf.Size() << " != " << NR << endl;
+           << fcf.Size() << " != " << NR_ << endl;
       continue;
     }
     if (fgzerocentre) {  // On veut avoir la direction angle=0 au milieu de la matrice
-      int milieu = (NR-1)/2;
-      if (NR%2 == 0) milieu++;
-      int decal = NR - milieu - 1;
+      int milieu = (NR_-1)/2;
+      if (NR_%2 == 0) milieu++;
+      int decal = NR_ - milieu - 1;
       for (int ir=0; ir<=milieu; ir++) 
-        cmplx_srcplane(jf, ir+decal) = fcf(ir);
+        cmplx_srcplane_(jf, ir+decal) = fcf(ir);
       //        rec_srcplane(jf, ir+decal) = myZmodule(fcf(ir));
 
-      for (int ir=milieu+1; ir<NR; ir++) 
-        cmplx_srcplane(jf, decal-(NR-ir)) = fcf(ir);
-      //        rec_srcplane(jf, decal-(NR-ir)) = myZmodule(fcf(ir));
+      for (int ir=milieu+1; ir<NR_; ir++) 
+        cmplx_srcplane_(jf, decal-(NR_-ir)) = fcf(ir);
+      //        rec_srcplane_(jf, decal-(NR_-ir)) = myZmodule(fcf(ir));
       
     }
-    else for (int ir=0; ir<NR; ir++) 
-      cmplx_srcplane(jf, ir) = fcf(ir);
+    else for (int ir=0; ir<NR_; ir++) 
+      cmplx_srcplane_(jf, ir) = fcf(ir);
     //      rec_srcplane(jf, ir) = myZmodule(fcf(ir));
 
-    if ( (PrtLev > 0) && (jf%pmod == 0)) 
+    if ( (PrtLev_ > 0) && (jf%pmod == 0)) 
       cout << " OK rec_srcplane(jf, ir) for jf=" << jf << endl;
   }

trunk/Cosmo/RadioBeam/mbeamcyl.h

@@ -31 +31 @@
   // ns = nb d'echantillons en temps de chaque paquet
   // posy = position spatiale en Y pour chaque cylindre (unite ou lambda=1 pour f=1)
- MultiBeamCyl(int nr=1024, int ns=4096, double posx=0., double posy=0.);
+  // posx = position spatiale en X pour chaque cylindre (unite ou lambda=1 pour f=1)
+ MultiBeamCyl(int nr=1024, int ns=4096, double posx=0., double posy=0. );
  ~MultiBeamCyl();
 
  // Niveau de print de debug
- inline void SetPrintLevel(int prl=0) { PrtLev=prl; }
+ inline void SetPrintLevel(int prl=0) { PrtLev_=prl; }
 
  // Reponse angulaire pour un recepteur (cylindre+antenne)
@@ -41 +42 @@
  inline double AngResponse(double angx, double angy) { return 1.; }
 
- inline double getCylinderYPos() { return posY; }
- inline double getCylinderXPos() { return posX; }
+ inline double getCylinderYPos() { return posY_; }
+ inline double getCylinderXPos() { return posX_; }
+
  // Specification de la frequence de base f0 et espacement des recepteurs
  // freq_vrai = freq_BRSourceGen + f0
  // pour f_vrai = 2, lambda = 0.5 et lambda/2 = 0.25
  inline void SetBaseFreqDa(double f0=2., double a=0.25) 
-   {  freq0 = f0;   Da = a; }
+   {  freq0_ = f0;   Da_ = a; }
 
  // frequences reduites (entre 0 ... 0.5) , angle en radian, amp ~ 1
@@ -55 +57 @@
 
  // Definition du sigma du bruit gaussien sur les echantillons
- inline void SetNoiseSigma(double sig=0.) {  signoise = sig; }
+ inline void SetNoiseSigma(double sig=0.) {  signoise_ = sig; }
  // Definition du sigma du jitter d'horloge (typ 0.01)
- inline void SetTimeJitter(double tjit=0.) { timejitter = tjit; } 
+ inline void SetTimeJitter(double tjit=0.) { timejitter_ = tjit; } 
  // Definition du sigma des offsets d'horloge entre recepteurs (typ 0.02)
- inline void SetTimeOffsetSigma(double tsig=0.) { toffsig = tsig; }
+ inline void SetTimeOffsetSigma(double tsig=0.) { toffsig_ = tsig; }
 
  // Definition du gain et sigma de fluctuations de gain 
@@ -79 +81 @@
  void ReconstructSourcePlane(bool fgzerocentre=true);
 
- inline  TMatrix< complex<r_4> > & getRecSrcPlane() { return cmplx_srcplane; }
+ inline  TMatrix< complex<r_4> > & getRecSrcPlane() { return cmplx_srcplane_; }
 
  //-------------- Variables membres ------------ 
- int NR, NS;  // nb recepteurs, nb d'echantillons ds chaque paquet
- double posY;    // Position selon Y (E-O) de la ligne de recepteurs
- double posX;    // Position selon X (N-S) de la ligne de recepteurs
+ int NR_, NS_;  // nb recepteurs, nb d'echantillons ds chaque paquet
+ double posX_;    // Position selon X (N-S) du premier recepteur du cylindre
+ double posY_;    // Position selon Y (E-O) de la ligne de recepteurs
 
- int PrtLev;   // Niveau de print de debug
+ int PrtLev_;   // Niveau de print de debug
 
- double Da; // distance entre recepteurs
- double freq0;  // frequence de base (freqvrai=f+freq0)
+ double Da_; // distance entre recepteurs
+ double freq0_;  // frequence de base (freqvrai=f+freq0)
  
- Vector texact;  // les temps exact
- Vector tjitt;   // temps avec jitter
- Vector toffset;  // Offset en temps entre recepteurs 
- double timejitter; // le sigma du jitter en temps des coups d'horloge
- double toffsig;    // sigma des offsets en temps
+ Vector texact_;  // les temps exact
+ Vector tjitt_;   // temps avec jitter
+ Vector toffset_;  // Offset en temps entre recepteurs 
+ double timejitter_; // le sigma du jitter en temps des coups d'horloge
+ double toffsig_;    // sigma des offsets en temps
  
- BRSourceGen* src;  // Les sources
- bool adfg; // if true, delete src 
+ BRSourceGen* src_;  // Les sources
+ bool adfg_; // if true, delete src 
 
- double signoise; // sigma du bruit additif (bruit ampli ... )
+ double signoise_; // sigma du bruit additif (bruit ampli ... )
  
- TVector<r_4> gain;  // Gain de chaque recepteur 
+ TVector<r_4> gain_;  // Gain de chaque recepteur 
 
- TVector<r_4> signal; // signal sans jitter en temps 
- TVector<r_4> sigjitt; // signal avec jitter en temps
+ TVector<r_4> signal_; // signal sans jitter en temps 
+ TVector<r_4> sigjitt_; // signal avec jitter en temps
 
- TVector< complex<r_4> > f_sig, f_sigjit;  // TF des vecteurs signal , sigjitt
+ TVector< complex<r_4> > f_sig_, f_sigjit_;  // TF des vecteurs signal , sigjitt
 
- TMatrix< complex<r_4> > cmplx_srcplane; // composantes complexe du plan-source reconstruit
- // TMatrix<r_4> rec_srcplane;  // Matrice plan source colonnes->angX, lignes->freq 
+ TMatrix< complex<r_4> > cmplx_srcplane_; // composantes complexe du plan-source reconstruit
+ // TMatrix<r_4> rec_srcplane_;  // Matrice plan source colonnes->angX, lignes->freq 
  };
 

trunk/Cosmo/RadioBeam/multicyl.cc

@@ -6 +6 @@
 #include "ctimer.h"
 #include "resusage.h"
+#include "datacards.h"
+
+static double cLight=0.3;       // in 1E9 m/s
+static double tClock = 2.; // should come from param file !!!!
+//static double cLight=1;       
+//static double tClock = 1.; 
 
 
@@ -12 +18 @@
 MultiCylinders::MultiCylinders(int nr, int ns)
 {
-  NR = nr;
-  NS = ns;
+  NR_ = nr;
+  NS_ = ns;
 
   SetPrintLevel(0);
@@ -22 +28 @@
   SetTimeOffsetSigma(0.);
   SetGains(1., 0., 0);
-  adfg = false; src = NULL;
+  adfg_ = false; src_ = NULL;
   SetSources(new BRSourceGen, true);
 }
 
+//-----------------------------------------------------------------------------
+MultiCylinders::MultiCylinders(char* fileName)
+{
+  adfg_ = false; src_ = NULL;
+  SetSources(new BRSourceGen, true);
+  
+//      read telescope parameters and fill variable members
+  DataCards dc;
+  dc.ReadFile(fileName);
+//      in old versions frequences were in units of 1/T = 0.5 GHz
+//      and distances 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.
+//      now f and d in real units 
+  double fUnit=1.;      // 0.5 GHz <=> T = 2 ns
+  double dUnit=1.;      // distance unit in m.
+  
+  NR_=dc.IParam("nAntenna");
+  NS_=dc.IParam("nSample");
+  PrtLev_=dc.IParam("printLevel");
+  Da_=dc.DParam("dAntenna")/dUnit;
+  freq0_=dc.DParam("freq0")/fUnit;
+  timejitter_=dc.DParam("sigmaTimeJitt");
+  toffsig_=dc.DParam("sigmaClockJitt");
+  signoise_=dc.DParam("noiseSigma");
+  gain_=dc.DParam("meanGain");
+  siggain_=dc.DParam("sigmaGain");
+  ngainzero_=dc.IParam("nDeadAntenna");
+  
+//  tClock=dc.DParam("tClock");
+  int nCyl=dc.IParam("nCyl");
+  for (int i=0; i<nCyl; i++){
+    double xCyl=dc.DParam("xCyl",i)/dUnit;
+    double yCyl=dc.DParam("yCyl",i)/dUnit;
+    AddCylinder(xCyl,yCyl);
+  }
+//  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;
+//  halfNY=dc.IParam("halfNY");
+//  NX=dc.IParam("NX");
+}
+
+//-----------------------------------------------------------------------------
 MultiCylinders::~MultiCylinders()
 {
-  if (adfg && src) delete src;
-  for(int k=0; k<mCyl.size(); k++) delete mCyl[k];
-}
-
+  if (adfg_ && src_) delete src_;
+  for(int k=0; k<(int)mCyl_.size(); k++) delete mCyl_[k];
+}
+
+//-----------------------------------------------------------------------------
 MultiBeamCyl& MultiCylinders::GetCylinder(int k)
 {
-  if ((k < 0) || (k >= mCyl.size()))
+  if ((k < 0) || (k >= (int)mCyl_.size())) {
+    cout <<"******************************************* k="<<k<<endl;
     throw RangeCheckError("MultiCylinders::GetCylinder(k) k<0 OR k>=NCyl");
-  return (*mCyl[k]);
-}
-
+  }
+  return (*mCyl_[k]);
+}
+
+//-----------------------------------------------------------------------------
 void MultiCylinders::SetSources(BRSourceGen* brs, bool ad)
 {
   if (brs == NULL) return;
-  if (adfg && src) delete src;
-  src = brs;  adfg=ad;
-}
-
-
+  if (adfg_ && src_) delete src_;
+  src_ = brs;  adfg_=ad;
+}
+
+
+//-----------------------------------------------------------------------------
 void MultiCylinders::ConfigureCylinders()
 {
-  cout << " MultiCylinders::ConfigureCylinders() with NCyl= " << mCyl.size() << endl;
-  for(int k=0; k<mCyl.size(); k++) {
-    mCyl[k]->SetPrintLevel(PrtLev);
-    mCyl[k]->SetBaseFreqDa(freq0, Da);
-    mCyl[k]->SetNoiseSigma(signoise);
-    mCyl[k]->SetTimeJitter(timejitter);
-    mCyl[k]->SetTimeOffsetSigma(toffsig);
-    mCyl[k]->SetGains(gain, siggain, ngainzero);
-    mCyl[k]->SetSources(src, false);
-  }
-}
-
-
+  cout << " MultiCylinders::ConfigureCylinders() with NCyl= " << mCyl_.size() << endl;
+  for(int k=0; k<(int)mCyl_.size(); k++) {
+    mCyl_[k]->SetPrintLevel(PrtLev_);
+    mCyl_[k]->SetBaseFreqDa(freq0_, Da_);
+    mCyl_[k]->SetNoiseSigma(signoise_);
+    mCyl_[k]->SetTimeJitter(timejitter_);
+    mCyl_[k]->SetTimeOffsetSigma(toffsig_);
+    mCyl_[k]->SetGains(gain_, siggain_, ngainzero_);
+    mCyl_[k]->SetSources(src_, false);
+  }
+}
+
+
+//-----------------------------------------------------------------------------
 void MultiCylinders::ReconstructCylinderPlaneS(bool fgzerocentre)
 {
   Timer tm("RecCylPlaneS");
   ResourceUsage resu;
-  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - NCyl= " << mCyl.size() 
+  cout << " MultiCylinders::ReconstructCylinderPlaneS()/Info - NCyl= " << mCyl_.size() 
        << " MemSize=" << resu.getMemorySize() << " kb" << endl;       
   ConfigureCylinders();
 
   char buff[128];
-  for(int k=0; k<mCyl.size(); k++) {
+  for(int k=0; k<(int)mCyl_.size(); k++) {
     cout << "---- Cyl[" << k << "] Calling MultiBeamCyl.ReconstructSourcePlane() ..." << endl;
-    mCyl[k]->ReconstructSourcePlane(fgzerocentre);
+    mCyl_[k]->ReconstructSourcePlane(fgzerocentre);
     sprintf(buff,"Cyl[%d].RecSrcPlane()",k);
     tm.Split(buff);
@@ -82 +139 @@
 }
 
-void MultiCylinders::ReconstructSourceBox(int halfny, double stepangy)
-{
+//-----------------------------------------------------------------------------
+void MultiCylinders::ReconstructSourceBox(int halfny, double stepangy, 
+        int nx, double stepangx)
+{
+  nx=256;
   ReconstructCylinderPlaneS(true);
   TMatrix< complex<r_4> > & mtx = GetCylinder(0).getRecSrcPlane(); 
-  // boite 3D X:angX, Y:angY , Z: freq
+// boite 3D X:angX, Y:angY , Z: freq
+//      if all cylinders at same x position NX is set to zero
+//      => x-size = mtx.NCols() = numbers of receptors per cylinders
+//              move that to readParam() ?
   sa_size_t sz[5] = {0,0,0,0,0};
-//  int nx=mtx.NCols(); // uncomment to go back to nxbin = nantenna
-  int nx=256;
   sz[0] = nx;
   sz[1] = halfny*2+1;
   sz[2] = mtx.NRows();
   TArray< complex<r_4> > & box = getRecSrcBox();
-  box.ReSize(3, sz);
-  //  box = complex< r_4 >( 0.f, 0.f );
+  box.ReSize(3, sz);            //  values initialized to zero ?
   cout << " MultiCylinders::ReconstructSourceBox(" << halfny << "," << stepangy
        << "=" << Angle(stepangy).ToArcMin() << " ) srcbox:" << endl;
@@ -102 +162 @@
 
 
-  int pmod = mtx.NRows()/10;
-  
-  double fstep = 1.0/(double)NS;  // pas en frequence, attention, on a vire la composante continu
+//  int pmod = mtx.NRows()/10;
+  
+  double fstep = 1.0/(double)NS_/tClock;  // pas en frequence, 
+                                // attention, on a vire la composante continu
+//  bool first=true;
   for(int kf=0; kf<mtx.NRows(); kf++) { // Loop over frequencies
-    double frq = (double)(kf+1.)*fstep + freq0;  // + 1 car f=0 a ete vire
-    
-    for(int lc=0; lc<mCyl.size(); lc++) {  // Loop over cylinders 
+    double frq = (double)(kf+1.)*fstep + freq0_;  // + 1 car f=0 a ete vire
+//            then up to frq =   (mtx.NRows() +1 ) * fstep            !!?
+//    cout<<"************"<<mCyl.size()
+    for(int lc=0; lc<(int)mCyl_.size(); lc++) {  // Loop over cylinders 
       MultiBeamCyl& cyl = GetCylinder(lc);
 
       TMatrix< complex<r_4> > & recp = cyl.getRecSrcPlane(); 
 
-      double facl = - 2*M_PI*frq*cyl.getCylinderYPos();  // attention au signe -
-//      double facl_x = - 2*M_PI*cyl.getCylinderXPos()/cyl.Da/(double)cyl.NR; 
-      double facl_x = - 2*M_PI*cyl.getCylinderXPos()/cyl.Da; 
-      double dphi;
+      double facl_y = - 2*M_PI*frq*cyl.getCylinderYPos()/cLight;  // attention signe -
+      double facl_x = - 2*M_PI*cyl.getCylinderXPos()/cyl.Da_; 
       complex< r_4 > phasefactor;
       int jyy = 0;
+      
       for(int jy=-halfny; jy<=halfny; jy++) {  // Loop over Y angular steps
-        for(int ix=0; ix<nx; ix++) {  // Loop over AngX directions
-          double dphi = facl * sin( (double)jy*stepangy )
-                        + facl_x*(double(ix)/double(nx)-1./2.);
-          phasefactor = complex< r_4 > ((r_4) cos(dphi) , (r_4) sin(dphi));
-          // sur recp : index ligne -> frequence  ,    index colonne -> angX ,
-          int ixx=(int)(ix*(double)cyl.NR/double(nx));
-          box(ix, jyy, kf) += recp(kf, ixx)*phasefactor;
-        }  //  
-        jyy++;
+      for(int ix=0; ix<nx; ix++) {  // Loop over AngX directions
+        double dphi = facl_y * sin( (double)jy*stepangy )
+                      + facl_x*(double(ix)/double(nx)-1./2.);
+        phasefactor = complex< r_4 > ((r_4) cos(dphi) , (r_4) sin(dphi));
+        // sur recp : index ligne -> frequence  ,    index colonne -> angX ,
+        int ixx=(int)(ix*(double)cyl.NR_/double(nx));
+        box(ix, jyy, kf) += recp(kf, ixx)*phasefactor;
+      }  //  
+      jyy++;
       } // End of Loop over Y angles
     } // End of loop over cylinders
+
     tm.Nop();
-    if ( (PrtLev>0) && (kf%pmod == 0) )   
+//    if ( (PrtLev_>0) && (kf%pmod == 0) )   
+    if ( (PrtLev_>0) && (kf%(mtx.NRows()/10) == 0) )   
       cout << " OK box(angx,angy, freq=kf) done for kf=" << kf 
-           << " / Max_kf=" << mtx.NRows() << endl;
+         << " / Max_kf=" << mtx.NRows() << endl;
   } // End of loop over over frequencies
 
@@ -140 +204 @@
 } 
 
+//-----------------------------------------------------------------------------
 inline float myZmodule(complex<r_4>& z) 
 {
@@ -145 +210 @@
 }
 
+//-----------------------------------------------------------------------------
 TMatrix< r_4 >  MultiCylinders::getRecXYSlice(int kfmin, int kfmax)
 {
@@ -158 +224 @@
   return(rmtx);
 }
+//-----------------------------------------------------------------------------
+TMatrix< r_4 >  MultiCylinders::getRecYXSlice(int kfmin, int kfmax)
+{
+  TArray< complex<r_4> > & box = getRecSrcBox();  
+  if ((kfmin < 0) || (kfmin >=  box.SizeZ()) || (kfmax < kfmin) || (kfmax >=  box.SizeZ()) )
+    throw RangeCheckError("MultiCylinders::getRecXYSlice(kfmin, kfmax)");
+  TMatrix< r_4> rmtx(box.SizeX(), box.SizeY());
+  for(int kf=kfmin; kf<=kfmax; kf++) {
+    for(int jy=0; jy<box.SizeY(); jy++) 
+      for(int ix=0; ix<box.SizeX(); ix++) 
+        rmtx(ix, jy) += myZmodule(box(ix, jy, kf));
+  }
+  return(rmtx);
+}

trunk/Cosmo/RadioBeam/multicyl.h

@@ -28 +28 @@
   //   longueur @ f=2 ~ 64 (256*lambda/2 = 256*0.25)
  MultiCylinders(int nr=256, int ns=1024);
+ MultiCylinders(char* filename);
  ~MultiCylinders();
 
  // Niveau de print de debug
- inline void SetPrintLevel(int prl=0) { PrtLev=prl; }
+ inline void SetPrintLevel(int prl=0) { PrtLev_=prl; }
 
- // Ajout d'un cylindre, en position posY
+ // Ajout d'un cylindre, en position posX, posY
  inline int AddCylinder(double posx, double posy)
-   {  mCyl.push_back( new MultiBeamCyl(NR, NS, posx, posy) ); return mCyl.size(); }
+   {  mCyl_.push_back( new MultiBeamCyl(NR_, NS_, posx, posy) ); return mCyl_.size(); }
 
  MultiBeamCyl & GetCylinder(int i);
@@ -41 +42 @@
  // Specification de la frequence de base f0 et espacement des recepteurs
  inline void SetBaseFreqDa(double f0=2., double a=0.25) 
-   {  freq0 = f0;   Da = a;  }
+   {  freq0_ = f0;   Da_ = a;  }
 
  // frequences reduites (entre 0 ... 0.5) , angle en radian, amp ~ 1
@@ -49 +50 @@
 
  // Definition du sigma du bruit gaussien sur les echantillons
- inline void SetNoiseSigma(double sig=0.) {  signoise = sig; }
+ inline void SetNoiseSigma(double sig=0.) {  signoise_ = sig; }
  // Definition du sigma du jitter d'horloge (typ 0.01)
- inline void SetTimeJitter(double tjit=0.) { timejitter = tjit; } 
+ inline void SetTimeJitter(double tjit=0.) { timejitter_ = tjit; } 
  // Definition du sigma des offsets d'horloge entre recepteurs (typ 0.02)
- inline void SetTimeOffsetSigma(double tsig=0.) { toffsig = tsig; }
+ inline void SetTimeOffsetSigma(double tsig=0.) { toffsig_ = tsig; }
 
  // Definition du gain et sigma de fluctuations de gain 
  // nzerogain: nb de recepteurs (choisis au hasard) avec gain mis a zero
  inline void SetGains(double g=1., double sigg=0., int nzerogain=0)
-   {  gain=g;  siggain=sigg;  ngainzero = nzerogain; }
+   {  gain_=g;  siggain_=sigg;  ngainzero_ = nzerogain; }
 
 
@@ -72 +73 @@
  // @f = 2 , lambda = 0.5 ===> posY <~ lambda/(2 sin resol) 
  //   ===> posY < ~ 20   
- void ReconstructSourceBox(int halfny=10, double stepangy=M_PI/300);
+// void ReconstructSourceBox(int halfny=10, double stepangy=M_PI/300);
+ void ReconstructSourceBox(int halfny, double stepangy, int ny, double stepangx);
 
 
@@ -78 +80 @@
  // avec la moyenne des modules entre kfmin <= k (selon z) <= kfmax 
  TMatrix< r_4 > getRecXYSlice(int kfmin, int kfmax);
+ TMatrix< r_4 > getRecYXSlice(int kfmin, int kfmax);
 
  // Acces a la boite 3D de sources reconstruite
- inline TArray< complex<r_4> > & getRecSrcBox() { return cmplx_srcbox; }
+ inline TArray< complex<r_4> > & getRecSrcBox() { return cmplx_srcbox_; }
 
  // Configure chaque cylindre a partir des parametres de la classe 
@@ -86 +89 @@
  void ConfigureCylinders();
 
+private:
  //-------------- Variables membres ------------ 
- int NR, NS;  // nb recepteurs, nb d'echantillons ds chaque paquet
+ int NR_, NS_;  // nb recepteurs, nb d'echantillons ds chaque paquet
 
- int PrtLev;   // Niveau de print de debug
+ int PrtLev_;   // Niveau de print de debug
 
- double Da; // distance entre recepteurs
- double freq0;  // frequence de base (freqvrai=f+freq0)
+ double Da_; // distance entre recepteurs
+ double freq0_;  // frequence de base (freqvrai=f+freq0)
  
- double timejitter; // le sigma du jitter en temps des coups d'horloge
- double toffsig;    // sigma des offsets en temps
+ double timejitter_; // le sigma du jitter en temps des coups d'horloge
+ double toffsig_;    // sigma des offsets en temps
  
- BRSourceGen* src;  // Les sources
- bool adfg; // if true, delete src 
+ BRSourceGen* src_;  // Les sources
+ bool adfg_; // if true, delete src 
 
- double signoise; // sigma du bruit additif (bruit ampli ... )
- double gain, siggain;
- int ngainzero;
+ double signoise_; // sigma du bruit additif (bruit ampli ... )
+ double gain_, siggain_;
+ int ngainzero_;
 
- vector<MultiBeamCyl *> mCyl;  // Les differents cylindres 
+ vector<MultiBeamCyl *> mCyl_;  // Les differents cylindres 
 
- TArray< complex<r_4> > cmplx_srcbox; // boite3D des sources (angX,Y,freq) reconstruit
+ TArray< complex<r_4> > cmplx_srcbox_; // boite3D des sources (angX,Y,freq) reconstruit
  };
 

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;
+}