Merge pull request #56 from project8/develop

Develop
project8 · Feb 22, 2019 · 47eb583 · 47eb583
2 parents 74350e1 + f30a1fe
commit 47eb583
Show file tree

Hide file tree

Showing 19 changed files with 109 additions and 74 deletions.
diff --git a/CMakeLists.txt b/CMakeLists.txt
@@ -6,7 +6,7 @@ cmake_minimum_required( VERSION 3.1 )
 
 # Define the project
 cmake_policy( SET CMP0048 NEW ) # version in project()
-project( locust_mc VERSION 1.9.7)
+project( locust_mc VERSION 1.9.8)
 
 list( APPEND CMAKE_MODULE_PATH ${PROJECT_SOURCE_DIR}/Scarab/cmake )
 include( PackageBuilder )

diff --git a/Dockerfile b/Dockerfile
@@ -3,7 +3,7 @@ FROM project8/p8compute_dependencies:v0.4.0 as locust_common
 ARG build_type=Release
 ENV LOCUST_BUILD_TYPE=$build_type
 
-ENV LOCUST_TAG=v1.9.7
+ENV LOCUST_TAG=v1.9.8
 ENV LOCUST_BUILD_PREFIX=/usr/local/p8/locust/$LOCUST_TAG
 
 RUN mkdir -p $LOCUST_BUILD_PREFIX &&\

diff --git a/Scarab b/Scarab
diff --git a/Source/Core/LMCEggWriter.cc b/Source/Core/LMCEggWriter.cc
@@ -168,13 +168,9 @@ namespace locust
         return true;
     }
 
-    bool EggWriter::WriteRecord( const Signal* aSignal )
+    bool EggWriter::WriteRecord( const Signal* aSignal, bool is_new_acq )
     {
 
-
-        static bool t_is_new_acq = true;
-
-
         if( f_state != kWriting && f_state != kPrepared )
         {
             LERROR( lmclog, "Egg file must be opened before writing records" );
@@ -208,8 +204,8 @@ namespace locust
         ++f_record_id;
         f_record_time += f_record_length;
 
-        bool t_return = f_stream->WriteRecord( true ); // pls had to edit false to true
-        t_is_new_acq = false;
+        bool t_return = f_stream->WriteRecord( is_new_acq );
+        //f_is_new_acq = false;
         return t_return;
     }
 

diff --git a/Source/Core/LMCEggWriter.hh b/Source/Core/LMCEggWriter.hh
@@ -60,7 +60,7 @@ namespace locust
         public:
             bool PrepareEgg( const RunLengthCalculator* aRLC, const Digitizer* aDig );
 
-            bool WriteRecord( const Signal* aSignal );
+            bool WriteRecord( const Signal* aSignal, bool is_new_acq );
 
             bool FinalizeEgg();
 

diff --git a/Source/Core/LMCPowerCombiner.cc b/Source/Core/LMCPowerCombiner.cc
@@ -23,9 +23,11 @@ namespace locust
 
 
 
-    double PowerCombiner::GetVoltageDamping(int njunctions)
+    double PowerCombiner::GetVoltageDamping(int NPatchesPerStrip, unsigned z_index)
     {
     	// test:  clobber voltage with 0.6 at each junction.
+      int njunctions = fabs((double)z_index - (double)NPatchesPerStrip/2.);
+      if (z_index >= NPatchesPerStrip/2) njunctions += 1;
       return pow(1.0, njunctions);
     }
 
@@ -39,18 +41,38 @@ namespace locust
         return dphi;
     }
 
