* add ability to set the target polarization in the simulation;

  separate settings are provided for the electron polarization
  and the nuclear polarization via new interactive commands
  /hdgeant4/targetNuclearPolarization <polx> <poly> <polz> and
  /hdgeant4/targetElectronPolarization <polx> <poly> <polz>. [rtj]
* add ability to change the diamond crystal orientation during
  the simulation via new interactive G4 interface command
  /hdgeant4/diamondAngles <thetax> <thetay> <thetaz>. [rtj]
* these new commands required that hooks be added in various
  places within the existing code for the getters and setters. [rtj]

Author: rjones30

GlueXBeamConversionProcess.cc

@@ -62,6 +62,8 @@ const double CobremsGeneration::me = 0.510998910e-3;
 const double CobremsGeneration::alpha = 7.2973525698e-3;
 const double CobremsGeneration::hbarc = 0.1973269718e-15;
 
+std::map<CobremsGeneration*, int> CobremsGeneration::CobremsGenerators;
+
 CobremsGeneration::CobremsGeneration(double Emax_GeV, double Epeak_GeV)
 {
    // Unique constructor for this class, initialize for the given
@@ -91,6 +93,8 @@ CobremsGeneration::CobremsGeneration(double Emax_GeV, double Epeak_GeV)
              << " primary coherent edge: " << Epeak_GeV << " GeV"
              << std::endl;
 #endif
+
+   CobremsGenerators[this] = 1;
 }
 
 void CobremsGeneration::updateTargetOrientation()
@@ -468,7 +472,10 @@ CobremsGeneration &CobremsGeneration::operator=(const CobremsGeneration &src)
    return *this;
 }
 
