From a0c32c3ed8a8926eb1902e590407ca1f0561b2ff Mon Sep 17 00:00:00 2001
From: Oliver Kirsebom <oliver.kirsebom@gmail.com>
Date: Fri, 2 Mar 2018 15:00:04 -0400
Subject: [PATCH] possibility to set coefficients of shape function

---
 include/VKPrimaryGeneratorAction.hh    |  6 +++
 include/VKPrimaryGeneratorMessenger.hh |  2 +
 mac/b12scan.mac                        |  2 +
 mac/f20scan.mac                        | 24 ++----------
 mac/forbidden.mac                      | 25 ++++++------
 mac/init.mac                           |  4 +-
 mac/test.mac                           | 53 +++++++++-----------------
 src/VKPrimaryGeneratorAction.cpp       | 15 +++++++-
 src/VKPrimaryGeneratorMessenger.cpp    | 22 +++++++++++
 9 files changed, 83 insertions(+), 70 deletions(-)

diff --git a/include/VKPrimaryGeneratorAction.hh b/include/VKPrimaryGeneratorAction.hh
index 64ff8ae..086120a 100644
--- a/include/VKPrimaryGeneratorAction.hh
+++ b/include/VKPrimaryGeneratorAction.hh
@@ -80,6 +80,10 @@ class VKPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     void SetBetaShapeFuncA5(double a5);
     void SetBetaShapeFuncA6(double a6);
     
+    // Set Z and A for calculation of corrections of beta spectrum
+    void SetBetaZ(int Z);
+    void SetBetaA(int A);
+    
     // Get coverage of sampling
     double GetCoverage() const;
     
@@ -134,6 +138,8 @@ class VKPrimaryGeneratorAction : public G4VUserPrimaryGeneratorAction
     TF1 * shapeFunc;
     G4bool fBetaSpectrumOn;
     double shapeFuncCoeff[6];
+    int fBetaZ;
+    int fBetaA;
     
     G4bool fDiskShapedSourceOn;
     G4double fSourceInnerDiameter;
diff --git a/include/VKPrimaryGeneratorMessenger.hh b/include/VKPrimaryGeneratorMessenger.hh
index 0cc2640..5ced949 100644
--- a/include/VKPrimaryGeneratorMessenger.hh
+++ b/include/VKPrimaryGeneratorMessenger.hh
@@ -70,6 +70,8 @@ class VKPrimaryGeneratorMessenger: public G4UImessenger
     G4UIcmdWithADouble *fShapeFuncA4Cmd;
     G4UIcmdWithADouble *fShapeFuncA5Cmd;
     G4UIcmdWithADouble *fShapeFuncA6Cmd;
+    G4UIcmdWithAnInteger *fBetaZCmd;
+    G4UIcmdWithAnInteger *fBetaACmd;
 };
 
 #endif
diff --git a/mac/b12scan.mac b/mac/b12scan.mac
index ff89a49..e48a43b 100644
--- a/mac/b12scan.mac
+++ b/mac/b12scan.mac
@@ -19,6 +19,8 @@
 
 # energy spectrum
 /VK/source/betaSpectrum true
+/VK/source/betaZ 6
+/VK/source/betaA 12
 /VK/source/betaEndPoint 13.37 MeV
 /VK/source/betaMinEnergy 0.00 MeV
 /VK/source/betaMaxEnergy 13.37 MeV
diff --git a/mac/f20scan.mac b/mac/f20scan.mac
index 7675d09..2333582 100644
--- a/mac/f20scan.mac
+++ b/mac/f20scan.mac
@@ -19,6 +19,8 @@
 
 # energy spectrum
 /VK/source/betaSpectrum true
+/VK/source/betaZ 10
+/VK/source/betaA 20
 /VK/source/betaEndPoint 5.391 MeV
 /VK/source/betaMinEnergy 0.00 MeV
 /VK/source/betaMaxEnergy 5.391 MeV
@@ -35,27 +37,7 @@
 
 # run simulations
 
