Added tests for crosssection and splitting biasing

01_element/test_drift.py

@@ -33,8 +33,6 @@ def test() :
                 [0,0,0,1,0,0],
                 [0,0,0,0,1,0],
                 [0,0,0,0,0,1]]
-    print('rounded matrix',pybdsim.Testing.round_matrix(rmat))
-    print('maximum matrix difference',pybdsim.Testing.max_matrix_diff(rmat,ref_rmat))    
 
     assert pybdsim.Testing.compare_matrix(rmat,ref_rmat)
 

01_element/test_quadrupole.py

@@ -46,12 +46,9 @@ def test() :
                     [                    0,                      0, sqrt(abs(k1))*sinh(a),                 cosh(a),0,0],
                     [                    0,                      0,                     0,                       0,1,0],
                     [                    0,                      0,                     0,                       0,0,0]]
-    print('rounded matrix')
     print(pybdsim.Testing.round_matrix(rmat,3))
-    print('rounded reference matrix')
     print(pybdsim.Testing.round_matrix(ref_rmat,3))
 
-    print('maximum matrix difference',pybdsim.Testing.max_matrix_diff(rmat,ref_rmat))    
 
     assert pybdsim.Testing.compare_matrix(rmat,ref_rmat)
 

01_element/test_rbend.py

@@ -40,11 +40,8 @@ def test() :
                 [ 0.,             0.,                0.,             0.,   1.,   0.],
                 [ 0.,             0.,                0.,             0.,  -0.,   0.]]
 
-    print('rounded bdsim matrix')
     print(pybdsim.Testing.round_matrix(rmat,3))
-    print('rounded reference matrix')
     print(pybdsim.Testing.round_matrix(ref_rmat,3))
-    print('maximum matrix difference')
     print(pybdsim.Testing.max_matrix_diff(rmat,ref_rmat))
 
     assert pybdsim.Testing.compare_matrix(rmat,ref_rmat,2e-3)

01_element/test_sbend.py

@@ -34,11 +34,8 @@ def test() :
                 [               0.,              0.,   0.,      1.,   0.,   0. ],
                 [               0.,              0.,   0.,      0.,   1.,   0. ],
                 [               0.,              0.,   0.,      0.,   0.,   0. ]]
-    print('rounded matrix')
     print(pybdsim.Testing.round_matrix(rmat,3))
-    print('rounded reference matrix')
     print(pybdsim.Testing.round_matrix(ref_rmat,3))
-    print('maximum matrix difference')
     print(pybdsim.Testing.max_matrix_diff(rmat,ref_rmat))
 
     assert pybdsim.Testing.compare_matrix(rmat,ref_rmat)

01_element/test_sextupole.py

@@ -28,7 +28,7 @@ def test() :
 
     do = pybdsim.DataPandas.REBDSIMOptics(optics_name)
     do_df = do.get_optics()
-    print(do_df)
+
 
     rmat = pybdsim.Analysis.CalculateTaylorMapOrder2(root_name,"d1.","s1.", average=True)
 
@@ -38,8 +38,7 @@ def test() :
     #            [0,0,0,1,0,0],
     #            [0,0,0,0,1,0],
     #            [0,0,0,0,0,1]]
-    print('rounded matrix')
-    print(pybdsim.Testing.round_matrix(rmat,3))
+
 
     #print('maximum matrix difference',pybdsim.Testing.max_matrix_diff(rmat,ref_rmat))
     #assert pybdsim.Testing.compare_matrix(rmat,ref_rmat)

04_options/crosssection_biasing.tpl

@@ -0,0 +1,22 @@
+d1: drift, l=1*m, aper1=10*cm;
+AlBlock: rcol, l=2*m, material="G4_Al";
+
+l0 : line = (d1,AlBlock);
+
+use, period=l0;
+
+sample, all;
+
+option, physicsList = "g4FTFP_BERT";
+
+piPlusBias:   xsecBias, particle="pi+", proc="Decay", xsecfact={{ BIAS_FACTOR }}, flag=1;
+piMinusBias:  xsecBias, particle="pi-", proc="Decay", xsecfact={{ BIAS_FACTOR }}, flag=1;
+
+option, biasForWorldVolume="piPlusBias piMinusBias";
+option, biasForWorldContents="piPlusBias piMinusBias";
+option, biasForWorldVacuum = "piPlusBias piMinusBias";
+
+beam, particle="pi+",
+      kineticEnergy=10*GeV;
+
+sampler: samplerplacement, z=3*m, aper1=5*m, aper2=5*m, shape="rectangular", partID={13,-13};

04_options/muon_splitting.tpl

@@ -0,0 +1,17 @@
+d1: drift, l=1*m, aper1=10*cm;
+AlBlock: rcol, l=2*m, material="G4_Al";
+
+l0 : line = (d1,AlBlock);
+
+use, period=l0;
+
+sample, all;
+
+option, physicsList = "g4FTFP_BERT";
+
+option, muonSplittingFactor={{ SPLITTING_FACTOR }};
+
+beam, particle="pi+",
+      kineticEnergy=10*GeV;
+
+sampler: samplerplacement, z=3*m, aper1=5*m, aper2=5*m, shape="rectangular", partID={13,-13};