-CobremsGeneration::~CobremsGeneration() { }
+CobremsGeneration::~CobremsGeneration()
+{
+   CobremsGenerators.erase(this);
+}
 
 double CobremsGeneration::CoherentEnhancement(double x)
 {

src/CobremsGeneration.cc

@@ -24,6 +24,7 @@
 
 #include <string>
 #include <vector>
+#include <map>
 
 #if BOOST_PYTHON_WRAPPING
 #include <boost/python.hpp>
@@ -36,6 +37,8 @@ class CobremsGeneration {
    CobremsGeneration &operator=(const CobremsGeneration &src);
    ~CobremsGeneration();
 
+   static std::map<CobremsGeneration*, int> CobremsGenerators;
+
    void setBeamEnergy(double Ebeam_GeV);
    void setBeamErms(double Erms_GeV);
    void setBeamEmittance(double emit_m_r);
@@ -274,8 +277,8 @@ inline void CobremsGeneration::setTargetThetaz(double thetaz) {
 }
 
 inline void CobremsGeneration::setTargetOrientation(double thetax,
-                                                   double thetay,
-                                                   double thetaz) {
+                                                    double thetay,
+                                                    double thetaz) {
    fTargetThetax = thetax;
    fTargetThetay = thetay;
    fTargetThetaz = thetaz;

src/CobremsGeneration.hh

@@ -134,7 +134,8 @@ GlueXBeamConversionProcess::GlueXBeamConversionProcess(const G4String &name,
    fAdaptiveSampler(0),
    isInitialised(false),
    fTargetZ(0),
-   fTargetA(0)
+   fTargetA(0),
+   fTargetPol{{0,0,0},{0,0,0}}
 {
    verboseLevel = 0;
 
@@ -1214,7 +1215,9 @@ void GlueXBeamConversionProcess::GenerateBetheHeitlerProcess(const G4Step &step)
    pOut.AllPol();
    nOut.AllPol();
    const TThreeVectorReal zero3vector(0,0,0);
-   nIn.SetPol(TThreeVectorReal(zero3vector));
+   nIn.SetPol(TThreeVectorReal(fTargetPol.nucl[0],
+                               fTargetPol.nucl[1],
+                               fTargetPol.nucl[2]));
    nIn.SetMom(TThreeVectorReal(zero3vector));
    const G4Track *track = step.GetTrack();
    LDouble_t kin = track->GetKineticEnergy()/GeV;
@@ -1706,7 +1709,9 @@ void GlueXBeamConversionProcess::GenerateTripletProcess(const G4Step &step)
    pOut.AllPol();
    eOut3.AllPol();
    const TThreeVectorReal zero3vector(0,0,0);
-   eIn.SetPol(TThreeVectorReal(zero3vector));
+   eIn.SetPol(TThreeVectorReal(fTargetPol.elec[0],
+                               fTargetPol.elec[1],
+                               fTargetPol.elec[2]));
    eIn.SetMom(TThreeVectorReal(zero3vector));
    const G4Track *track = step.GetTrack();
    LDouble_t kin = track->GetKineticEnergy()/GeV;

src/GlueXBeamConversionProcess.cc

@@ -85,6 +85,11 @@ class GlueXBeamConversionProcess: public G4VEmProcess
    void prepareAdaptiveSampler();
 
    void setConverterMaterial(double Z, double A);
+   void getConverterMaterial(double &Z, double &A) const;
+
+ public:
+   void setTargetPolarization(double nucl[3], double elec[3]);
+   void getTargetPolarization(double nucl[3], double elec[3]) const;
 
  private:
    GlueXBeamConversionProcess() = delete;
@@ -101,12 +106,44 @@ class GlueXBeamConversionProcess: public G4VEmProcess
 
    double fTargetZ;
    double fTargetA;
+   struct target_polarization_t {
+      double nucl[3];
+      double elec[3];
+   } fTargetPol;
 };
 
+inline void GlueXBeamConversionProcess::getConverterMaterial(double &Z, double &A) const
+{
+   Z = fTargetZ;
+   A = fTargetA;
+}
+
 inline void GlueXBeamConversionProcess::setConverterMaterial(double Z, double A)
 {
    fTargetZ = Z;
    fTargetA = A;
 }
 
+inline void GlueXBeamConversionProcess::getTargetPolarization(double nucl[3],
+                                                              double elec[3]) const
+{
+   nucl[0] = fTargetPol.nucl[0];
+   nucl[1] = fTargetPol.nucl[1];
+   nucl[2] = fTargetPol.nucl[2];
+   elec[0] = fTargetPol.elec[0];
+   elec[1] = fTargetPol.elec[1];
+   elec[2] = fTargetPol.elec[2];
+}
+
+inline void GlueXBeamConversionProcess::setTargetPolarization(double nucl[3],
+                                                              double elec[3])
+{
+   fTargetPol.nucl[0] = nucl[0];
+   fTargetPol.nucl[1] = nucl[1];
+   fTargetPol.nucl[2] = nucl[2];
+   fTargetPol.elec[0] = elec[0];
+   fTargetPol.elec[1] = elec[1];
+   fTargetPol.elec[2] = elec[2];
+}
+
 #endif

src/GlueXBeamConversionProcess.hh

@@ -12,7 +12,12 @@
 #include "G4UIcmdWithAString.hh"
 #include "G4UIcmdWithADoubleAndUnit.hh"
 #include "G4UIcmdWithoutParameter.hh"
+#include "G4UIcmdWith3Vector.hh"
 #include "G4GeometryManager.hh"
+#include "G4RunManager.hh"
+#include "GlueXPhysicsList.hh"
+#include "GlueXPrimaryGeneratorAction.hh"
+#include "GlueXBeamConversionProcess.hh"
 
 GlueXDetectorMessenger::GlueXDetectorMessenger(GlueXDetectorConstruction* myDet)
 :myDetector(myDet)
@@ -32,7 +37,7 @@ GlueXDetectorMessenger::GlueXDetectorMessenger(GlueXDetectorConstruction* myDet)
 
   StepMaxCmd = new G4UIcmdWithADoubleAndUnit("/hdgeant4/det/stepMax",this);
   StepMaxCmd->SetGuidance("Define a step max");
-  StepMaxCmd->SetParameterName("stepMax",false);
+  StepMaxCmd->SetParameterName("stepMax",true,true);
   StepMaxCmd->SetUnitCategory("Length");
   StepMaxCmd->AvailableForStates(G4State_Idle);
 
@@ -45,6 +50,30 @@ GlueXDetectorMessenger::GlueXDetectorMessenger(GlueXDetectorConstruction* myDet)
   CloseGeomCmd = new G4UIcmdWithoutParameter("/hdgeant4/closeGeometry",this);
   CloseGeomCmd->SetGuidance("Explicitly close the geometry.");
   CloseGeomCmd->AvailableForStates(G4State_Idle);
+
+  RadiatorAnglesCmd = new G4UIcmdWith3Vector("/hdgeant4/diamondAngles",this);
+  RadiatorAnglesCmd->SetGuidance("Control the diamond radiator orientation.");
+  RadiatorAnglesCmd->SetGuidance("Arguments are theta_x, theta_y, theta_z");
+  RadiatorAnglesCmd->SetGuidance(" in radians, applied in the order theta_x,");
+  RadiatorAnglesCmd->SetGuidance(" followed by theta_y, followed by theta_z.");
+  RadiatorAnglesCmd->SetParameterName("theta_x","theta_y","theta_z",true,true);
+  RadiatorAnglesCmd->AvailableForStates(G4State_Idle);
+
+  TargetNuclPolCmd = new G4UIcmdWith3Vector("/hdgeant4/targetNuclearPolarization",this);
+  TargetNuclPolCmd->SetGuidance("Control the primary target nuclear polarization.");
+  TargetNuclPolCmd->SetGuidance("Arguments are pol_x, pol_y, pol_z");
+  TargetNuclPolCmd->SetGuidance(" in fractions, normalized to a unit vector");
+  TargetNuclPolCmd->SetGuidance(" for 100% polarization.");
+  TargetNuclPolCmd->SetParameterName("pol_x","pol_y","pol_z",true,true);
+  TargetNuclPolCmd->AvailableForStates(G4State_Idle);
+
+  TargetElecPolCmd = new G4UIcmdWith3Vector("/hdgeant4/targetElectronPolarization",this);
+  TargetElecPolCmd->SetGuidance("Control the primary target electron polarization.");
+  TargetElecPolCmd->SetGuidance("Arguments are pol_x, pol_y, pol_z");
+  TargetElecPolCmd->SetGuidance(" in fractions, normalized to a unit vector");
+  TargetElecPolCmd->SetGuidance(" for 100% polarization.");
+  TargetElecPolCmd->SetParameterName("pol_x","pol_y","pol_z",true,true);
+  TargetElecPolCmd->AvailableForStates(G4State_Idle);
 }
 
 GlueXDetectorMessenger::~GlueXDetectorMessenger()
@@ -53,6 +82,10 @@ GlueXDetectorMessenger::~GlueXDetectorMessenger()
   delete StepMaxCmd;
   delete detDir;
   delete hdgeant4Dir;
+  delete OpenGeomCmd;
+  delete RadiatorAnglesCmd;
+  delete TargetNuclPolCmd;
+  delete TargetElecPolCmd;
 }
 
 void GlueXDetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
@@ -62,11 +95,121 @@ void GlueXDetectorMessenger::SetNewValue(G4UIcommand* command,G4String newValue)
   }
   else if (command == StepMaxCmd) {
     myDetector->SetMaxStep(StepMaxCmd->GetNewDoubleValue(newValue));
+    std::cout << "MaxStep is now " << StepMaxCmd->GetCurrentValue() << std::endl;
   }
   else if (command == OpenGeomCmd) {
     G4GeometryManager::GetInstance()->OpenGeometry();
   }
   else if (command == CloseGeomCmd) {
     G4GeometryManager::GetInstance()->CloseGeometry();
   }
+  else if (command == RadiatorAnglesCmd) {
+    G4ThreeVector angles(RadiatorAnglesCmd->GetNew3VectorValue(newValue));
+    std::map<CobremsGeneration*, int>::iterator igen;
+    int assigned(0);
+    for (igen = CobremsGeneration::CobremsGenerators.begin();
+         igen != CobremsGeneration::CobremsGenerators.end();
+         igen++)
+    {
+       igen->first->setTargetOrientation(angles[0], angles[1], angles[2]);
+       std::cout << "radiator angles are now " 
+                 << RadiatorAnglesCmd->GetCurrentValue() << std::endl;
+       assigned += igen->second;
+    }
+    if (!assigned) {
+       std::cerr << "CobremsGeneration has not been started yet, you must start"
+                 << " up the photon beam generator before attempting to set"
+                 << " the radiator orientation." << std::endl;
+    }
+  }
+  else if (command == TargetNuclPolCmd) {
+    G4ThreeVector polar(RadiatorAnglesCmd->GetNew3VectorValue(newValue));
+    const GlueXPhysicsList *phy = dynamic_cast<const GlueXPhysicsList*>
+       (G4RunManager::GetRunManager()->GetUserPhysicsList());
+    GlueXBeamConversionProcess *con = (GlueXBeamConversionProcess*)
+                                      phy->getBeamConversionProcess();
+    if (con != 0 && polar.mag() <= 1) {
+       double nucl[3];
+       double elec[3];
+       con->getTargetPolarization(nucl, elec);
+       nucl[0] = polar[0];
+       nucl[1] = polar[1];
+       nucl[2] = polar[2];
+       con->setTargetPolarization(nucl, elec);
+    }
+    else {
+       std::cerr << "polarization vector normalization >1, assignment failed." 
+                 << std::endl;
+    }
+    std::cout << "target nuclear polarization vector is now " 
+              << TargetNuclPolCmd->GetCurrentValue() << std::endl;
+  }
+  else if (command == TargetElecPolCmd) {
+    G4ThreeVector polar(RadiatorAnglesCmd->GetNew3VectorValue(newValue));
+    const GlueXPhysicsList *phy = dynamic_cast<const GlueXPhysicsList*>
+       (G4RunManager::GetRunManager()->GetUserPhysicsList());
+    GlueXBeamConversionProcess *con = (GlueXBeamConversionProcess*)
+                                      phy->getBeamConversionProcess();
+    if (con != 0 && polar.mag() <= 1) {
+       double nucl[3];
+       double elec[3];
+       con->getTargetPolarization(nucl, elec);
+       elec[0] = polar[0];
+       elec[1] = polar[1];
+       elec[2] = polar[2];
+       con->setTargetPolarization(nucl, elec);
+    }
+    else {
+       std::cerr << "polarization vector normalization >1, assignment failed."
+                 << std::endl;
+    }
+    std::cout << "target electron polarization vector is now " 
+              << TargetElecPolCmd->GetCurrentValue() << std::endl;
+  }
+}
+
+G4String GlueXDetectorMessenger::GetCurrentValue(G4UIcommand* command)
+{
+  G4String cv;
+  if (command == StepMaxCmd) {
+    cv = StepMaxCmd->ConvertToString(myDetector->GetMaxStep(cm));
+  }
+  else if (command == RadiatorAnglesCmd) {
+    G4ThreeVector angles;
+    std::map<CobremsGeneration*, int>::iterator igen;
+    for (igen = CobremsGeneration::CobremsGenerators.begin();
+         igen != CobremsGeneration::CobremsGenerators.end();
+         igen++)
+    {
+       angles[0] = igen->first->getTargetThetax();
+       angles[1] = igen->first->getTargetThetay();
+       angles[2] = igen->first->getTargetThetaz();
+       cv = RadiatorAnglesCmd->ConvertToString(angles);
+    }
+  }
+  else if (command == TargetNuclPolCmd) {
+    const GlueXPhysicsList *phy = dynamic_cast<const GlueXPhysicsList*>
+                   (G4RunManager::GetRunManager()->GetUserPhysicsList());
+    if (phy != 0) {
+      GlueXBeamConversionProcess *con = (GlueXBeamConversionProcess*)
+                                        phy->getBeamConversionProcess();
+      double nucl[3];
+      double elec[3];
+      con->getTargetPolarization(nucl, elec);
+      cv = TargetNuclPolCmd->ConvertToString(G4ThreeVector(nucl[0], nucl[1], nucl[2]));
+    }
+  }
+  else if (command == TargetElecPolCmd) {
+    const GlueXPhysicsList *phy = dynamic_cast<const GlueXPhysicsList*>
+                   (G4RunManager::GetRunManager()->GetUserPhysicsList());
+    if (phy != 0) {
+      GlueXBeamConversionProcess *con = (GlueXBeamConversionProcess*)
+                                        phy->getBeamConversionProcess();
+      double nucl[3];
+      double elec[3];
+      con->getTargetPolarization(nucl, elec);
+      cv = TargetElecPolCmd->ConvertToString(G4ThreeVector(elec[0], elec[1], elec[2]));
+    }
+  }
+  return cv;
 }