-/VK/field/maxbfield 0.1504 tesla
-/VK/output/openFile output/f20scan/15mm_y5mm_v1/b1500.root                
-/run/beamOn 500000
-
 /VK/field/maxbfield 0.2000 tesla
-/VK/output/openFile output/f20scan/15mm_y5mm_v1/b2000.root                
+/VK/output/openFile output/f20scan/15mm_y5mm/b2000.root                
 /run/beamOn 500000
 
-/VK/field/maxbfield 0.2504 tesla
-/VK/output/openFile output/f20scan/15mm_y5mm_v1/b2500.root                
-/run/beamOn 500000
-
-/VK/field/maxbfield 0.3000 tesla
-/VK/output/openFile output/f20scan/15mm_y5mm_v1/b3000.root                
-/run/beamOn 500000
-
-/VK/field/maxbfield 0.3504 tesla
-/VK/output/openFile output/f20scan/15mm_y5mm_v1/b3500.root                
-/run/beamOn 500000
-
-/VK/source/betaMinEnergy 2.5 MeV
-/VK/field/maxbfield 0.4000 tesla
-/VK/output/openFile output/f20scan/15mm_y5mm_v1/b4000.root                
-/run/beamOn 500000
diff --git a/mac/forbidden.mac b/mac/forbidden.mac
index 12801ce..472ecd9 100644
--- a/mac/forbidden.mac
+++ b/mac/forbidden.mac
@@ -7,21 +7,31 @@
 # source dimensions
 /VK/source/diskShapedSource true
 /VK/source/innerDiameter 0 cm
-/VK/source/outerDiameter 9 mm
+/VK/source/outerDiameter 15 mm
+/VK/source/catcherFoilPosition 0 0 -299.5 mm
 
 # radiation
 /gps/particle e-
-/gps/position 0 0 -299.5 mm
+/gps/position 0 5 -299.5 mm
 /gps/ang/type iso
 /gps/ang/mintheta  90 deg
 /gps/ang/maxtheta 180 deg
 
 # energy spectrum
 /VK/source/betaSpectrum true
+/VK/source/betaZ 10
+/VK/source/betaA 20
 /VK/source/betaEndPoint 7.025 MeV
 /VK/source/betaMinEnergy 0.00 MeV
 /VK/source/betaMaxEnergy 7.025 MeV
 
+##/VK/source/betaShapeFuncA1  1.0
+##/VK/source/betaShapeFuncA2  80.
+##/VK/source/betaShapeFuncA3 -2.0
+##/VK/source/betaShapeFuncA4 -0.28
+##/VK/source/betaShapeFuncA5 -0.009
+##/VK/source/betaShapeFuncA6  0.0080
+
 # output 
 /VK/output/zeroSuppression true
 /VK/storeTrajectoryData false
@@ -29,14 +39,7 @@
 
 # run simulation
 
-#/VK/source/betaMinEnergy 5.00 MeV
-##/VK/field/maxbfield 0.48 tesla
-#/VK/field/maxbfield 0.543 tesla
-#/VK/output/openFile output/forbidden_cut5000_b543_k10.root 
-#/run/beamOn 1000000
-
 /VK/source/betaMinEnergy 5.00 MeV
-#/VK/field/maxbfield 0.50 tesla
-/VK/field/maxbfield 0.566 tesla
-/VK/output/openFile output/forbidden_cut5000_b566_k10.root
+/VK/field/maxbfield 0.50 tesla
+/VK/output/openFile output/forbidden_cut5000_b50_shape.root 
 /run/beamOn 1000000
diff --git a/mac/init.mac b/mac/init.mac
index d1c2391..91638e5 100644
--- a/mac/init.mac
+++ b/mac/init.mac
@@ -10,8 +10,8 @@
 /VK/physics/opticalPhysics false
 
 # setup
-/VK/detector/volumes/betaDetector true
-/VK/detector/volumes/betaDetectorV2 false
+/VK/detector/volumes/betaDetector false
+/VK/detector/volumes/betaDetectorV2 true
 /VK/detector/volumes/magnet true
 /VK/detector/volumes/catcherFoil true
 /VK/detector/volumes/shield true
