diff --git a/.github/workflows/build.yaml b/.github/workflows/build.yaml
index 592aedda..0fb64758 100644
--- a/.github/workflows/build.yaml
+++ b/.github/workflows/build.yaml
@@ -3,7 +3,7 @@ name: Build and Test Locust
on:
pull_request:
push:
- branches: [master, develop]
+ branches: [master, develop, feature/pileupChecks]
tags: ['*']
workflow_dispatch:
diff --git a/Config/CMakeLists.txt b/Config/CMakeLists.txt
index 29b8cd6c..dfe6cd98 100644
--- a/Config/CMakeLists.txt
+++ b/Config/CMakeLists.txt
@@ -5,9 +5,11 @@ configure_file( JustKass.xml.in JustKass.xml )
configure_file( JustKassFieldMap.xml.in JustKassFieldMap.xml )
configure_file( LocustKass_FreeSpace_Template.xml.in LocustKass_FreeSpace_Template.xml )
configure_file( LocustKass_Cavity_CCA.xml.in LocustKass_Cavity_CCA.xml )
+configure_file( LocustKass_Cavity_CCA_Pileup.xml.in LocustKass_Cavity_CCA_Pileup.xml )
configure_file( LocustKass_Cavity_1GHz.xml.in LocustKass_Cavity_1GHz.xml )
configure_file( LocustWaveguideTemplate.json.in LocustWaveguideTemplate.json )
configure_file( LocustCavityCCA.json.in LocustCavityCCA.json )
+configure_file( LocustCavityCCA_Pileup.json.in LocustCavityCCA_Pileup.json )
configure_file( LocustCavityCCA_ModeMapTest.json.in LocustCavityCCA_ModeMapTest.json )
configure_file( LocustCavityRectangular.json.in LocustCavityRectangular.json )
configure_file( LocustCavity1GHz.json.in LocustCavity1GHz.json )
@@ -27,6 +29,7 @@ set (EXAMPLE_CONFIGFILES
${CMAKE_CURRENT_BINARY_DIR}/JustKassFieldMap.xml
${CMAKE_CURRENT_BINARY_DIR}/LocustKass_FreeSpace_Template.xml
${CMAKE_CURRENT_BINARY_DIR}/LocustKass_Cavity_CCA.xml
+ ${CMAKE_CURRENT_BINARY_DIR}/LocustKass_Cavity_CCA_Pileup.xml
${CMAKE_CURRENT_BINARY_DIR}/LocustKass_Cavity_1GHz.xml
FreeSpaceGeometry.xml
CavityGeometry_VTinyCoil.xml
@@ -34,6 +37,7 @@ set (EXAMPLE_CONFIGFILES
${CMAKE_CURRENT_BINARY_DIR}/LocustWaveguideTemplate.json
${CMAKE_CURRENT_BINARY_DIR}/LocustCavityCCA.json
${CMAKE_CURRENT_BINARY_DIR}/LocustCavityCCA_ModeMapTest.json
+ ${CMAKE_CURRENT_BINARY_DIR}/LocustCavityCCA_Pileup.json
${CMAKE_CURRENT_BINARY_DIR}/LocustCavityRectangular.json
${CMAKE_CURRENT_BINARY_DIR}/LocustCavity1GHz.json
${CMAKE_CURRENT_BINARY_DIR}/LocustFreeSpaceTemplate.json
diff --git a/Config/LocustCavityCCA_Pileup.json.in b/Config/LocustCavityCCA_Pileup.json.in
new file mode 100644
index 00000000..44e3cb26
--- /dev/null
+++ b/Config/LocustCavityCCA_Pileup.json.in
@@ -0,0 +1,69 @@
+{
+ "generators":
+ [
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "cavity-signal",
+ "lpf-fft",
+ "decimate-signal",
+ "digitizer"
+ ],
+
+ "cavity-signal":
+ {
+ "transmitter": "kass-current",
+ "accumulate-truth-info": true,
+ "cavity-radius": 0.007,
+ "cavity-length": 0.1,
+ "back-reaction": "true",
+ "dho-cavity-frequency": 25.9e9,
+ "dho-time-resolution": 9.0e-11,
+ "dho-threshold-factor": 0.01,
+ "event-spacing-samples": 10,
+ "track-length": 0.5e-6,
+ "ks-starting-xpos-min": -0.005,
+ "ks-starting-xpos-max": 0.005,
+ "ks-starting-ypos-min": -0.005,
+ "ks-starting-ypos-max": 0.005,
+ "ks-starting-zpos-min": 0.0,
+ "ks-starting-zpos-max": 0.0,
+ "ks-starting-energy-min": 18600.0,
+ "ks-starting-energy-max": 18600.0,
+ "ks-starting-pitch-min": 87.0,
+ "ks-starting-pitch-max": 90.0,
+ "rectangular-waveguide": false,
+ "voltage-check": "false",
+ "lo-frequency": 25.9602e9,
+ "xml-filename": "${CMAKE_INSTALL_PREFIX}/config/LocustKass_Cavity_CCA_Pileup.xml"
+ },
+
+ "simulation":
+ {
+ "egg-filename": "${CMAKE_INSTALL_PREFIX}/output/locust_mc.egg",
+ "n-records": 1,
+ "record-size": 81920,
+ "acquisition-rate": 301,
+ "n-channels": 1
+ },
+
+ "gaussian-noise":
+ {
+ "noise-floor-psd": 2.76e-22,
+ "domain": "time"
+ },
+
+ "digitizer":
+ {
+ "v-range": 0.08,
+ "v-offset": -0.04
+ }
+
+}
+
diff --git a/Config/LocustKass_Cavity_CCA_Pileup.xml.in b/Config/LocustKass_Cavity_CCA_Pileup.xml.in
new file mode 100644
index 00000000..6f7d2e16
--- /dev/null
+++ b/Config/LocustKass_Cavity_CCA_Pileup.xml.in
@@ -0,0 +1,512 @@
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
diff --git a/Source/Fields/LMCCylindricalCavity.cc b/Source/Fields/LMCCylindricalCavity.cc
index 02aa52bc..000b06bb 100644
--- a/Source/Fields/LMCCylindricalCavity.cc
+++ b/Source/Fields/LMCCylindricalCavity.cc
@@ -88,6 +88,8 @@ namespace locust
fIntermediateFile = aParam["intermediate-file"]().as_bool();
}
+ scarab::path dataDir = aParam.get_value( "data-dir", ( TOSTRING(PB_DATA_INSTALL_DIR) ) );
+
if( aParam.has( "upload-modemap-filename" ) ) // import "realistic" test mode map
{
fFieldCore = new ModeMapCavity();
@@ -96,29 +98,19 @@ namespace locust
LERROR(lmclog,"There was a problem configuring the mode map.");
exit(-1);
}
- scarab::path dataDir = aParam.get_value( "data-dir", ( TOSTRING(PB_DATA_INSTALL_DIR) ) );
if (!fFieldCore->ReadModeMapTE_E((dataDir / aParam["upload-modemap-filename"]().as_string()).string(), aParam))
{
LERROR(lmclog,"There was a problem uploading the mode map.");
exit(-1);
}
- fFieldCore->ReadBesselZeroes((dataDir / "BesselZeros.txt").string(), 0 );
- fFieldCore->ReadBesselZeroes((dataDir / "BesselPrimeZeros.txt").string(), 1 );
- SetNormFactorsTE( SetUnityNormFactors(GetNModes()));
- SetNormFactorsTM( SetUnityNormFactors(GetNModes()));
+ SetNormFactors(SetUnityNormFactors(GetNModes(), 0)); // Temporary quick normalization factors of 1.0
}
else // otherwise default to ideal Pozar mode map
{
fFieldCore = new PozarCylindricalCavity();
- scarab::path dataDir = aParam.get_value( "data-dir", ( TOSTRING(PB_DATA_INSTALL_DIR) ) );
- fFieldCore->ReadBesselZeroes((dataDir / "BesselZeros.txt").string(), 0 );
- fFieldCore->ReadBesselZeroes((dataDir / "BesselPrimeZeros.txt").string(), 1 );
- SetNormFactorsTE(CalculateNormFactors(GetNModes(),1));
- SetNormFactorsTM(CalculateNormFactors(GetNModes(),0));
- CheckNormalization(GetNModes()); // E fields integrate to 1.0 for both TE and TM modes.
+ SetNormFactors(SetUnityNormFactors(GetNModes(), 0)); // Temporary quick normalization factors of 1.0
}
-
if( PlotModeMaps() )
{
double zSlice = 0.0;
@@ -129,47 +121,12 @@ namespace locust
PrintModeMaps(GetNModes(), zSlice, thetaSlice);
}
- return true;
- }
-
+ fFieldCore->ReadBesselZeroes((dataDir / "BesselZeros.txt").string(), 0 );
+ fFieldCore->ReadBesselZeroes((dataDir / "BesselPrimeZeros.txt").string(), 1 );
+ SetNormFactors(CalculateNormFactors(GetNModes(), 0)); // Calculate the realistic normalization factors.
+ CheckNormalization(GetNModes(), 0); // E fields integration over volume
- std::vector>> CylindricalCavity::CalculateNormFactors(int nModes, bool bTE)
- {
-
- LPROG(lmclog, "Calculating mode normalization factors ... " );
-
- std::vector>> aModeNormFactor;
- aModeNormFactor.resize(nModes);
-
- for (unsigned m=0; m CylindricalCavity::GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool teMode)
+ std::vector CylindricalCavity::GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool bTE)
{
double tR = tKassParticleXP[0];
double tTheta = tKassParticleXP[1];
double tZ = tKassParticleXP[2];
std::vector tField;
- double normFactor;
- if(teMode)
+ double normFactor = GetNormFactors()[bTE][l][m][n];
+
+ if(bTE)
{
tField = fFieldCore->TE_E(GetDimR(),2.*LMCConst::Pi(),GetDimL(),l,m,n,tR,tTheta,tZ,includeOtherPols);
- normFactor = GetNormFactorsTE()[l][m][n];
}
else
{
tField = fFieldCore->TM_E(GetDimR(),2.*LMCConst::Pi(),GetDimL(),l,m,n,tR,tTheta,tZ,includeOtherPols);
- normFactor = GetNormFactorsTM()[l][m][n];
}
+
auto it = tField.begin();
while (it != tField.end())
{
@@ -448,55 +405,6 @@ namespace locust
}
- void CylindricalCavity::CheckNormalization(int nModes)
- {
-
- printf("\n \\int{|E_xlm|^2 dV} = \\mu / \\epsilon \\int{|H_xlm|^2 dV} ?\n\n");
-
- for (int l=0; l GetDopplerFrequency(int l, int m, int n, std::vector tKassParticleXP);
virtual std::vector GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool teMode);
- virtual std::vector>> CalculateNormFactors(int nModes, bool bTE);
virtual std::vector GetTE_E(int l, int m, int n, double r, double theta, double z, bool includeOtherPols);
virtual std::vector GetTM_E(int l, int m, int n, double r, double theta, double z, bool includeOtherPols);
virtual double CalculateDotProductFactor(int l, int m, int n, std::vector tKassParticleXP, std::vector anE_normalized, double tThisEventNSamples);
virtual double GetDotProductFactor(std::vector tKassParticleXP, std::vector anE_normalized, bool IntermediateFile);
virtual bool InVolume(std::vector tKassParticleXP);
- virtual void CheckNormalization(int nModes);
virtual void PrintModeMaps(int nModes, double zSlice, double thetaSlice);
void PrintModeMapsLongSlice(int nModes, double thetaSlice);
virtual std::vector GetFieldAtProbe(int l, int m, int n, bool includeOtherPols, std::vector tKassParticleXP, bool teMode);
diff --git a/Source/Fields/LMCField.cc b/Source/Fields/LMCField.cc
index 925fef95..7a5f7666 100644
--- a/Source/Fields/LMCField.cc
+++ b/Source/Fields/LMCField.cc
@@ -135,58 +135,127 @@ namespace locust
return tModeSet;
}
-
- std::vector>> Field::SetUnityNormFactors(int nModes)
+ std::vector>>> Field::CalculateNormFactors(int nModes, bool bWaveguide)
{
- LPROG(lmclog, "Setting mode normalization factors to 1.0 ... " );
- std::vector>> tNormFactor;
- tNormFactor.resize(nModes);
+ LPROG(lmclog, "Calculating mode normalization factors ... " );
- for (unsigned m=0; m>>> aModeNormFactor;
+ aModeNormFactor.resize(2);
- for (unsigned l=0; l> tModeSet = ModeSelect(bWaveguide, 0);
- }
+ for (int mu=0; mu>> Field::GetNormFactorsTE()
+ std::vector>>> Field::SetUnityNormFactors(int nModes, bool bWaveguide)
{
- return fModeNormFactorTE;
+
+ std::vector>>> aModeNormFactor;
+ aModeNormFactor.resize(2);
+
+ for (int bTE=0; bTE<2; bTE++)
+ {
+ aModeNormFactor[bTE].resize(nModes);
+ for (unsigned l=0; l>> aNormFactor)
+
+ void Field::CheckNormalization(int nModes, bool bWaveguide)
{
- fModeNormFactorTE = aNormFactor;
+
+ printf("\n \\int{|E_xlm|^2 dV} = \\mu / \\epsilon \\int{|H_xlm|^2 dV} ?\n\n");
+
+ std::vector> tModeSet = ModeSelect(bWaveguide, 0);
+
+ for (int mu=0; mu>> Field::GetNormFactorsTM()
+
+ std::vector>>> Field::GetNormFactors()
{
- return fModeNormFactorTM;
+ return fModeNormFactor;
}
- void Field::SetNormFactorsTM(std::vector>> aNormFactor)
+ void Field::SetNormFactors(std::vector>>> aNormFactor)
{
- fModeNormFactorTM = aNormFactor;
+ fModeNormFactor = aNormFactor;
}
std::vector>> Field::GetAvgDotProductFactor()
diff --git a/Source/Fields/LMCField.hh b/Source/Fields/LMCField.hh
index f6d68e07..2924e8fb 100644
--- a/Source/Fields/LMCField.hh
+++ b/Source/Fields/LMCField.hh
@@ -24,8 +24,6 @@ namespace locust
@details
Available configuration options:
- "n-modes" : int -- [2] Range of l, m, and n indices used to configure available mode normalizations.
- Of the available normalized modes, modes to be included in the simulation itself are identified in the function
- LMCFieldCalculator::ModeSelect().
- "n-pixels" : int -- [100] Number of pixels used in each dimension of the mode field definitions.
- "plot-mode-maps": bool -- [false] Option to print all normalized (see n-modes above) mode maps to
2D histograms in Root files, for inspection.
@@ -88,23 +86,21 @@ namespace locust
virtual std::vector GetDopplerFrequency(int l, int m, int n, std::vector tKassParticleXP) {return {0.};};
virtual std::vector GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool teMode) {return {0.};};
double NormalizedEFieldMag(std::vector field);
- virtual std::vector>> CalculateNormFactors(int nModes, bool bTE) {return {{{0.}}};};
- std::vector>> SetUnityNormFactors(int nModes);
+ std::vector>>> CalculateNormFactors(int nModes, bool bWaveguide);
+ std::vector>>> SetUnityNormFactors(int nModes, bool bWaveguide);
virtual std::vector GetTE_E(int l, int m, int n, double r, double theta, double z, bool includeOtherPols) {return {0.};};
virtual std::vector GetTM_E(int l, int m, int n, double r, double theta, double z, bool includeOtherPols) {return {0.};};
virtual double CalculateDotProductFactor(int l, int m, int n, std::vector tKassParticleXP, std::vector aTE_E_normalized, double tThisEventNSamples) {return {0.};};
virtual double GetDotProductFactor(std::vector tKassParticleXP, std::vector aTE_E_normalized, bool IntermediateFile) {return {0.};};
virtual bool InVolume(std::vector tKassParticleXP){return false;};
- virtual void CheckNormalization(int nModes){};
+ void CheckNormalization(int nModes, bool bWaveguide);
virtual void PrintModeMaps(int nModes, double zSlice, double thetaSlice){};
virtual std::vector GetFieldAtProbe(int l, int m, int n, bool includeOtherPols, std::vector tKassParticleXP, bool teMode){return {0.};};
virtual double ScaleEPoyntingVector(double fcyc){return 0.;};
std::vector> ModeSelect(bool bWaveguide, bool bNormCheck);
- std::vector>> GetNormFactorsTE();
- void SetNormFactorsTE(std::vector>> aNormFactor);
- std::vector>> GetNormFactorsTM();
- void SetNormFactorsTM(std::vector>> aNormFactor);
+ std::vector>>> GetNormFactors();
+ void SetNormFactors(std::vector>>> aNormFactor);
std::vector>> GetAvgDotProductFactor();
void SetAvgDotProductFactor(std::vector>> aFactor);
double GetCentralFrequency();
@@ -134,8 +130,7 @@ namespace locust
private:
int fNModes;
- std::vector>> fModeNormFactorTE; // 3D vector [n-modes][n-modes][n-modes].
- std::vector>> fModeNormFactorTM; // 3D vector [n-modes][n-modes][n-modes].
+ std::vector>>> fModeNormFactor; // 4D vector [2][n-modes][n-modes][n-modes].
double fCentralFrequency;
int fnPixels;
double fR; // Cylindrical cavity dimenions.
diff --git a/Source/Fields/LMCModeMapCavity.cc b/Source/Fields/LMCModeMapCavity.cc
index 8cf7492f..6d80212d 100644
--- a/Source/Fields/LMCModeMapCavity.cc
+++ b/Source/Fields/LMCModeMapCavity.cc
@@ -116,15 +116,29 @@ namespace locust
std::string token;
std::stringstream ss(lineContent);
int wordCount = 0;
+ double r, theta, z;
int i,j,k;
- double Erho,Etheta,Ez;
+ double Erho,Etheta;
+ double Ex, Ey, Ez;
while (ss >> token)
{
- if (wordCount == 0) i = (int)((std::stod(token)-fDim1_min)/(fDim1_max-fDim1_min)*(fnPixel1)); // var1 position
- else if (wordCount == 1) j = (int)((std::stod(token)-fDim2_min)/(fDim2_max-fDim2_min)*(fnPixel2)); // var2 position
- else if (wordCount == 2) k = (int)((std::stod(token)-fDim3_min)/(fDim3_max-fDim3_min)*(fnPixel3)); // var3 position
- else if (wordCount == 3) Erho = std::stod(token); // mode E field value
- else if (wordCount == 4) Etheta = std::stod(token); // mode E field value
+ if (wordCount == 0)
+ {
+ r = std::stod(token);
+ i = (int)((r-fDim1_min)/(fDim1_max-fDim1_min)*(fnPixel1)); // var1 position
+ }
+ else if (wordCount == 1)
+ {
+ theta = std::stod(token);
+ j = (int)((theta-fDim2_min)/(fDim2_max-fDim2_min)*(fnPixel2)); // var2 position
+ }
+ else if (wordCount == 2)
+ {
+ z = std::stod(token);
+ k = (int)((z-fDim3_min)/(fDim3_max-fDim3_min)*(fnPixel3)); // var3 position
+ }
+ else if (wordCount == 3) Ex = std::stod(token); // mode E field value
+ else if (wordCount == 4) Ey = std::stod(token); // mode E field value
else if (wordCount == 5) Ez = std::stod(token); // mode E field value
else
{
@@ -134,11 +148,28 @@ namespace locust
++wordCount;
}
+ if ((i==fnPixel1) or (j==fnPixel2) or (k==fnPixel3))
+ {
+ continue;
+ }
if ((i>=fnPixel1) or (j>=fnPixel2) or (k>=fnPixel3))
{
LERROR(lmclog,"Imported mode map dimensions don't agree with those in \"" << aFilename <<".\" Double check dim[1,2,3]-max.");
return false;
}
+
+ //Must convert E field from cartesian coordinates to cylindrical coordinates
+ if(r<1.e-10)
+ {
+ Erho = 0.;
+ Etheta = 0.;
+ }
+ else
+ {
+ Erho = ((Ex * r*cos(theta)) + Ey * r*sin(theta)) / r;
+ Etheta = ((Ey * r*cos(theta)) - Ex * r*sin(theta)) / r;
+ }
+
std::vector E_input = {Erho,Etheta,Ez};
fModeMapTE_E[i][j][k] = E_input;
// printf("read var1 is %g, var2 is %g, E is %g\n", fModeMapTE_E.back()[0], fModeMapTE_E.back()[1], fModeMapTE_E.back()[2]);
diff --git a/Source/Fields/LMCRectangularCavity.cc b/Source/Fields/LMCRectangularCavity.cc
index 8e36e218..09275388 100644
--- a/Source/Fields/LMCRectangularCavity.cc
+++ b/Source/Fields/LMCRectangularCavity.cc
@@ -96,17 +96,15 @@ namespace locust
LERROR(lmclog,"There was a problem uploading the mode map.");
exit(-1);
}
- SetNormFactorsTE( SetUnityNormFactors(GetNModes()));
- SetNormFactorsTM( SetUnityNormFactors(GetNModes()));
}
else
{
fFieldCore = new PozarRectangularCavity();
- SetNormFactorsTE(CalculateNormFactors(GetNModes(),1));
- SetNormFactorsTM(CalculateNormFactors(GetNModes(),0));
- CheckNormalization(GetNModes()); // E fields integrate to 1.0 for both TE and TM modes.
}
+ SetNormFactors(CalculateNormFactors(GetNModes(), 0));
+ CheckNormalization(GetNModes(), 0); // E fields integrate to 1.0 for both TE and TM modes.
+
if( PlotModeMaps() )
{
@@ -119,46 +117,6 @@ namespace locust
return true;
}
- std::vector>> RectangularCavity::CalculateNormFactors(int nModes, bool bTE)
- {
-
- LPROG(lmclog, "Calculating mode normalization factors ... " );
-
- std::vector>> aModeNormFactor;
- aModeNormFactor.resize(nModes);
-
- for (unsigned m=0; m RectangularCavity::GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool teMode)
+ std::vector RectangularCavity::GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool bTE)
{
double tR = tKassParticleXP[0];
double tTheta = tKassParticleXP[1];
@@ -303,17 +261,18 @@ namespace locust
double tX = tR*cos(tTheta);
double tY = tR*sin(tTheta);
std::vector tField;
- double normFactor;
- if(teMode)
+
+ double normFactor = GetNormFactors()[bTE][l][m][n];
+
+ if(bTE)
{
tField = fFieldCore->TE_E(GetDimX(),GetDimY(),GetDimL(),l,m,n,tX,tY,tZ,0);
- normFactor = GetNormFactorsTE()[l][m][n];
}
else
{
tField = fFieldCore->TM_E(GetDimX(),GetDimY(),GetDimL(),l,m,n,tX,tY,tZ,0);
- normFactor = GetNormFactorsTM()[l][m][n];
}
+
auto it = tField.begin();
while (it != tField.end())
{
@@ -400,56 +359,6 @@ namespace locust
}
- void RectangularCavity::CheckNormalization(int nModes)
- {
-
- printf("\n \\int{|E_xlm|^2 dV} = \\mu / \\epsilon \\int{|H_xlm|^2 dV} ?\n\n");
-
- for (int l=0; l GetDopplerFrequency(int l, int m, int n, std::vector tKassParticleXP);
virtual std::vector GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool teMode);
- virtual std::vector>> CalculateNormFactors(int nModes, bool bTE);
virtual std::vector GetTE_E(int l, int m, int n, double r, double theta, double z, bool includeOtherPols);
virtual std::vector GetTM_E(int l, int m, int n, double r, double theta, double z, bool includeOtherPols);
virtual double CalculateDotProductFactor(int l, int m, int n, std::vector tKassParticleXP, std::vector anE_normalized, double tThisEventNSamples);
virtual double GetDotProductFactor(std::vector tKassParticleXP, std::vector anE_normalized, bool IntermediateFile);
virtual bool InVolume(std::vector tKassParticleXP);
- virtual void CheckNormalization(int nModes);
virtual void PrintModeMaps(int nModes, double zSlice);
void PrintModeMapsLongSlice(int nModes, double xSlice);
virtual std::vector GetFieldAtProbe(int l, int m, int n, bool includeOtherPols, std::vector tKassParticleXP, bool teMode);
diff --git a/Source/Fields/LMCRectangularWaveguide.cc b/Source/Fields/LMCRectangularWaveguide.cc
index 7412d214..eb5b7c27 100644
--- a/Source/Fields/LMCRectangularWaveguide.cc
+++ b/Source/Fields/LMCRectangularWaveguide.cc
@@ -59,10 +59,8 @@ namespace locust
fFieldCore = new PozarRectangularWaveguide();
- SetNormFactorsTE(CalculateNormFactors(GetNModes(),1));
- SetNormFactorsTM(CalculateNormFactors(GetNModes(),0));
-
- CheckNormalization(GetNModes()); // E fields integrate to 1.0 for both TE and TM modes.
+ SetNormFactors(CalculateNormFactors(GetNModes(), 1));
+ CheckNormalization(GetNModes(), 1); // E fields integrate to 1.0 for both TE and TM modes.
if( PlotModeMaps() )
{
@@ -223,7 +221,7 @@ namespace locust
double tY = tKassParticleXP[0] * sin(tKassParticleXP[1]);
double fcyc = tKassParticleXP[7];
std::vector tTE_E_electron = fFieldCore->TE_E(GetDimX(),GetDimY(),m,n,tX,tY,fcyc);
- double normFactor = GetNormFactorsTE()[l][m][n];
+ double normFactor = GetNormFactors()[bTE][l][m][n];
auto it = tTE_E_electron.begin();
while (it != tTE_E_electron.end())
@@ -282,100 +280,7 @@ namespace locust
return unitJdotE;
}
- std::vector>> RectangularWaveguide::CalculateNormFactors(int nModes, bool bTE)
- {
-
- LPROG( lmclog, "Calculating mode normalization factors ... " );
-
- std::vector>> aModeNormFactor;
- aModeNormFactor.resize(nModes);
-
- for (unsigned m=0; m tKassParticleXP)
{
@@ -430,16 +335,9 @@ namespace locust
TH2D* hTEy = new TH2D(hname_y, haxis_y, nbins, -GetDimX()/2., GetDimX()/2., nbins, -GetDimY()/2., GetDimY()/2.);
TH2D* hTEx = new TH2D(hname_x, haxis_x, nbins, -GetDimX()/2., GetDimX()/2., nbins, -GetDimY()/2., GetDimY()/2.);
- double normFactor = 1.0;
+ double normFactor = GetNormFactors()[bTE][l][m][n];
int a = 0;
- if (bTE)
- {
- normFactor = GetNormFactorsTE()[l][m][n];
- }
- else
- {
- normFactor = GetNormFactorsTM()[l][m][n];
- }
+
for (unsigned i=0; i GetDopplerFrequency(int l, int m, int n, std::vector tKassParticleXP);
virtual std::vector GetNormalizedModeField(int l, int m, int n, std::vector tKassParticleXP, bool includeOtherPols, bool bTE);
- virtual std::vector>> CalculateNormFactors(int nModes, bool bTE);
virtual double CalculateDotProductFactor(int l, int m, int n, std::vector tKassParticleXP, std::vector anE_normalized, double tThisEventNSamples);
virtual double GetDotProductFactor(std::vector tKassParticleXP, std::vector aTE_E_normalized, bool IntermediateFile);
virtual bool InVolume(std::vector tKassParticleXP);
- virtual void CheckNormalization(int nModes);
virtual void PrintModeMaps(int nModes, double zSlice, double thetaSlice);
double GetGroupVelocity(int m, int n, double fcyc);
double ScaleEPoyntingVector(double fcyc);
diff --git a/Source/Generators/LMCLowPassFilterFFTGenerator.cc b/Source/Generators/LMCLowPassFilterFFTGenerator.cc
index 76b662a0..69c7dce7 100644
--- a/Source/Generators/LMCLowPassFilterFFTGenerator.cc
+++ b/Source/Generators/LMCLowPassFilterFFTGenerator.cc
@@ -103,6 +103,8 @@ namespace locust
bool LowPassFilterFFTGenerator::DoGenerateTime( Signal* aSignal )
{
+ LWARN(lmclog,"Calculating LPF-FFT ...\n");
+
const unsigned nchannels = fNChannels;
double CutoffFreq = fThreshold * fAcquisitionRate/2. * 1.e6; // Hz
@@ -120,8 +122,6 @@ namespace locust
if (variableWindowSize > 1)
{
- LWARN(lmclog,"Calculating LPF-FFT ...\n");
-
fftw_complex *SignalComplex;
SignalComplex = (fftw_complex*)fftw_malloc( sizeof(fftw_complex) * variableWindowSize );
fftw_complex *FFTComplex;
diff --git a/Source/IO/LMCEvent.cc b/Source/IO/LMCEvent.cc
index d8b279c1..453436f9 100644
--- a/Source/IO/LMCEvent.cc
+++ b/Source/IO/LMCEvent.cc
@@ -30,6 +30,7 @@ namespace locust
fStartFrequencies.push_back( aTrack.StartFrequency );
fEndFrequencies.push_back( aTrack.EndFrequency );
fAvgFrequencies.push_back( aTrack.AvgFrequency );
+ fOutputAvgFrequencies.push_back( aTrack.OutputAvgFrequency );
fAvgAxialFrequencies.push_back( aTrack.AvgAxialFrequency );
fTrackPowers.push_back( aTrack.TrackPower );
fStartTimes.push_back( aTrack.StartTime );
diff --git a/Source/IO/LMCEvent.hh b/Source/IO/LMCEvent.hh
index 8559057c..2202e2d6 100644
--- a/Source/IO/LMCEvent.hh
+++ b/Source/IO/LMCEvent.hh
@@ -38,6 +38,7 @@ namespace locust
std::vector fStartFrequencies;
std::vector fEndFrequencies;
std::vector fAvgFrequencies;
+ std::vector fOutputAvgFrequencies;
std::vector fAvgAxialFrequencies;
std::vector fTrackPowers;
std::vector fStartTimes;
diff --git a/Source/IO/LMCRootTreeWriter.cc b/Source/IO/LMCRootTreeWriter.cc
index 8926a1da..a0de1899 100644
--- a/Source/IO/LMCRootTreeWriter.cc
+++ b/Source/IO/LMCRootTreeWriter.cc
@@ -61,6 +61,7 @@ namespace locust
aTree->Branch("StartFrequencies", "std::vector", &anEvent->fStartFrequencies);
aTree->Branch("EndFrequencies", "std::vector", &anEvent->fEndFrequencies);
aTree->Branch("AvgFrequencies", "std::vector", &anEvent->fAvgFrequencies);
+ aTree->Branch("OutputAvgFrequencies", "std::vector", &anEvent->fOutputAvgFrequencies);
aTree->Branch("AvgAxialFrequencies", "std::vector", &anEvent->fAvgAxialFrequencies);
aTree->Branch("StartTimes", "std::vector", &anEvent->fStartTimes);
aTree->Branch("EndTimes", "std::vector", &anEvent->fEndTimes);
diff --git a/Source/IO/LMCTrack.cc b/Source/IO/LMCTrack.cc
index 50063209..390cdd54 100644
--- a/Source/IO/LMCTrack.cc
+++ b/Source/IO/LMCTrack.cc
@@ -31,7 +31,9 @@ namespace locust
OutputStartFrequency = -99.;
StartFrequency = -99.;
EndFrequency = -99.;
- AvgFrequency = -99.;
+ AvgFrequency = 0.;
+ OutputAvgFrequency = 0.;
+ AvgAxialFrequency = 0.;
LOFrequency = -99.;
TrackPower = -99.;
Slope = -99.;
diff --git a/Source/IO/LMCTrack.hh b/Source/IO/LMCTrack.hh
index 96a6edfa..353bd01f 100644
--- a/Source/IO/LMCTrack.hh
+++ b/Source/IO/LMCTrack.hh
@@ -39,6 +39,7 @@ namespace locust
double StartFrequency = -99.;
double EndFrequency = -99.;
double AvgFrequency = -99.;
+ double OutputAvgFrequency = -99.;
double AvgAxialFrequency = -99.;
double LOFrequency = -99.;
double TrackPower = -99.;
diff --git a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc
index 83c9dbf3..78612273 100644
--- a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc
+++ b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.cc
@@ -13,7 +13,6 @@ namespace locust
fNewParticleHistory(),
fFieldCalculator( NULL ),
fPitchAngle( -99. ),
- fLMCTrackID( -2 ),
fT0trapMin( 0. ),
fNCrossings( 0 ),
fSampleIndex( 0 ),
@@ -25,7 +24,6 @@ namespace locust
fNewParticleHistory(),
fFieldCalculator( NULL ),
fPitchAngle( aCopy.fPitchAngle ),
- fLMCTrackID( aCopy.fLMCTrackID ),
fT0trapMin( aCopy.fT0trapMin ),
fNCrossings( aCopy.fNCrossings ),
fSampleIndex( aCopy.fSampleIndex ),
@@ -64,7 +62,10 @@ namespace locust
bool CyclotronRadiationExtractor::UpdateTrackProperties( Kassiopeia::KSParticle &aFinalParticle, unsigned index, bool bStart )
{
- double tTime = index / fInterface->aRunParameter->fSamplingRateMHz / 1.e6 / fInterface->aRunParameter->fDecimationFactor;
+ double tLOfrequency = fInterface->aRunParameter->fLOfrequency; // Hz
+ double tSamplingRate = fInterface->aRunParameter->fSamplingRateMHz; // MHz
+ double tOffset = -tLOfrequency + tSamplingRate * 1.e6 / 2.; // Hz
+ double tTime = index / tSamplingRate / 1.e6 / fInterface->aRunParameter->fDecimationFactor;
#ifdef ROOT_FOUND
if (bStart)
{
@@ -75,13 +76,13 @@ namespace locust
fInterface->aTrack.Radius = pow(tX*tX + tY*tY, 0.5);
fInterface->aTrack.RadialPhase = calcOrbitPhase(tX, tY);
fInterface->aTrack.StartingEnergy_eV = LMCConst::kB_eV() / LMCConst::kB() * aFinalParticle.GetKineticEnergy();
- std::this_thread::sleep_for(std::chrono::milliseconds(1000));
}
else
{
fInterface->aTrack.EndTime = tTime;
unsigned nElapsedSamples = index - fStartingIndex;
fInterface->aTrack.AvgFrequency = ( fInterface->aTrack.AvgFrequency * nElapsedSamples + aFinalParticle.GetCyclotronFrequency() ) / ( nElapsedSamples + 1);
+ fInterface->aTrack.OutputAvgFrequency = fInterface->aTrack.AvgFrequency + tOffset;
fInterface->aTrack.TrackLength = tTime - fInterface->aTrack.StartTime;
fInterface->aTrack.Slope = (fInterface->aTrack.EndFrequency - fInterface->aTrack.StartFrequency) / (fInterface->aTrack.TrackLength);
}
@@ -213,11 +214,13 @@ namespace locust
if (!fInterface->fDoneWithSignalGeneration) // if Locust is still acquiring voltages.
{
-
- if ( aFinalParticle.GetIndexNumber() > fLMCTrackID ) // check for new track
+ // Check for either a new track, or a new event:
+ if ( ( fInterface->fNewTrackStarting ) || ( aFinalParticle.GetParentTrackId() < 0 ) )
{
- fLMCTrackID = aFinalParticle.GetIndexNumber();
- fPitchAngle = -99.; // new electron needs central pitch angle reset.
+ std::this_thread::sleep_for(std::chrono::milliseconds(1000));
+ fInterface->fNewTrackStarting = false;
+ aFinalParticle.SetParentTrackId( aFinalParticle.GetParentTrackId() + 1 );
+ fPitchAngle = -99.; // new track needs central pitch angle reset.
double dt = aFinalParticle.GetTime() - anInitialParticle.GetTime();
fFieldCalculator->SetNFilterBinsRequired( dt );
UpdateTrackProperties( aFinalParticle, fInterface->fSampleIndex, 1 );
diff --git a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh
index f1160ed0..3581145a 100644
--- a/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh
+++ b/Source/Kassiopeia/LMCCyclotronRadiationExtractor.hh
@@ -57,7 +57,6 @@ namespace locust
private:
std::deque fNewParticleHistory;
double fPitchAngle;
- int fLMCTrackID;
double fT0trapMin;
int fNCrossings;
FieldCalculator* fFieldCalculator;
diff --git a/Source/Kassiopeia/LMCKassLocustInterface.cc b/Source/Kassiopeia/LMCKassLocustInterface.cc
index a31f5d7a..23cab136 100644
--- a/Source/Kassiopeia/LMCKassLocustInterface.cc
+++ b/Source/Kassiopeia/LMCKassLocustInterface.cc
@@ -39,7 +39,8 @@ namespace locust
fbWaveguide( false ),
fSampleIndex( 0 ),
fTriggerConfirm( 100000 ),
- fFastRecordLength( 0 )
+ fFastRecordLength( 0 ),
+ fNewTrackStarting( false )
{}
KLInterfaceBootstrapper::KLInterfaceBootstrapper() :
diff --git a/Source/Kassiopeia/LMCKassLocustInterface.hh b/Source/Kassiopeia/LMCKassLocustInterface.hh
index 7e5c3287..2ef3aad2 100644
--- a/Source/Kassiopeia/LMCKassLocustInterface.hh
+++ b/Source/Kassiopeia/LMCKassLocustInterface.hh
@@ -71,6 +71,7 @@ namespace locust
unsigned fSampleIndex;
int fTriggerConfirm;
int fFastRecordLength;
+ bool fNewTrackStarting;
#ifdef ROOT_FOUND
Event* anEvent;
diff --git a/Source/Kassiopeia/LMCRunPause.cc b/Source/Kassiopeia/LMCRunPause.cc
index 16b3e994..043225eb 100644
--- a/Source/Kassiopeia/LMCRunPause.cc
+++ b/Source/Kassiopeia/LMCRunPause.cc
@@ -24,8 +24,31 @@ namespace locust
RunPause::RunPause() :
fToolbox(KToolbox::GetInstance()),
KSComponent(),
+ fLocustMaxTimeTerminator( nullptr ),
+ fLocustMaxRTerminator( nullptr ),
+ fBox( nullptr ),
+ fKGSpace( nullptr ),
+ fSurface( nullptr ),
+ fLocustTermDeath( nullptr ),
+ fCommand( nullptr ),
+ fKSSpace( nullptr ),
+ fGenDirectionComposite( nullptr ),
+ fThetaGenerator( nullptr ),
+ fPhiGenerator( nullptr ),
+ fGenPositionComposite( nullptr ),
+ fPositionXGenerator( nullptr ),
+ fPositionYGenerator( nullptr ),
+ fPositionZGenerator( nullptr ),
+ fGenEnergyComposite( nullptr ),
+ fEnergyGenerator( nullptr ),
+ fTimeGenerator( nullptr ),
+ fGenTimeComposite( nullptr ),
+ fGenPidComposite( nullptr ),
+ fGenerator( nullptr ),
fMinTrackLengthFraction(0.1),
- fConfigurationComplete( false ),
+ fConfigurationComplete( false ),
+ fEventCounter( 0 ),
+ fMaxEvents( 1 ),
fInterface( KLInterfaceBootstrapper::get_instance()->GetInterface() )
{
}
@@ -33,8 +56,31 @@ namespace locust
RunPause::RunPause( const RunPause& aCopy ) :
fToolbox(KToolbox::GetInstance()),
KSComponent(),
+ fLocustMaxTimeTerminator( nullptr ),
+ fLocustMaxRTerminator( nullptr ),
+ fBox( nullptr ),
+ fKGSpace( nullptr ),
+ fSurface( nullptr ),
+ fLocustTermDeath( nullptr ),
+ fCommand( nullptr ),
+ fKSSpace( nullptr ),
+ fGenDirectionComposite( nullptr ),
+ fThetaGenerator( nullptr ),
+ fPhiGenerator( nullptr ),
+ fGenPositionComposite( nullptr ),
+ fPositionXGenerator( nullptr ),
+ fPositionYGenerator( nullptr ),
+ fPositionZGenerator( nullptr ),
+ fGenEnergyComposite( nullptr ),
+ fEnergyGenerator( nullptr ),
+ fTimeGenerator( nullptr ),
+ fGenTimeComposite( nullptr ),
+ fGenPidComposite( nullptr),
+ fGenerator( nullptr ),
fMinTrackLengthFraction(0.1),
- fConfigurationComplete( false ),
+ fConfigurationComplete( false ),
+ fEventCounter( 0 ),
+ fMaxEvents( 1 ),
fInterface( aCopy.fInterface )
{
}
@@ -66,13 +112,20 @@ namespace locust
LPROG(lmclog,"RunPause class did not need to be configured.");
return true;
}
+
fConfigurationComplete = true;
+ LPROG(lmclog,"RunPause has been configured.");
}
+
return true;
}
bool RunPause::Configure( const scarab::param_node& aParam )
{
+ if( aParam.has( "n-max-events" ) )
+ {
+ fMaxEvents = aParam["n-max-events"]().as_double();
+ }
if ( fToolbox.IsInitialized() )
{
@@ -121,105 +174,115 @@ namespace locust
bool RunPause::AddGenerator( const scarab::param_node& aParam )
{
- if ( aParam.has( "ks-starting-xpos-min" ) && aParam.has( "ks-starting-xpos-max" )
- && aParam.has( "ks-starting-ypos-min" ) && aParam.has( "ks-starting-ypos-max" )
- && aParam.has( "ks-starting-zpos-min" ) && aParam.has( "ks-starting-zpos-max" )
- && aParam.has( "ks-starting-energy-min" ) && aParam.has( "ks-starting-energy-max" )
- && aParam.has( "ks-starting-pitch-min" ) && aParam.has( "ks-starting-pitch-max" ) )
+ if (!fToolbox.HasKey("gen_project8"))
{
- KRandom::GetInstance().SetSeed(time(NULL));
- auto tGen = fToolbox.GetAll();
- for (unsigned i=0; iIsActivated()) && (tGen[i]->GetName()!="root_generator") )
+ KRandom::GetInstance().SetSeed(time(NULL));
+ auto tGen = fToolbox.GetAll();
+ for (unsigned i=0; iGetName());
- fToolbox.Get("root_generator")->ClearGenerator(tGen[i]);
- fToolbox.Remove(tGen[i]->GetName());
+ if ( (tGen[i]->IsActivated()) && (tGen[i]->GetName()!="root_generator") )
+ {
+ LPROG(lmclog,"Clearing " << tGen[i]->GetName() << " from KSRoot ... ");
+ fToolbox.Get("root_generator")->ClearGenerator(tGen[i]);
+ tGen[i]->RemoveTags(tGen[i]->GetTags());
+ tGen[i]->Deactivate();
+ tGen[i]->Deinitialize();
+ fToolbox.Remove(tGen[i]->GetName());
+ }
}
- }
- Kassiopeia::KSGenValueFix* tGenPidComposite = new Kassiopeia::KSGenValueFix();
- tGenPidComposite->SetValue(11); // electron
+ if ( fGenPidComposite == nullptr ) fGenPidComposite = new Kassiopeia::KSGenValueFix();
+ fGenPidComposite->SetValue(11); // electron
- Kassiopeia::KSGenTimeComposite* tGenTimeComposite = new Kassiopeia::KSGenTimeComposite();
- Kassiopeia::KSGenValueUniform* tTimeGenerator = new Kassiopeia::KSGenValueUniform();
- tTimeGenerator->SetValueMin(0.);
- tTimeGenerator->SetValueMax(0.);
- tGenTimeComposite->SetTimeValue(tTimeGenerator);
+ if ( fGenTimeComposite == nullptr ) fGenTimeComposite = new Kassiopeia::KSGenTimeComposite();
+ if ( fTimeGenerator == nullptr ) fTimeGenerator = new Kassiopeia::KSGenValueUniform();
+ fTimeGenerator->SetValueMin(0.);
+ fTimeGenerator->SetValueMax(0.);
+ fGenTimeComposite->SetTimeValue(fTimeGenerator);
+
+ if ( fGenEnergyComposite == nullptr ) fGenEnergyComposite = new Kassiopeia::KSGenEnergyComposite();
+ if ( fEnergyGenerator == nullptr ) fEnergyGenerator = new Kassiopeia::KSGenValueUniform();
+ if ( aParam.has( "ks-starting-energy-min" ) && ( aParam.has( "ks-starting-energy-max" ) ) )
+ {
+ fEnergyGenerator->SetValueMin( aParam["ks-starting-energy-min"]().as_double() ); // eV
+ fEnergyGenerator->SetValueMax( aParam["ks-starting-energy-max"]().as_double() ); // eV
+ }
+ else
+ {
+ fEnergyGenerator->SetValueMin( 18600. ); // eV
+ fEnergyGenerator->SetValueMax( 18600. ); // eV
+ }
+
+ fGenEnergyComposite->SetEnergyValue(fEnergyGenerator);
+
+ if ( fGenPositionComposite == nullptr ) fGenPositionComposite = new Kassiopeia::KSGenPositionRectangularComposite();
+ fGenPositionComposite->SetOrigin(GetKGWorldSpace()->GetOrigin());
+ if ( fPositionXGenerator == nullptr ) fPositionXGenerator = new Kassiopeia::KSGenValueUniform();
+ if ( fPositionYGenerator == nullptr ) fPositionYGenerator = new Kassiopeia::KSGenValueUniform();
+ if ( fPositionZGenerator == nullptr ) fPositionZGenerator = new Kassiopeia::KSGenValueUniform();
+ fPositionXGenerator->SetValueMin( aParam["ks-starting-xpos-min"]().as_double() ); // meters
+ fPositionXGenerator->SetValueMax( aParam["ks-starting-xpos-max"]().as_double() );
+ fPositionYGenerator->SetValueMin( aParam["ks-starting-ypos-min"]().as_double() );
+ fPositionYGenerator->SetValueMax( aParam["ks-starting-ypos-max"]().as_double() );
+ fPositionZGenerator->SetValueMin( aParam["ks-starting-zpos-min"]().as_double() );
+ fPositionZGenerator->SetValueMax( aParam["ks-starting-zpos-max"]().as_double() );
+ fGenPositionComposite->SetXValue(fPositionXGenerator);
+ fGenPositionComposite->SetYValue(fPositionYGenerator);
+ fGenPositionComposite->SetZValue(fPositionZGenerator);
+
+ if ( fGenDirectionComposite == nullptr ) fGenDirectionComposite = new Kassiopeia::KSGenDirectionSphericalComposite();
+ if ( fThetaGenerator == nullptr ) fThetaGenerator = new Kassiopeia::KSGenValueAngleSpherical();
+ fThetaGenerator->SetAngleMin( aParam["ks-starting-pitch-min"]().as_double() );
+ fThetaGenerator->SetAngleMax( aParam["ks-starting-pitch-max"]().as_double() );
+ if ( fPhiGenerator == nullptr ) fPhiGenerator = new Kassiopeia::KSGenValueUniform();
+ if ( aParam.has( "ks-starting-phi-min" ) && ( aParam.has( "ks-starting-phi-max" ) ) )
+ {
+ fPhiGenerator->SetValueMin( aParam["ks-starting-phi-min"]().as_double() );
+ fPhiGenerator->SetValueMax( aParam["ks-starting-phi-max"]().as_double() );
+ }
+ else
+ {
+ fPhiGenerator->SetValueMin( 0. );
+ fPhiGenerator->SetValueMax( 0. );
+ }
+ fGenDirectionComposite->SetPhiValue(fPhiGenerator);
+ fGenDirectionComposite->SetThetaValue(fThetaGenerator);
+
+ if ( fGenerator == nullptr ) fGenerator = new Kassiopeia::KSGenGeneratorComposite();
+
+ fGenerator->SetPid(fGenPidComposite);
+ fGenerator->AddCreator(fGenPositionComposite);
+ fGenerator->AddCreator(fGenEnergyComposite);
+ fGenerator->AddCreator(fGenDirectionComposite);
+ fGenerator->AddCreator(fGenTimeComposite);
+ fGenerator->SetName("gen_project8");
+ fGenerator->Initialize();
+ fGenerator->Activate();
- Kassiopeia::KSGenEnergyComposite* tGenEnergyComposite = new Kassiopeia::KSGenEnergyComposite();
- Kassiopeia::KSGenValueUniform* tEnergyGenerator = new Kassiopeia::KSGenValueUniform();
- if ( aParam.has( "ks-starting-energy-min" ) && ( aParam.has( "ks-starting-energy-max" ) ) )
- {
- tEnergyGenerator->SetValueMin( aParam["ks-starting-energy-min"]().as_double() ); // eV
- tEnergyGenerator->SetValueMax( aParam["ks-starting-energy-max"]().as_double() ); // eV
- }
- else
- {
- tEnergyGenerator->SetValueMin( 18600. ); // eV
- tEnergyGenerator->SetValueMax( 18600. ); // eV
- }
- tGenEnergyComposite->SetEnergyValue(tEnergyGenerator);
-
- Kassiopeia::KSGenPositionRectangularComposite* tGenPositionComposite = new Kassiopeia::KSGenPositionRectangularComposite();
- tGenPositionComposite->SetOrigin(GetKGWorldSpace()->GetOrigin());
- Kassiopeia::KSGenValueUniform* tPositionXGenerator = new Kassiopeia::KSGenValueUniform();
- Kassiopeia::KSGenValueUniform* tPositionYGenerator = new Kassiopeia::KSGenValueUniform();
- Kassiopeia::KSGenValueUniform* tPositionZGenerator = new Kassiopeia::KSGenValueUniform();
- tPositionXGenerator->SetValueMin( aParam["ks-starting-xpos-min"]().as_double() ); // meters
- tPositionXGenerator->SetValueMax( aParam["ks-starting-xpos-max"]().as_double() );
- tPositionYGenerator->SetValueMin( aParam["ks-starting-ypos-min"]().as_double() );
- tPositionYGenerator->SetValueMax( aParam["ks-starting-ypos-max"]().as_double() );
- tPositionZGenerator->SetValueMin( aParam["ks-starting-zpos-min"]().as_double() );
- tPositionZGenerator->SetValueMax( aParam["ks-starting-zpos-max"]().as_double() );
- tGenPositionComposite->SetXValue(tPositionXGenerator);
- tGenPositionComposite->SetYValue(tPositionYGenerator);
- tGenPositionComposite->SetZValue(tPositionZGenerator);
-
- Kassiopeia::KSGenDirectionSphericalComposite* tGenDirectionComposite = new Kassiopeia::KSGenDirectionSphericalComposite();
- Kassiopeia::KSGenValueUniform* tThetaGenerator = new Kassiopeia::KSGenValueUniform();
- tThetaGenerator->SetValueMin( aParam["ks-starting-pitch-min"]().as_double() );
- tThetaGenerator->SetValueMax( aParam["ks-starting-pitch-max"]().as_double() );
- Kassiopeia::KSGenValueUniform* tPhiGenerator = new Kassiopeia::KSGenValueUniform();
- if ( aParam.has( "ks-starting-phi-min" ) && ( aParam.has( "ks-starting-phi-max" ) ) )
- {
- tPhiGenerator->SetValueMin( aParam["ks-starting-phi-min"]().as_double() );
- tPhiGenerator->SetValueMax( aParam["ks-starting-phi-max"]().as_double() );
- }
- else
- {
- tPhiGenerator->SetValueMin( 0. );
- tPhiGenerator->SetValueMax( 0. );
- }
- tGenDirectionComposite->SetPhiValue(tPhiGenerator);
- tGenDirectionComposite->SetThetaValue(tThetaGenerator);
-
- fGenerator = new Kassiopeia::KSGenGeneratorComposite();
- fGenerator->SetPid(tGenPidComposite);
- fGenerator->AddCreator(tGenPositionComposite);
- fGenerator->AddCreator(tGenEnergyComposite);
- fGenerator->AddCreator(tGenDirectionComposite);
- fGenerator->AddCreator(tGenTimeComposite);
- fGenerator->SetName("gen_project8");
- fGenerator->Initialize();
- fGenerator->Activate();
-
- if (!fToolbox.HasKey("gen_project8"))
- {
fToolbox.Add(fGenerator);
fToolbox.Get("root_generator")->SetGenerator(fGenerator);
+ LPROG(lmclog,"\"gen-project8\" has just been added to the KToolbox.");
}
+ else
+ {
+ LERROR(lmclog,"To configure starting e- kinematics, all of these parameters are needed: "
+ "ks-starting-xpos-min, ks-starting-xpos-max, ks-starting-ypos-min, ks-starting-ypos-max, "
+ "ks-starting-zpos-min, ks-starting-zpos-max, ks-starting-pitch-min, ks-starting-pitch-max,"
+ "ks-starting-energy-min, ks-starting-energy-max ");
+ return false;
+ }
}
-
else
{
- LERROR(lmclog,"To configure starting e- kinematics, all of these parameters are needed: "
- "ks-starting-xpos-min, ks-starting-xpos-max, ks-starting-ypos-min, ks-starting-ypos-max, "
- "ks-starting-zpos-min, ks-starting-zpos-max, ks-starting-pitch-min, ks-starting-pitch-max,"
- "ks-starting-energy-min, ks-starting-energy-max ");
- return false;
+ LPROG(lmclog,"\"gen-project8\" is already in the KToolbox.");
}
return true;
@@ -228,22 +291,28 @@ namespace locust
bool RunPause::AddMaxRTerminator( const scarab::param_node& aParam )
{
- /* Remove any existing KSTermMaxR objects */
- auto tMaxR = fToolbox.GetAll();
- for (unsigned i=0; i("root_terminator")->RemoveTerminator(tMaxR[i]);
- fToolbox.Remove(tMaxR[i]->GetName());
- }
-
if (!fToolbox.HasKey("ksmax-r-project8"))
{
- fLocustMaxRTerminator = new Kassiopeia::KSTermMaxR();
+ /* Remove any existing KSTermMaxR objects */
+ auto tMaxR = fToolbox.GetAll();
+ for (unsigned i=0; i("root_terminator")->RemoveTerminator(tMaxR[i]);
+ fToolbox.Remove(tMaxR[i]->GetName());
+ }
+
+ if ( fLocustMaxRTerminator == nullptr ) fLocustMaxRTerminator = new Kassiopeia::KSTermMaxR();
fLocustMaxRTerminator->SetName("ksmax-r-project8");
fLocustMaxRTerminator->SetMaxR( aParam["cavity-radius"]().as_double() );
fLocustMaxRTerminator->Initialize();
+ fLocustMaxRTerminator->Activate();
fToolbox.Add(fLocustMaxRTerminator);
fToolbox.Get("root_terminator")->AddTerminator(fLocustMaxRTerminator);
+ LPROG(lmclog,"\"ksmax-r-project8\" has just been added to the KToolbox.");
+ }
+ else
+ {
+ LPROG(lmclog,"\"ksmax-r-project8\" is already in the KToolbox.");
}
return true;
}
@@ -267,40 +336,39 @@ namespace locust
bool RunPause::AddMaxTimeTerminator( const scarab::param_node& aParam )
{
-
- /* Remove any existing KSTermMaxTime objects */
- auto tMaxTime = fToolbox.GetAll();
- for (unsigned i=0; i("root_terminator")->RemoveTerminator(tMaxTime[i]);
- fToolbox.Remove(tMaxTime[i]->GetName());
- }
-
if (!fToolbox.HasKey("ksmax-time-project8"))
{
+ /* Remove any existing KSTermMaxTime objects */
+ auto tMaxTime = fToolbox.GetAll();
+ for (unsigned i=0; i("root_terminator")->RemoveTerminator(tMaxTime[i]);
+ fToolbox.Remove(tMaxTime[i]->GetName());
+ }
+
double tMaxTrackLength = 0.;
- fLocustMaxTimeTerminator = new Kassiopeia::KSTermMaxTime();
+ if ( fLocustMaxTimeTerminator == nullptr ) fLocustMaxTimeTerminator = new Kassiopeia::KSTermMaxTime();
if ( aParam.has( "min-track-length-fraction" ) )
- {
+ {
fMinTrackLengthFraction = aParam["min-track-length-fraction"]().as_double();
LPROG(lmclog,"Setting minimum track length fraction to " << fMinTrackLengthFraction);
}
if ( aParam.has( "track-length" ) )
- {
- tMaxTrackLength = aParam["track-length"]().as_double();
- fLocustMaxTimeTerminator->SetTime( tMaxTrackLength );
- }
- else
- {
+ {
+ tMaxTrackLength = aParam["track-length"]().as_double();
+ fLocustMaxTimeTerminator->SetTime( tMaxTrackLength );
+ }
+ else
+ {
LERROR(lmclog,"Parameter \"track-length\" is required to set the maximum track length.");
return false;
- }
+ }
- if ( aParam.has( "random-track-length" ) )
- {
+ if ( aParam.has( "random-track-length" ) )
+ {
if ( aParam["random-track-length"]().as_bool() == true)
{
scarab::param_node default_setting;
@@ -312,57 +380,63 @@ namespace locust
fLocustMaxTimeTerminator->SetTime( tRandomTime );
LPROG(lmclog,"Randomizing the track length to " << tRandomTime);
}
- }
+ }
fLocustMaxTimeTerminator->SetName("ksmax-time-project8");
fLocustMaxTimeTerminator->Initialize();
fLocustMaxTimeTerminator->Activate();
fToolbox.Add(fLocustMaxTimeTerminator);
fToolbox.Get("root_terminator")->AddTerminator(fLocustMaxTimeTerminator);
+ LPROG(lmclog,"\"ksmax-time-project8\" has just been added to the KToolbox.");
+ }
+ else
+ {
+ LPROG(lmclog,"\"ksmax-time-project8\" is already in the KToolbox.");
}
return true;
}
bool RunPause::AddWaveguideTerminator( const scarab::param_node& aParam )
{
- if ( aParam.has( "waveguide-y" ) && aParam.has( "waveguide-z" ) )
+ if (!fToolbox.HasKey("waveguide_surfaces_project8"))
{
- fBox = new KGeoBag::KGBoxSpace();
- fBox->XA(-aParam["waveguide-x"]().as_double()/2.);
- fBox->XB(aParam["waveguide-x"]().as_double()/2.);
- fBox->YA(-aParam["waveguide-y"]().as_double()/2.);
- fBox->YB(aParam["waveguide-y"]().as_double()/2.);
- fBox->ZA(-aParam["waveguide-z"]().as_double()/2.);
- fBox->ZB(aParam["waveguide-z"]().as_double()/2.);
- fBox->SetTag("waveguide_box");
-
- fKGSpace = GetKGWorldSpace();
- KGeoBag::KGSpace* tKGSpace = new KGeoBag::KGSpace();
- tKGSpace->Volume(std::shared_ptr(fBox));
- fKGSpace->GetChildSpaces()->at(0)->AddChildSpace(tKGSpace);
-
- fSurface = new Kassiopeia::KSGeoSurface();
- fSurface->SetName("waveguide_surfaces_project8");
- auto it = tKGSpace->GetBoundaries()->begin();
- while (it != tKGSpace->GetBoundaries()->end())
- {
- fSurface->AddContent(*it);
- *it++;
- }
- if (!fToolbox.HasKey("waveguide_surfaces_project8"))
+ if ( aParam.has( "waveguide-y" ) && aParam.has( "waveguide-z" ) )
{
+ if ( fBox == nullptr ) fBox = new KGeoBag::KGBoxSpace();
+ fBox->XA(-aParam["waveguide-x"]().as_double()/2.);
+ fBox->XB(aParam["waveguide-x"]().as_double()/2.);
+ fBox->YA(-aParam["waveguide-y"]().as_double()/2.);
+ fBox->YB(aParam["waveguide-y"]().as_double()/2.);
+ fBox->ZA(-aParam["waveguide-z"]().as_double()/2.);
+ fBox->ZB(aParam["waveguide-z"]().as_double()/2.);
+ fBox->SetTag("waveguide_box");
+
+ if ( fKGSpace == nullptr ) fKGSpace = GetKGWorldSpace();
+ KGeoBag::KGSpace* tKGSpace = new KGeoBag::KGSpace();
+ tKGSpace->Volume(std::shared_ptr(fBox));
+ fKGSpace->GetChildSpaces()->at(0)->AddChildSpace(tKGSpace);
+
+ if ( fSurface == nullptr ) fSurface = new Kassiopeia::KSGeoSurface();
+ fSurface->SetName("waveguide_surfaces_project8");
+ auto it = tKGSpace->GetBoundaries()->begin();
+ while (it != tKGSpace->GetBoundaries()->end())
+ {
+ fSurface->AddContent(*it);
+ *it++;
+ }
+
fToolbox.Add(fSurface);
fLocustTermDeath = new Kassiopeia::KSTermDeath();
fLocustTermDeath->Initialize();
+ fLocustTermDeath->Activate();
fToolbox.Add(fLocustTermDeath);
fLocustTermDeath->SetName("waveguide_death");
fCommand = fToolbox.Get("root_terminator")->Command("add_terminator", fLocustTermDeath);
- fKSSpace = GetKSWorldSpace();
+ if ( fKSSpace == nullptr ) fKSSpace = GetKSWorldSpace();
fKSSpace->AddSurface(fSurface);
fSurface->AddCommand(fCommand);
}
- return true;
}
else
{
@@ -370,6 +444,8 @@ namespace locust
return false;
}
+ return true;
+
}
@@ -426,11 +502,6 @@ namespace locust
return true;
}
}
-//.
-//.
-// Between these two functions there will be an nevents interrupt in KSRoot.
-//.
-//.
bool RunPause::DeleteLocalKassObjects()
{
@@ -457,8 +528,6 @@ namespace locust
fLocustTermDeath = NULL;
fCommand = NULL;
fBox = NULL;
- fKGSpace = NULL;
- fKSSpace = NULL;
fSurface = NULL;
LPROG(lmclog,"Removing fSurface from KSRoot ...");
}
@@ -466,10 +535,21 @@ namespace locust
// fGenerator
if (fToolbox.HasKey("gen_project8"))
{
- fKGSpace = NULL;
- fKSSpace = NULL;
+ fToolbox.Get("root_generator")->ClearGenerator(fGenerator);
+ fToolbox.Remove(fGenerator->GetName());
fGenerator = NULL;
- LPROG(lmclog,"Leaving fGenerator in KSRoot, but setting LMC pointers to Kass objects to NULL ...");
+ LPROG(lmclog,"Removing fGenerator from KToolbox ... ");
+
+ auto tGen = fToolbox.GetAll();
+ for (unsigned i=0; iIsActivated()) && (tGen[i]->GetName()!="root_generator") )
+ {
+ LPROG(lmclog,"Replacing KSRoot generator with " << tGen[i]->GetName());
+ fToolbox.Get("root_generator")->SetGenerator(tGen[i]);
+ }
+ }
+
}
return true;
@@ -477,12 +557,7 @@ namespace locust
bool RunPause::ExecutePostRunModification(Kassiopeia::KSRun & aRun)
{
- // No interrupt has happened yet in KSRoot. Run still in progress.
-
- // Do not delete anything, local or otherwise, that is used by the KToolbox, in particular
- // if running multiple events or multiple tracks. Keep this next line commented out for now.
-// DeleteLocalKassObjects();
-
+ // No interrupt has happened yet in KSRoot. Run still in progress.
return true;
}
diff --git a/Source/Kassiopeia/LMCRunPause.hh b/Source/Kassiopeia/LMCRunPause.hh
index aa36998f..bb0beb65 100644
--- a/Source/Kassiopeia/LMCRunPause.hh
+++ b/Source/Kassiopeia/LMCRunPause.hh
@@ -21,6 +21,7 @@
#include "KSRootGenerator.h"
#include "KSGenGeneratorComposite.h"
#include "KSGenValueUniform.h"
+#include "KSGenValueAngleSpherical.h"
#include "KSGenValueFix.h"
#include "KSGenEnergyComposite.h"
#include "KSGenPositionRectangularComposite.h"
@@ -79,11 +80,25 @@ namespace locust
Kassiopeia::KSCommand* fCommand;
Kassiopeia::KSGeoSpace* fKSSpace;
Kassiopeia::KSGenGeneratorComposite* fGenerator;
+ Kassiopeia::KSGenDirectionSphericalComposite* fGenDirectionComposite;
+ Kassiopeia::KSGenValueAngleSpherical* fThetaGenerator;
+ Kassiopeia::KSGenValueUniform* fPhiGenerator;
+ Kassiopeia::KSGenPositionRectangularComposite* fGenPositionComposite;
+ Kassiopeia::KSGenValueUniform* fPositionXGenerator;
+ Kassiopeia::KSGenValueUniform* fPositionYGenerator;
+ Kassiopeia::KSGenValueUniform* fPositionZGenerator;
+ Kassiopeia::KSGenValueUniform* fEnergyGenerator;
+ Kassiopeia::KSGenEnergyComposite* fGenEnergyComposite;
+ Kassiopeia::KSGenTimeComposite* fGenTimeComposite;
+ Kassiopeia::KSGenValueUniform* fTimeGenerator;
+ Kassiopeia::KSGenValueFix* fGenPidComposite;
std::shared_ptr< BaseDistribution> fTrackLengthDistribution;
DistributionInterface fDistributionInterface;
double fMinTrackLengthFraction;
bool fConfigurationComplete;
+ int fEventCounter;
+ int fMaxEvents;
diff --git a/Source/Kassiopeia/LMCTrackHold.cc b/Source/Kassiopeia/LMCTrackHold.cc
index db166902..09856b63 100644
--- a/Source/Kassiopeia/LMCTrackHold.cc
+++ b/Source/Kassiopeia/LMCTrackHold.cc
@@ -64,6 +64,7 @@ namespace locust
bool TrackHold::ExecutePreTrackModification(Kassiopeia::KSTrack &aTrack)
{
fInterface->aTrack.Initialize();
+ fInterface->fNewTrackStarting = true;
double tTime = aTrack.GetInitialParticle().GetTime();
double tPitchAngle = aTrack.GetInitialParticle().GetPolarAngleToB();
@@ -73,11 +74,13 @@ namespace locust
bool TrackHold::ExecutePostTrackModification(Kassiopeia::KSTrack &aTrack)
{
- fInterface->anEvent->AddTrack( fInterface->aTrack );
+ if ( aTrack.GetTotalSteps() > 0)
+ {
+ fInterface->anEvent->AddTrack( fInterface->aTrack );
+ }
- double tPitchAngle = aTrack.GetFinalParticle().GetPolarAngleToB();
double tTime = aTrack.GetFinalParticle().GetTime();
- LWARN(lmclog,"LMCTrack " << fTrackCounter << " is complete at Kass time " << tTime << ", with instantaneous pitch angle " << tPitchAngle);
+ LWARN(lmclog,"LMCTrack " << fTrackCounter << " is complete at Kass time " << tTime << " with total steps " << aTrack.GetTotalSteps() );
fTrackCounter += 1;
return true;
}