-  double PowerCombiner::GetCenterFedLinePhaseCorr(int NPatchesPerStrip, unsigned z_index, double DopplerFrequency, double PatchSpacing)
-    {
-        int njunctions = fabs((double)z_index - (double)NPatchesPerStrip/2.);
-        if (z_index >= NPatchesPerStrip/2) njunctions += 1;
-
-        double D = PatchSpacing; // m.  18.0 keV 90 degree electron, lambda in kapton.
-        double c_n = LMCConst::C()/1.704;  // speed of light in Kapton.
-        double lambda = c_n/DopplerFrequency;
-        double dphi = 2.*LMCConst::Pi() * D * njunctions / lambda;
-	//	printf("2PI*D/lambda is %g\n", 2.*LMCConst::Pi()*D/lambda); getchar();
-        return dphi;
-    }
+    double PowerCombiner::GetCenterFedLinePhaseCorr(int NPatchesPerStrip, unsigned z_index, double DopplerFrequency, double PatchSpacing)
+      {
+          int njunctions = fabs((double)z_index - (double)NPatchesPerStrip/2.);
+          if (z_index >= NPatchesPerStrip/2) njunctions += 1;
+//          printf("z_index is %d and njunctions is %d\n", z_index, njunctions); getchar();
+	  njunctions -= 1; //now, njunctions is the number of junctions between current patch and amplifier
+          double D = PatchSpacing; // m.  18.0 keV 90 degree electron, lambda in kapton.
+          double c_n = LMCConst::C()/1.704;  // speed of light in Kapton.
+          double lambda = c_n/DopplerFrequency;
+          double dphi = 2.*LMCConst::Pi() * D * njunctions / lambda;
+  	//	printf("2PI*D/lambda is %g\n", 2.*LMCConst::Pi()*D/lambda); getchar();
+          return dphi;
+      }
+
+
+    double PowerCombiner::GetCenterFedUnitCellDamping(int NPatchesPerStrip, unsigned z_index)
+      {
+          int njunctions = fabs((double)z_index - (double)NPatchesPerStrip/2.);
+          if (z_index >= NPatchesPerStrip/2) njunctions += 1;
+          double dampingfactor = 0.;
+
+          if (njunctions>1)
+            {
+	      dampingfactor = 0.38*0.66*pow(0.87, njunctions-1.); // patch, amplifier, njunction-1 regular junctions.
+            }
+          else
+            {
+	    dampingfactor = 0.38*0.66; //patch and amplifier, no junction
+            }
+
+          return dampingfactor;
+      }
 
 
 

diff --git a/Source/Core/LMCPowerCombiner.hh b/Source/Core/LMCPowerCombiner.hh
@@ -24,7 +24,8 @@ namespace locust
 
             double GetLinePhaseCorr(unsigned z_index, double DopplerFrequency);
       double GetCenterFedLinePhaseCorr(int NPatchesPerStrip, unsigned z_index, double DopplerFrequency, double PatchSpacing);
-            double GetVoltageDamping(int njunctions);
+            double GetVoltageDamping(int NPatchesPerStrip, unsigned z_index);
+            double GetCenterFedUnitCellDamping(int NPatchesPerStrip, unsigned z_index);
 
 
         private:

diff --git a/Source/Core/LMCSimulationController.cc b/Source/Core/LMCSimulationController.cc
@@ -116,6 +116,7 @@ namespace locust
         }
 
         bool IQStream = true;  // get rid of this.
+        bool isNewAcquisition = true;
         unsigned nRecords = fRunLengthCalc.GetNRecords();
         unsigned recordSize = fRunLengthCalc.GetRecordSize();
         unsigned nchannels = fRunLengthCalc.GetNChannels();