04_options/test_crosssection_biasing.py

@@ -0,0 +1,26 @@
+import pytest
+import pybdsim
+import os
+
+def test() :
+
+    os.chdir(os.path.dirname(__file__))
+
+    base_name     = "crosssection_biasing"
+    template_name = base_name+".tpl"
+    gmad_name     = base_name+".gmad"
+    root_name     = base_name+".root"
+
+    data = {'BIAS_FACTOR': '10'}
+
+    pybdsim.Run.RenderGmadJinjaTemplate(template_name,gmad_name,data)
+    pybdsim.Run.Bdsim(gmad_name,base_name,ngenerate=100,seed=1)
+
+    d = pybdsim.DataPandas.BDSIMOutput(root_name)
+    s = d.get_sampler("sampler.")
+    n_muons = len(s['x'])
+    weight_sum = sum(s['weight'])
+
+    # count varies with G4 version but unbiased muon count is ~12
+    assert(20 < n_muons < 50)
+    assert(6 < weight_sum < 18)

04_options/test_muon_splitting.py

@@ -0,0 +1,26 @@
+import pytest
+import pybdsim
+import os
+
+def test() :
+
+    os.chdir(os.path.dirname(__file__))
+
+    base_name     = "muon_splitting"
+    template_name = base_name+".tpl"
+    gmad_name     = base_name+".gmad"
+    root_name     = base_name+".root"
+
+    data = {'SPLITTING_FACTOR': '30'}
+
+    pybdsim.Run.RenderGmadJinjaTemplate(template_name,gmad_name,data)
+    pybdsim.Run.Bdsim(gmad_name,base_name,ngenerate=100,seed=1)
+
+    d = pybdsim.DataPandas.BDSIMOutput(root_name)
+    s = d.get_sampler("sampler.")
+    n_muons = len(s['x'])
+    weight_sum = sum(s['weight'])
+
+    # count varies with G4 version but unbiased muon count is ~12
+    assert(160 < n_muons < 500)
+    assert(6 < weight_sum < 18)

08_processes/laserwire_compton_cumulative.tpl

@@ -0,0 +1,27 @@
+
+lb1: laserflux, wavelength=532.0*nm,
+m2=1.2,
+pulseDuration=71.0*ps,
+pulseEnergy={{ LASER_ENERGY }}*mJ,
+w0=2.3*um;
+
+laser1: laserwire, l=0.01*m, wireLength=0.1*m, laserOffsetTheta=-pi/2, laserOffsetPhi=0,laserOffsetY=0.0*m, laserOffsetZ=0*m, laserOffsetX=0, laserBeam="lb1";
+
+lat : line=(laser1);
+use, lat;
+
+beam, particle="e-",
+      energy={{ BEAM_ENERGY }}*GeV,
+      X0=0.0*m,
+      Xp0=0.0,
+      Y0=0.0*m,
+      Yp0=0.0, 
+      distrType="gauss",
+      sigmaX=119*um,
+      sigmaY=1.07*um,
+      sigmaE=0,
+      sigmaT=0;
+option, 
+physicsList="laser_cumulative_compton_scattering";
+
+sample, range=laser1;

08_processes/laserwire_compton_multiStep.tpl

@@ -0,0 +1,27 @@
+
+lb1: laserflux, wavelength=532.0*nm,
+m2=1.2,
+pulseDuration=71.0*ps,
+pulseEnergy={{ LASER_ENERGY }}*mJ,
+w0=2.3*um;
+
+laser1: laserwire, l=0.01*m, wireLength=0.1*m, laserOffsetTheta=-pi/2, laserOffsetPhi=0,laserOffsetY=0.0*m, laserOffsetZ=0*m, laserOffsetX=0, laserBeam="lb1";
+
+lat : line=(laser1);
+use, lat;
+
+beam, particle="e-",
+      energy={{ BEAM_ENERGY }}*GeV,
+      X0=0.0*m,
+      Xp0=0.0,
+      Y0=0.0*m,
+      Yp0=0.0, 
+      distrType="gauss",
+      sigmaX=119*um,
+      sigmaY=1.07*um,
+      sigmaE=0,
+      sigmaT=0;
+option, 
+physicsList="laser_compton_scattering", scaleFactorLaser=1000;
+
+sample, range=laser1;

08_processes/laserwire_photodetachment_cumulative.tpl

@@ -0,0 +1,26 @@
+lb1: laserflux, wavelength=1064.0*nm,
+m2=1.2,
+pulseDuration=106.0*ns,
+pulseEnergy={{ LASER_ENERGY }}*mJ,
+w0=50*mm;
+
+laser1: laserwire, l=0.01*m, wireLength=0.1*m, laserOffsetTheta=-pi/2, laserOffsetPhi=0,laserOffsetY=0.0*m, laserOffsetZ=0*m, laserOffsetX=0, laserBeam="lb1";
+
+lat : line=(laser1);
+use, lat;
+
+beam, particle="ion 1 1 -1",
+      energy={{ BEAM_ENERGY }}*MeV,
+      X0=0.0*m,
+      Xp0=0.0,
+      Y0=0.0*m,
+      Yp0=0.0, 
+      distrType="gauss",
+      sigmaX=2*mm,
+      sigmaY=4*mm,
+      sigmaE=0,
+      sigmaT=0;
+option, 
+physicsList="laser_cumulative_photo_detachment";
+
+sample, range=laser1;