src/GlueXDetectorMessenger.cc

@@ -17,6 +17,7 @@ class G4UIdirectory;
 class G4UIcmdWithAString;
 class G4UIcmdWithADoubleAndUnit;
 class G4UIcmdWithoutParameter;
+class G4UIcmdWith3Vector;
 
 class GlueXDetectorMessenger: public G4UImessenger
 {
@@ -27,6 +28,7 @@ class GlueXDetectorMessenger: public G4UImessenger
    ~GlueXDetectorMessenger();
 
     void SetNewValue(G4UIcommand*, G4String);
+    G4String GetCurrentValue(G4UIcommand*);
 
   private:
     GlueXDetectorConstruction* myDetector;
@@ -37,6 +39,9 @@ class GlueXDetectorMessenger: public G4UImessenger
     G4UIcmdWithADoubleAndUnit* StepMaxCmd;    
     G4UIcmdWithoutParameter*   OpenGeomCmd;
     G4UIcmdWithoutParameter*   CloseGeomCmd;
+    G4UIcmdWith3Vector*        RadiatorAnglesCmd;
+    G4UIcmdWith3Vector*        TargetNuclPolCmd;
+    G4UIcmdWith3Vector*        TargetElecPolCmd;
 };
 
 #endif

src/GlueXDetectorMessenger.hh