@@ -125,6 +126,7 @@ namespace locust
         for( unsigned record = 0; record < nRecords; ++record )
         {
             LINFO( lmclog, "Simulating record " << record );
+            LDEBUG( lmclog, "first generator: "<<fFirstGenerator->GetName());
             Signal* simulatedSignal = fFirstGenerator->Run( recordSize );
             if( simulatedSignal == NULL )
             {
@@ -172,12 +174,13 @@ namespace locust
             }
 
 
-            if( ! fEggWriter.WriteRecord( simulatedSignal ) )
+            if( ! fEggWriter.WriteRecord( simulatedSignal, isNewAcquisition ) )
             {
                 LERROR( lmclog, "Something went wrong while writing record " << record );
                 delete simulatedSignal;
                 return false;
             }
+            isNewAcquisition = false;
 
             // temporarily, immediately cleanup
 	    delete simulatedSignal;

diff --git a/Source/Generators/LMCGaussianNoiseGenerator.cc b/Source/Generators/LMCGaussianNoiseGenerator.cc
@@ -86,7 +86,7 @@ namespace locust
 
     void GaussianNoiseGenerator::SetMean( double aMean )
     {
-        SetMeanAndSigma( aMean, fSigma );
+        fMean = aMean;
         return;
     }
 
@@ -97,13 +97,13 @@ namespace locust
 
     void GaussianNoiseGenerator::SetSigma( double aSigma )
     {
-        SetMeanAndSigma( fMean, aSigma );
+        fSigma = aSigma;
         return;
     }
 
-    void GaussianNoiseGenerator::SetMeanAndSigma( double aMean, double aSigma )
+    void GaussianNoiseGenerator::SetMeanAndSigma( double aMean, double aSigma, double aSampledSigma )
     {
-        fNormDist = std::normal_distribution< double >( aMean, aSigma );
+        fNormDist = std::normal_distribution< double >( aMean, aSampledSigma );
         fUniDist = std::uniform_real_distribution< double >(0.,360.);
         fMean = aMean;
         fSigma = aSigma;
@@ -143,21 +143,23 @@ namespace locust
     bool GaussianNoiseGenerator::DoGenerateTime( Signal* aSignal )
     {
 
-    	SetMeanAndSigma( fMean, fSigma * sqrt(fAcquisitionRate * 1.e6) );
+    	SetMeanAndSigma( fMean, fSigma, fSigma * sqrt(fAcquisitionRate * 1.e6) );
 
         double gain=1.;
         const unsigned nchannels = fNChannels;
         double phi = 0.;  // voltage phase
-        double mag = 0.;  // voltage mag
+        double mag_r = 0.;  // voltage mag
+        double mag_i = 0.;
 
         for (int ch=0; ch<nchannels; ch++)
         {
             for( unsigned index = 0; index < aSignal->TimeSize(); ++index )
             {
                 phi = fUniDist( fRNG );
-                mag = fNormDist( fRNG );
-                aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][0] += gain*sqrt(50.)* mag * cos(phi*LMCConst::Pi()/180.);
-                aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][1] += gain*sqrt(50.)* mag * sin(phi*LMCConst::Pi()/180.);
+                mag_r = fNormDist( fRNG ) * sqrt(0.5);
+                mag_i = fNormDist( fRNG ) * sqrt(0.5);
+                aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][0] += gain*sqrt(50.)* mag_r;
+                aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][1] += gain*sqrt(50.)* mag_i;
                 //	    printf("noise signal is %g\n", aSignal->SignalTimeComplex()[ch*aSignal->TimeSize() + index][1]); getchar();
             }
         }

diff --git a/Source/Generators/LMCGaussianNoiseGenerator.hh b/Source/Generators/LMCGaussianNoiseGenerator.hh
@@ -57,7 +57,7 @@ namespace locust
             double GetSigma() const;
             void SetSigma( double aSigma );
 
-            void SetMeanAndSigma( double aMean, double aSigma );
+            void SetMeanAndSigma( double aMean, double aSigma, double aSampledSigma);
 
             Signal::State GetDomain() const;
             void SetDomain( Signal::State aDomain );

diff --git a/Source/Generators/LMCHighPassFilterFFTGenerator.cc b/Source/Generators/LMCHighPassFilterFFTGenerator.cc
@@ -21,7 +21,7 @@ namespace locust
     HighPassFilterFFTGenerator::HighPassFilterFFTGenerator( const std::string& aName ) :
         Generator( aName ),
         fDoGenerateFunc( &HighPassFilterFFTGenerator::DoGenerateTime ),
-        fThreshold ( 1.0e6 ) // 1 MHz
+        fThreshold ( 2.0e6 ) // MHz
     {
         fRequiredSignalState = Signal::kTime;
     }