08_processes/laserwire_photodetachment_multiStep.tpl

@@ -0,0 +1,26 @@
+lb1: laserflux, wavelength=1064.0*nm,
+m2=1.2,
+pulseDuration=106.0*ns,
+pulseEnergy={{ LASER_ENERGY }}*mJ,
+w0=50*mm;
+
+laser1: laserwire, l=0.01*m, wireLength=0.1*m, laserOffsetTheta=-pi/2, laserOffsetPhi=0,laserOffsetY=0.0*m, laserOffsetZ=0*m, laserOffsetX=0, laserBeam="lb1";
+
+lat : line=(laser1);
+use, lat;
+
+beam, particle="ion 1 1 -1",
+      energy={{ BEAM_ENERGY }}*MeV,
+      X0=0.0*m,
+      Xp0=0.0,
+      Y0=0.0*m,
+      Yp0=0.0, 
+      distrType="gauss",
+      sigmaX=2*mm,
+      sigmaY=4*mm,
+      sigmaE=0,
+      sigmaT=0;
+option, 
+physicsList="laser_photo_detachment", scaleFactorLaser=1000;
+
+sample, range=laser1;

08_processes/test_laserwire_compton_cumulative.py

@@ -0,0 +1,37 @@
+import pytest
+import pybdsim
+import os
+
+def test() :
+
+    os.chdir(os.path.dirname(__file__))
+
+    base_name     = "laserwire_compton_cumulative"
+    template_name = base_name+".tpl"
+    gmad_name     = base_name+".gmad"
+    root_name     = base_name+".root"
+
+    data = {
+        'LASER_ENERGY': '150',
+        'BEAM_ENERGY' : '1.3'
+    }
+
+    pybdsim.Run.RenderGmadJinjaTemplate(template_name,gmad_name,data)
+    pybdsim.Run.Bdsim(gmad_name,base_name,10000,1)
+
+    d = pybdsim.Data.Load(root_name)
+    samplerData = pybdsim.Data.SamplerData(d,'laser1') 
+    weights=samplerData.data.get("weight")
+    partid=samplerData.data.get("partID")
+    wpho=[]
+
+    for i in range(len(partid)):
+        if partid[i]==22:
+            wpho.append(weights[i])
+
+
+    Npho=sum(wpho)
+    ref_Npho=0.626774271968543
+
+
+    assert (Npho==pytest.approx(ref_Npho,abs=1e-3))

08_processes/test_laserwire_compton_multiStep.py

@@ -0,0 +1,33 @@
+import pytest
+import pybdsim
+import os
+
+def test() :
+
+    os.chdir(os.path.dirname(__file__))
+
+    base_name     = "laserwire_compton_multiStep"
+    template_name = base_name+".tpl"
+    gmad_name     = base_name+".gmad"
+    root_name     = base_name+".root"
+
+    l  = 2.0 
+    data = {
+        'LASER_ENERGY': '150',
+        'BEAM_ENERGY' : '1.3'
+    }
+
+    pybdsim.Run.RenderGmadJinjaTemplate(template_name,gmad_name,data)
+    pybdsim.Run.Bdsim(gmad_name,base_name,10000,1)
+
+    d = pybdsim.Data.Load(root_name)
+    samplerData = pybdsim.Data.SamplerData(d,'laser1') 
+    partid=samplerData.data.get('partID')
+    Npho=0
+    for i in range(len(partid)):
+        if partid[i]==22:
+            Npho+=1
+
+    ref_Npho=637
+
+    assert (Npho==ref_Npho)

08_processes/test_laserwire_photodetachment_cumulative.py

@@ -0,0 +1,35 @@
+import pytest
+import pybdsim
+import os
+
+def test() :
+
+    os.chdir(os.path.dirname(__file__))
+
+    base_name     = "laserwire_photodetachment_cumulative"
+    template_name = base_name+".tpl"
+    gmad_name     = base_name+".gmad"
+    root_name     = base_name+".root"
+
+    data = {
+        'LASER_ENERGY': '67.4',
+        'BEAM_ENERGY' : '0.942022+0.003'
+    }
+
+    pybdsim.Run.RenderGmadJinjaTemplate(template_name,gmad_name,data)
+    pybdsim.Run.Bdsim(gmad_name,base_name,10000,1)
+
+    d = pybdsim.Data.Load(root_name)
+    samplerData = pybdsim.Data.SamplerData(d,'laser1') 
+    weights=samplerData.data.get("weight")
+    partid=samplerData.data.get("partID")
+
+    we=[]
+    for i in range(len(weights)):
+        if partid[i]==11:
+            we.append(weights[i])
+    Ne=sum(we)
+    ref_Ne=5.843840659450507
+
+
+    assert (Ne==ref_Ne)