diff --git a/mac/test.mac b/mac/test.mac
index ae817c1..704fc81 100644
--- a/mac/test.mac
+++ b/mac/test.mac
@@ -12,50 +12,33 @@
 
 # radiation
 /gps/particle e-
-/gps/position 0 7 -299.5 mm
+/gps/position 0 5 -299.5 mm
 /gps/ang/type iso
 /gps/ang/mintheta  90 deg
 /gps/ang/maxtheta 180 deg
 
 # energy spectrum
 /VK/source/betaSpectrum true
-/VK/source/betaEndPoint 5.391 MeV
+/VK/source/betaZ 10
+/VK/source/betaA 20
+/VK/source/betaEndPoint 7.025 MeV
 /VK/source/betaMinEnergy 0.00 MeV
-/VK/source/betaMaxEnergy 5.391 MeV
-/VK/source/betaShapeFuncA1 0.9674
-/VK/source/betaShapeFuncA2 -0.0021
-/VK/source/betaShapeFuncA3 0.01186
-/VK/source/betaShapeFuncA4 0.000333
-/VK/source/betaShapeFuncA5 0
-/VK/source/betaShapeFuncA6 0
+/VK/source/betaMaxEnergy 7.025 MeV
 
-# output 
-/VK/output/zeroSuppression true
-/VK/storeTrajectoryData false
-
-# run simulations
-
-/VK/field/maxbfield 0.1504 tesla
-/VK/output/openFile output/f20scan/15mm_y7mm_v1/b1500.root                
-/run/beamOn 500000
+#/VK/source/betaShapeFuncA1  1.0
+#/VK/source/betaShapeFuncA2  80.
+#/VK/source/betaShapeFuncA3 -2.0
+#/VK/source/betaShapeFuncA4 -0.28
+#/VK/source/betaShapeFuncA5 -0.009
+#/VK/source/betaShapeFuncA6  0.0080
 
-/VK/field/maxbfield 0.2000 tesla
-/VK/output/openFile output/f20scan/15mm_y7mm_v1/b2000.root                
-/run/beamOn 500000
-
-/VK/field/maxbfield 0.2504 tesla
-/VK/output/openFile output/f20scan/15mm_y7mm_v1/b2500.root                
-/run/beamOn 500000
+# output 
+/VK/output/zeroSuppression false
+/VK/storeTrajectoryData true
 
-/VK/field/maxbfield 0.3000 tesla
-/VK/output/openFile output/f20scan/15mm_y7mm_v1/b3000.root                
-/run/beamOn 500000
 
-/VK/field/maxbfield 0.3504 tesla
-/VK/output/openFile output/f20scan/15mm_y7mm_v1/b3500.root                
-/run/beamOn 500000
+# run simulation
 
-/VK/source/betaMinEnergy 2.5 MeV
-/VK/field/maxbfield 0.4000 tesla
-/VK/output/openFile output/f20scan/15mm_y7mm_v1/b4000.root                
-/run/beamOn 500000
+/VK/field/maxbfield 0.48 tesla
+/VK/output/openFile output/test1.root 
+/run/beamOn 10000
diff --git a/src/VKPrimaryGeneratorAction.cpp b/src/VKPrimaryGeneratorAction.cpp
index d4db072..5f7870e 100644
--- a/src/VKPrimaryGeneratorAction.cpp
+++ b/src/VKPrimaryGeneratorAction.cpp
@@ -56,6 +56,8 @@ VKPrimaryGeneratorAction::VKPrimaryGeneratorAction()
    fBetaEKinEndPoint(7.025*MeV),
    fBetaEKinMin(fBetaLowCut),
    fBetaEKinMax(fBetaEKinEndPoint),