@@ -67,10 +67,14 @@ namespace locust
 
     bool HighPassFilterFFTGenerator::DoGenerateTime( Signal* aSignal )
     {
+
         const unsigned nchannels = fNChannels;
-        double CutoffFreq = fThreshold;
-        int nwindows = 10;
-        int windowsize = aSignal->TimeSize()/nwindows;
+
+
+        double CutoffFreq = fThreshold; // Hz
+        double FastNyquist = fAcquisitionRate/2. * 1.e6 * aSignal->DecimationFactor();
+        int nwindows = 80;
+        int windowsize = aSignal->DecimationFactor()*aSignal->TimeSize()/nwindows;
 
 
         fftw_complex *SignalComplex;
@@ -90,20 +94,23 @@ namespace locust
                 // Construct complex voltage.
                 for( unsigned index = 0; index < windowsize; ++index )
                 {
-                    SignalComplex[index][0] = aSignal->SignalTimeComplex()[ ch*aSignal->TimeSize() + nwin*windowsize + index ][0];
-                    SignalComplex[index][1] = aSignal->SignalTimeComplex()[ ch*aSignal->TimeSize() + nwin*windowsize + index ][1];
+                    SignalComplex[index][0] = aSignal->LongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + nwin*windowsize + index ][0];
+                    SignalComplex[index][1] = aSignal->LongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + nwin*windowsize + index ][1];
+                    //if (index==20000) {printf("signal 20000 is %g\n", aSignal->SignalTime()[index]); getchar();}
                 }
 
                 fftw_execute(ForwardPlan);
 
+//                printf("windowsize/2 is %g and CutoffFreq is %g and FastNyquist is %g\n", windowsize/2, CutoffFreq, FastNyquist); getchar();
+
                 // High Pass FilterFFT
                 for( unsigned index = 0; index < windowsize; ++index )
                 {
-                    if ((index < CutoffFreq/(fAcquisitionRate*1.e6)*windowsize)||(index>0.5*windowsize))
-                    {
-                        FFTComplex[index][0] = 0.;
-                        FFTComplex[index][1] = 0.;
-                    }
+                if ((index < windowsize/2. * CutoffFreq/FastNyquist)||(index > windowsize/2.))
+                  {
+                  FFTComplex[index][0] = 0.;
+                  FFTComplex[index][1] = 0.;
+                  }
                 }
 
                 fftw_execute(ReversePlan);
@@ -113,16 +120,16 @@ namespace locust
                 for( unsigned index = 0; index < windowsize; ++index )
                 {
                     // normalize
-                    aSignal->SignalTimeComplex()[ ch*aSignal->TimeSize() + nwin*windowsize + index ][0] = SignalComplex[index][0]/norm;
-                    aSignal->SignalTimeComplex()[ ch*aSignal->TimeSize() + nwin*windowsize + index ][1] = SignalComplex[index][1]/norm;
+                    aSignal->LongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + nwin*windowsize + index ][0] = SignalComplex[index][0]/norm;
+                    aSignal->LongSignalTimeComplex()[ ch*aSignal->TimeSize()*aSignal->DecimationFactor() + nwin*windowsize + index ][1] = SignalComplex[index][1]/norm;
+                    //if (index>=20000) {printf("filtered signal is %g\n", aSignal->SignalTime()[index]); getchar();}
                 }
             }  // nwin
         }  // NCHANNELS
 
         delete [] SignalComplex;
         delete [] FFTComplex;
 
-
         return true;
     }
 

diff --git a/Source/Generators/LMCHighPassFilterFFTGenerator.hh b/Source/Generators/LMCHighPassFilterFFTGenerator.hh
@@ -20,7 +20,9 @@ namespace locust
      @brief Apply high pass filter to signal in frequency domain
 
      @details
