Merge pull request #44992 from lowette/sexaquark-simulation

Integration of sexaquark simulation code (for EXO-24-003)

SimG4Core/CustomPhysics/data/particles_sexaq_1p5_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     1.5                  # sexaq
+ -1020000020     1.5                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/data/particles_sexaq_1p7_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     1.7                  # sexaq
+ -1020000020     1.7                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/data/particles_sexaq_1p85_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     1.85                  # sexaq
+ -1020000020     1.85                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/data/particles_sexaq_1p8_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     1.8                  # sexaq
+ -1020000020     1.8                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/data/particles_sexaq_1p9_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     1.9                  # sexaq
+ -1020000020     1.9                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/data/particles_sexaq_2_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     2.0                  # sexaq
+ -1020000020     2.0                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/data/particles_sexaq_2p1_GeV.txt

@@ -0,0 +1,5 @@
+Block MASS   #
+#  PDG code    mass                 particle
+ 1020000020     2.1                  # sexaq
+ -1020000020     2.1                  # anti_sexaq
+Block

SimG4Core/CustomPhysics/interface/CMSAntiSQ.h

@@ -0,0 +1,23 @@
+#ifndef CMSAntiSQ_h
+#define CMSAntiSQ_h 1
+
+#include "globals.hh"
+#include "G4ios.hh"
+#include "G4ParticleDefinition.hh"
+
+// ######################################################################
+// ###                      ANTI-SEXAQUARK                            ###
+// ######################################################################
+
+class CMSAntiSQ : public G4ParticleDefinition {
+private:
+  static CMSAntiSQ* theInstance;
+  CMSAntiSQ() {}
+  ~CMSAntiSQ() {}
+
+public:
+  static CMSAntiSQ* Definition(double mass);
+  static CMSAntiSQ* AntiSQ(double mass);
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSQ.h

@@ -0,0 +1,24 @@
+
+#ifndef CMSSQ_h
+#define CMSSQ_h 1
+
+#include "globals.hh"
+#include "G4ios.hh"
+#include "G4ParticleDefinition.hh"
+
+// ######################################################################
+// ###                         SEXAQUARK                              ###
+// ######################################################################
+
+class CMSSQ : public G4ParticleDefinition {
+private:
+  static CMSSQ* theInstance;
+  CMSSQ() {}
+  ~CMSSQ() {}
+
+public:
+  static CMSSQ* Definition(double mass);
+  static CMSSQ* SQ(double mass);
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSQInelasticCrossSection.h

@@ -0,0 +1,30 @@
+
+#ifndef CMSSQInelasticCrossSection_h
+#define CMSSQInelasticCrossSection_h
+
+#include "globals.hh"
+#include "G4VCrossSectionDataSet.hh"
+
+class G4NistManager;
+class CMSSQ;
+class CMSAntiSQ;
+
+class CMSSQInelasticCrossSection : public G4VCrossSectionDataSet {
+public:
+  CMSSQInelasticCrossSection(double mass);
+
+  ~CMSSQInelasticCrossSection();
+
+  virtual G4bool IsElementApplicable(const G4DynamicParticle* aPart, G4int Z, const G4Material*);
+
+  virtual G4double GetElementCrossSection(const G4DynamicParticle*, G4int Z, const G4Material*);
+
+  G4double GetSQCrossSection(G4double kineticEnergy, G4int Z);
+
+private:
+  G4NistManager* nist;
+  CMSSQ* theSQ;
+  CMSAntiSQ* theAntiSQ;
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSQLoopProcess.h

@@ -0,0 +1,43 @@
+#ifndef CMSSQLoopProcess_h
+#define CMSSQLoopProcess_h 1
+
+#include "G4VContinuousProcess.hh"
+#include "globals.hh"
+#include "G4Track.hh"
+#include "G4ParticleChange.hh"
+
+class G4Step;
+class G4ParticleDefinition;
+
+class CMSSQLoopProcess : public G4VContinuousProcess {
+public:
+  CMSSQLoopProcess(const G4String& name = "SQLooper", G4ProcessType type = fUserDefined);
+  virtual ~CMSSQLoopProcess();
+
+public:
+  virtual G4VParticleChange* AlongStepDoIt(const G4Track&, const G4Step&);
+  virtual G4double AlongStepGetPhysicalInteractionLength(const G4Track& track,
+                                                         G4double previousStepSize,
+                                                         G4double currentMinimumStep,
+                                                         G4double& proposedSafety,
+                                                         G4GPILSelection* selection);
+  virtual void StartTracking(G4Track* aTrack);
+
+protected:
+  virtual G4double GetContinuousStepLimit(const G4Track& track,
+                                          G4double previousStepSize,
+                                          G4double currentMinimumStep,
+                                          G4double& currentSafety);
+
+private:
+  CMSSQLoopProcess(CMSSQLoopProcess&);
+  CMSSQLoopProcess& operator=(const CMSSQLoopProcess& right);
+
+protected:
+  G4ParticleChange* fParticleChange;
+
+private:
+  G4ThreeVector posini;
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSQLoopProcessDiscr.h

@@ -0,0 +1,42 @@
+#ifndef CMSSQLoopProcessDiscr_h
+#define CMSSQLoopProcessDiscr_h 1
+
+#include "G4VDiscreteProcess.hh"
+#include "globals.hh"
+#include "G4Track.hh"
+#include "G4ParticleChange.hh"
+#include "G4ParticleChangeForTransport.hh"
+#include "CMSSQ.h"
+#include "CMSAntiSQ.h"
+
+class G4Step;
+class G4ParticleDefinition;
+
+class CMSSQLoopProcessDiscr : public G4VDiscreteProcess {
+public:
+  CMSSQLoopProcessDiscr(double mass, const G4String& name = "SQLooper", G4ProcessType type = fUserDefined);
+  virtual ~CMSSQLoopProcessDiscr();
+
+public:
+  virtual G4VParticleChange* PostStepDoIt(const G4Track&, const G4Step&);
+  virtual G4double PostStepGetPhysicalInteractionLength(const G4Track& track,
+                                                        G4double previousStepSize,
+                                                        G4ForceCondition* condition);
+  virtual G4double GetMeanFreePath(const G4Track&, G4double, G4ForceCondition*);
+  void SetTimeLimit(G4double);
+  virtual void StartTracking(G4Track* aTrack);
+
+private:
+  CMSSQLoopProcessDiscr(CMSSQLoopProcessDiscr&);
+  CMSSQLoopProcessDiscr& operator=(const CMSSQLoopProcessDiscr& right);
+
+protected:
+  G4ParticleChange* fParticleChange;
+  double GenMass;
+
+private:
+  G4ThreeVector posini;
+  G4double globaltimeini;
+};
+
+#endif

SimG4Core/CustomPhysics/interface/CMSSQNeutronAnnih.h

@@ -0,0 +1,30 @@
+
+#ifndef CMSSQNeutronAnnih_h
+#define CMSSQNeutronAnnih_h 1
+
+#include "globals.hh"
+#include "G4HadronicInteraction.hh"
+#include "G4HadProjectile.hh"
+#include "G4Nucleus.hh"
+#include "G4IonTable.hh"
+
+class G4ParticleDefinition;
+
+class CMSSQNeutronAnnih : public G4HadronicInteraction {
+public:
+  CMSSQNeutronAnnih(double mass);
+
+  ~CMSSQNeutronAnnih() override;
+
+  G4double momDistr(G4double x_in);
+
+  G4HadFinalState* ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus) override;
+
+private:
+  G4ParticleDefinition* theSQ;
+  G4ParticleDefinition* theK0S;
+  G4ParticleDefinition* theAntiL;
+  G4ParticleDefinition* theProton;
+};
+
+#endif

SimG4Core/CustomPhysics/src/CMSAntiSQ.cc

@@ -0,0 +1,43 @@
+
+#include "SimG4Core/CustomPhysics/interface/CMSAntiSQ.h"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleTable.hh"
+
+#include "G4PhaseSpaceDecayChannel.hh"
+#include "G4DecayTable.hh"
+
+// ######################################################################
+// ###                       ANTI-SEXAQUARK                           ###
+// ######################################################################
+
+CMSAntiSQ* CMSAntiSQ::theInstance = 0;
+
+CMSAntiSQ* CMSAntiSQ::Definition(double mass) {
+  if (theInstance != 0)
+    return theInstance;
+  const G4String name = "anti_sexaq";
+  // search in particle table]
+  G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
+  G4ParticleDefinition* anInstance = pTable->FindParticle(name);
+  if (anInstance == 0) {
+    // create particle
+    //
+    //    Arguments for constructor are as follows
+    //               name             mass          width         charge
+    //             2*spin           parity  C-conjugation
+    //          2*Isospin       2*Isospin3       G-parity
+    //               type    lepton number  baryon number   PDG encoding
+    //             stable         lifetime    decay table
+    //             shortlived      subType    anti_encoding
+
+    anInstance = new G4ParticleDefinition(
+        name, mass, 0, 0.0, 0, +1, 0, 0, 0, 0, "baryon", 0, -2, -1020000020, true, -1.0, nullptr, false, "sexaq");
+  }
+  theInstance = reinterpret_cast<CMSAntiSQ*>(anInstance);
+  return theInstance;
+}
+
+CMSAntiSQ* CMSAntiSQ::AntiSQ(double mass) {
+  return Definition(mass * GeV);  // will use correct mass if instance exists
+}

SimG4Core/CustomPhysics/src/CMSSQ.cc

@@ -0,0 +1,43 @@
+
+#include "SimG4Core/CustomPhysics/interface/CMSSQ.h"
+#include "G4PhysicalConstants.hh"
+#include "G4SystemOfUnits.hh"
+#include "G4ParticleTable.hh"
+
+#include "G4PhaseSpaceDecayChannel.hh"
+#include "G4DecayTable.hh"
+
+// ######################################################################
+// ###                          SEXAQUARK                             ###
+// ######################################################################
+
+CMSSQ* CMSSQ::theInstance = 0;
+
+CMSSQ* CMSSQ::Definition(double mass) {
+  if (theInstance != 0)
+    return theInstance;
+  const G4String name = "sexaq";
+  // search in particle table]
+  G4ParticleTable* pTable = G4ParticleTable::GetParticleTable();
+  G4ParticleDefinition* anInstance = pTable->FindParticle(name);
+  if (anInstance == 0) {
+    // create particle
+    //
+    //    Arguments for constructor are as follows
+    //               name             mass          width         charge
+    //             2*spin           parity  C-conjugation
+    //          2*Isospin       2*Isospin3       G-parity
+    //               type    lepton number  baryon number   PDG encoding
+    //             stable         lifetime    decay table
+    //             shortlived      subType    anti_encoding
+
+    anInstance = new G4ParticleDefinition(
+        name, mass, 0, 0.0, 0, +1, 0, 0, 0, 0, "baryon", 0, +2, 1020000020, true, -1.0, nullptr, false, "sexaq");
+  }
+  theInstance = reinterpret_cast<CMSSQ*>(anInstance);
+  return theInstance;
+}
+
+CMSSQ* CMSSQ::SQ(double mass) {
+  return Definition(mass * GeV);  // will use correct mass if instance exists
+}

SimG4Core/CustomPhysics/src/CMSSQInelasticCrossSection.cc

@@ -0,0 +1,42 @@
+
+#include "G4SystemOfUnits.hh"
+#include "G4DynamicParticle.hh"
+#include "G4NistManager.hh"
+
+#include "SimG4Core/CustomPhysics/interface/CMSSQ.h"
+#include "SimG4Core/CustomPhysics/interface/CMSAntiSQ.h"
+#include "SimG4Core/CustomPhysics/interface/CMSSQInelasticCrossSection.h"
+
+CMSSQInelasticCrossSection::CMSSQInelasticCrossSection(double mass) : G4VCrossSectionDataSet("SQ-neutron") {
+  nist = G4NistManager::Instance();
+  theSQ = CMSSQ::SQ(mass);
+  theAntiSQ = CMSAntiSQ::AntiSQ(mass);
+}
+
+CMSSQInelasticCrossSection::~CMSSQInelasticCrossSection() {}
+
+G4bool CMSSQInelasticCrossSection::IsElementApplicable(const G4DynamicParticle* aPart, G4int Z, const G4Material*) {
+  return true;
+}
+
+G4double CMSSQInelasticCrossSection::GetElementCrossSection(const G4DynamicParticle* aPart,
+                                                            G4int Z,
+                                                            const G4Material*) {
+  // return zero for particle instead of antiparticle
+  // sexaquark interaction with matter expected really tiny
+  if (aPart->GetDefinition() != theAntiSQ)
+    return 0;
+
+  //I don't want to interact on hydrogen
+  if (Z <= 1) {
+    return 0.0;
+  }
+
+  // get the atomic weight (to estimate nr neutrons)
+  G4double A = nist->GetAtomicMassAmu(Z);
+
+  // put the X section low for the antiS to get a flat interaction rate,
+  // but also make it scale with the number of neutrons in the material
+  // because we are going to interact on neutrons, not on protons
+  return (100. * (A - (G4double)Z) / (G4double)Z) * millibarn;
+}