@@ -53,6 +53,10 @@ class GlueXPhysicsList: public G4VModularPhysicsList
    virtual void DoProcessReordering();
    virtual void CheckProcessOrdering();
 
+   G4VEmProcess *getBeamConversionProcess() const;
+   G4VEmProcess *getBernardConversionProcess() const;
+   G4OpticalPhysics *getOpticalPhysicsProcess() const;
+
  protected:
    GlueXUserOptions *fOptions;
 
@@ -72,4 +76,19 @@ class GlueXPhysicsList: public G4VModularPhysicsList
 #endif
 };
 
+inline G4VEmProcess *GlueXPhysicsList::getBeamConversionProcess() const
+{
+   return fBeamConversion;
+}
+
+inline G4VEmProcess *GlueXPhysicsList::getBernardConversionProcess() const
+{
+   return fBernardConversion;
+}
+
+inline G4OpticalPhysics *GlueXPhysicsList::getOpticalPhysicsProcess() const
+{
+   return fOpticalPhysics;
+}
+
 #endif

src/GlueXPhysicsList.hh

@@ -569,7 +569,7 @@ GlueXPrimaryGeneratorAction::~GlueXPrimaryGeneratorAction()
    }
 }
 
-const CobremsGeneration *GlueXPrimaryGeneratorAction::GetCobremsGeneration() const
+const CobremsGeneration *GlueXPrimaryGeneratorAction::GetCobremsGeneration()
 {
    return fCobremsGeneration;
 }