-      This is a straightforward filter in frequency space.  It assumed that the signal has been sampled at fAcquisitionRate, and not fAcquisitionRate*aSignal->DecimationFactor().  The initial purpose of this filter is to exclude the DC bin.
+      This is a filter in frequency space.  It assumes that the signal has been
+      sampled at fAcquisitionRate*aSignal->DecimationFactor(), and needs to be
+      followed by the decimation generator.
       Available configuration options:
       - "threshold": double -- threshold below which signals are suppressed (Hz).
 

diff --git a/Source/Generators/LMCKassSignalGenerator.cc b/Source/Generators/LMCKassSignalGenerator.cc
@@ -28,6 +28,7 @@ namespace locust
         fLO_Frequency( 0.),
         gxml_filename("blank.xml"),
         gpitchangle_filename("blank.xml"),
+		fTruth( 0 ),
         phi_t1(0.),
         phi_t2(0.),
         phiLO_t(0.),
@@ -57,7 +58,10 @@ namespace locust
             gpitchangle_filename = aParam->get_value< std::string >( "pitchangle-filename" );
         }
 
-
+        if( aParam->has( "truth" ) )
+        {
+            fTruth = aParam->get_value< bool >( "truth" );
+        }
 
         return true;
     }
@@ -382,21 +386,13 @@ namespace locust
             {
                 if (ReceivedKassReady()) fPreEventInProgress = true;
                 printf("LMC says it ReceivedKassReady()\n");                
-                /*		if (Project8Phase == 2)
-                        {
-                        if ((index-StartEventTimer) > 1e6) 
-                        {
-                        break;  // if breaking after just one event
-                        }
-                        }
-                        */
-            }
+             }
 
             if (fPreEventInProgress)
             {
                 PreEventCounter += 1;
 
-                if (PreEventCounter > NPreEventSamples)  // finished noise samples.  Start event.
+                if (((!fTruth)&&(PreEventCounter > NPreEventSamples))||((fTruth)&&(PreEventCounter > NPreEventSamples)&&(index%(8192*aSignal->DecimationFactor())==0)  ))// finished pre-samples.  Start event.
                 {
                     fPreEventInProgress = false;  // reset.
                     fEventInProgress = true;

diff --git a/Source/Generators/LMCKassSignalGenerator.hh b/Source/Generators/LMCKassSignalGenerator.hh
@@ -52,6 +52,7 @@ namespace locust
 
         private:
             double fLO_Frequency;  // typically defined by a parameter in json file.
+            bool fTruth; // parameter in json file.  default is false.
             bool DoGenerate( Signal* aSignal );
             void* DriveAntenna(int PreEventCounter, unsigned index, Signal* aSignal, FILE *fp);
             int FindNode(double tNew) const;

diff --git a/Source/Generators/LMCPatchSignalGenerator.cc b/Source/Generators/LMCPatchSignalGenerator.cc
@@ -215,7 +215,7 @@ namespace locust
         	// assume 2PI delay between junctions, so we don't calculated phase mismatches.
         	// instead calculate damping on voltage amplitude:
             int njunctions = (int)fabs(z_index - fNPatchesPerStrip/2);
-            VoltageAmplitude *= aPowerCombiner.GetVoltageDamping(njunctions);
+            VoltageAmplitude *= aPowerCombiner.GetVoltageDamping(fNPatchesPerStrip, z_index);
         }
 
 	//        if (VoltageAmplitude>0.) {printf("voltageamplitude is %g\n", VoltageAmplitude); getchar();}
+1 −1		VERSION
+37 −0		documentation/cpp/Scarab_LibraryOverview.rst
+1 −0		documentation/cpp/requirements.txt
+26 −14		library/CMakeLists.txt
+9 −77		library/authentication/authentication.cc
+3 −3		library/authentication/authentication.hh
+0 −4		library/param/codec/yaml/param_yaml.cc
+6 −1		library/param/configurator.cc
+2 −1		library/param/param.cc
+1 −0		library/param/param.hh
+14 −2		library/test/CMakeLists.txt
+112 −0		library/test/test_digital.cc
+61 −18		library/utility/digital.hh