+   fBetaZ(10),
+   fBetaA(20),
    fPosCentre(G4ThreeVector(0,0,0))   
 {
     // create messenger 
@@ -260,7 +262,7 @@ void VKPrimaryGeneratorAction::InitBetaSpectrum()
     shapeFunc = new TF1("shapeFunc", Form("%f + %f/x + %f*x + %f*pow(x,2) + %f*pow(x,3) + %f*pow(x,4)", shapeFuncCoeff[0], shapeFuncCoeff[1], shapeFuncCoeff[2], shapeFuncCoeff[3], shapeFuncCoeff[4], shapeFuncCoeff[5]), 0.5, 1E3);
 
     // allowed shape with correct Fermi function
-    fBetaSpectrum = new VKBetaSpectrum(6, 12, fBetaEKinEndPoint, shapeFunc);
+    fBetaSpectrum = new VKBetaSpectrum(fBetaZ, fBetaA, fBetaEKinEndPoint, shapeFunc);
     int npar = 0;
     fBetaSpectrumFunc = new TF1("betaSpectrum", fBetaSpectrum, &VKBetaSpectrum::Evaluate, fBetaLowCut/MeV, fBetaEKinEndPoint/MeV, npar);
 
@@ -327,6 +329,17 @@ void VKPrimaryGeneratorAction::SetBetaShapeFuncA6(double a6)
     shapeFuncCoeff[5] = a6;
 }
 
+void VKPrimaryGeneratorAction::SetBetaZ(int Z)      
+{ 
+    fBetaZ = Z;     
+}
+
+void VKPrimaryGeneratorAction::SetBetaA(int A)      
+{ 
+    fBetaA = A;     
+}
+
+
 double VKPrimaryGeneratorAction::GetBetaCoverage() const
 {
     double coverage = 1;
diff --git a/src/VKPrimaryGeneratorMessenger.cpp b/src/VKPrimaryGeneratorMessenger.cpp
index 70d26e9..aee4dfa 100644
--- a/src/VKPrimaryGeneratorMessenger.cpp
+++ b/src/VKPrimaryGeneratorMessenger.cpp
@@ -157,6 +157,18 @@ VKPrimaryGeneratorMessenger::VKPrimaryGeneratorMessenger(VKPrimaryGeneratorActio
     fShapeFuncA6Cmd -> SetParameterName("shapeFuncA0", false);
     fShapeFuncA6Cmd -> AvailableForStates(G4State_PreInit, G4State_Idle);
     fShapeFuncA6Cmd -> SetToBeBroadcasted(false);    
+
+    fBetaZCmd = new G4UIcmdWithAnInteger("/VK/source/betaZ", this);
+    fBetaZCmd -> SetGuidance("Set Z of beta-decay daughter nucleus");
+    fBetaZCmd -> SetParameterName("Z", false);
+    fBetaZCmd -> AvailableForStates(G4State_PreInit, G4State_Idle);
+    fBetaZCmd -> SetToBeBroadcasted(false);    
+
+    fBetaACmd = new G4UIcmdWithAnInteger("/VK/source/betaA", this);
+    fBetaACmd -> SetGuidance("Set A of beta-decay nucleus");
+    fBetaACmd -> SetParameterName("Z", false);
+    fBetaACmd -> AvailableForStates(G4State_PreInit, G4State_Idle);
+    fBetaACmd -> SetToBeBroadcasted(false);    
 }
 
 
@@ -180,6 +192,8 @@ VKPrimaryGeneratorMessenger::~VKPrimaryGeneratorMessenger()
     delete fShapeFuncA4Cmd;
     delete fShapeFuncA5Cmd;
     delete fShapeFuncA6Cmd;
+    delete fBetaZCmd;
+    delete fBetaACmd;
 }
 
 
@@ -257,4 +271,12 @@ void VKPrimaryGeneratorMessenger::SetNewValue(G4UIcommand* command, G4String new
     {
         fPrimGenAct -> SetBetaShapeFuncA6(fShapeFuncA6Cmd -> GetNewDoubleValue(newValue));
     }
+    else if (command == fBetaZCmd)
+    {
+        fPrimGenAct -> SetBetaZ(fBetaZCmd -> GetNewIntValue(newValue));
+    }
+    else if (command == fBetaACmd)
+    {
+        fPrimGenAct -> SetBetaA(fBetaACmd -> GetNewIntValue(newValue));
+    }
 }
-- 
GitLab