diff --git a/Makefile b/Makefile index 74aa717b6..0477f544e 100644 --- a/Makefile +++ b/Makefile @@ -97,8 +97,8 @@ physics-neutrino-scattering-modes: FORCE cd ${GENIE}/src/Physics/DeepInelastic/EventGen && $(MAKE) && \ cd ${GENIE}/src/Physics/Diffractive/XSection && $(MAKE) && \ cd ${GENIE}/src/Physics/Diffractive/EventGen && $(MAKE) && \ - cd ${GENIE}/src/Physics/GlashowResonance/XSection && $(MAKE) && \ - cd ${GENIE}/src/Physics/GlashowResonance/EventGen && $(MAKE) && \ + cd ${GENIE}/src/Physics/HELepton/XSection && $(MAKE) && \ + cd ${GENIE}/src/Physics/HELepton/EventGen && $(MAKE) && \ cd ${GENIE}/src/Physics/InverseBetaDecay/XSection && $(MAKE) && \ cd ${GENIE}/src/Physics/InverseBetaDecay/EventGen && $(MAKE) && \ cd ${GENIE}/src/Physics/Multinucleon/XSection && $(MAKE) && \ @@ -393,9 +393,9 @@ endif mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/Diffractive mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/Diffractive/XSection mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/Diffractive/EventGen - mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/GlashowResonance - mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/GlashowResonance/XSection - mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/GlashowResonance/EventGen + mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/HELepton + mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/HELepton/XSection + mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/HELepton/EventGen mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/Hadronization mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/HadronTensors mkdir ${GENIE_INC_INSTALLATION_PATH}/Physics/HadronTransport @@ -465,8 +465,8 @@ copy-install-files: FORCE cd ${GENIE}/src/Physics/DeepInelastic/EventGen && $(MAKE) install && \ cd ${GENIE}/src/Physics/Diffractive/XSection && $(MAKE) install && \ cd ${GENIE}/src/Physics/Diffractive/EventGen && $(MAKE) install && \ - cd ${GENIE}/src/Physics/GlashowResonance/XSection && $(MAKE) install && \ - cd ${GENIE}/src/Physics/GlashowResonance/EventGen && $(MAKE) install && \ + cd ${GENIE}/src/Physics/HELepton/XSection && $(MAKE) install && \ + cd ${GENIE}/src/Physics/HELepton/EventGen && $(MAKE) install && \ cd ${GENIE}/src/Physics/Hadronization && $(MAKE) install && \ cd ${GENIE}/src/Physics/HadronTransport && $(MAKE) install && \ cd ${GENIE}/src/Physics/HadronTensors && $(MAKE) install && \ @@ -527,8 +527,8 @@ purge: FORCE cd ${GENIE}/src/Physics/DeepInelastic/EventGen && $(MAKE) purge && \ cd ${GENIE}/src/Physics/Diffractive/XSection && $(MAKE) purge && \ cd ${GENIE}/src/Physics/Diffractive/EventGen && $(MAKE) purge && \ - cd ${GENIE}/src/Physics/GlashowResonance/XSection && $(MAKE) purge && \ - cd ${GENIE}/src/Physics/GlashowResonance/EventGen && $(MAKE) purge && \ + cd ${GENIE}/src/Physics/HELepton/XSection && $(MAKE) purge && \ + cd ${GENIE}/src/Physics/HELepton/EventGen && $(MAKE) purge && \ cd ${GENIE}/src/Physics/Hadronization && $(MAKE) purge && \ cd ${GENIE}/src/Physics/HadronTensors && $(MAKE) purge && \ cd ${GENIE}/src/Physics/HadronTransport && $(MAKE) purge && \ @@ -590,8 +590,8 @@ clean-files: FORCE cd ${GENIE}/src/Physics/DeepInelastic/EventGen && $(MAKE) clean && \ cd ${GENIE}/src/Physics/Diffractive/XSection && $(MAKE) clean && \ cd ${GENIE}/src/Physics/Diffractive/EventGen && $(MAKE) clean && \ - cd ${GENIE}/src/Physics/GlashowResonance/XSection && $(MAKE) clean && \ - cd ${GENIE}/src/Physics/GlashowResonance/EventGen && $(MAKE) clean && \ + cd ${GENIE}/src/Physics/HELepton/XSection && $(MAKE) clean && \ + cd ${GENIE}/src/Physics/HELepton/EventGen && $(MAKE) clean && \ cd ${GENIE}/src/Physics/Hadronization && $(MAKE) clean && \ cd ${GENIE}/src/Physics/HadronTensors && $(MAKE) clean && \ cd ${GENIE}/src/Physics/HadronTransport && $(MAKE) clean && \ @@ -678,8 +678,8 @@ endif cd ${GENIE}/src/Physics/DeepInelastic/EventGen && $(MAKE) distclean && \ cd ${GENIE}/src/Physics/Diffractive/XSection && $(MAKE) distclean && \ cd ${GENIE}/src/Physics/Diffractive/EventGen && $(MAKE) distclean && \ - cd ${GENIE}/src/Physics/GlashowResonance/XSection && $(MAKE) distclean && \ - cd ${GENIE}/src/Physics/GlashowResonance/EventGen && $(MAKE) distclean && \ + cd ${GENIE}/src/Physics/HELepton/XSection && $(MAKE) distclean && \ + cd ${GENIE}/src/Physics/HELepton/EventGen && $(MAKE) distclean && \ cd ${GENIE}/src/Physics/Hadronization && $(MAKE) distclean && \ cd ${GENIE}/src/Physics/HadronTransport && $(MAKE) distclean && \ cd ${GENIE}/src/Physics/HadronTensors && $(MAKE) distclean && \ diff --git a/config/EventGenerator.xml b/config/EventGenerator.xml index 8f99389c2..5fe64420a 100644 --- a/config/EventGenerator.xml +++ b/config/EventGenerator.xml @@ -687,7 +687,7 @@ XSecModel alg Yes Cross section model used at the thread 5 genie::InitialStateAppender/Default genie::VertexGenerator/Default - genie::HEDISKinematicsGenerator/Default + genie::HEDISKinematicsGenerator/Default genie::HEDISGenerator/Default genie::UnstableParticleDecayer/BeforeHadronTransport genie::HEDISInteractionListGenerator/CC-Default @@ -698,7 +698,7 @@ XSecModel alg Yes Cross section model used at the thread 5 genie::InitialStateAppender/Default genie::VertexGenerator/Default - genie::HEDISKinematicsGenerator/Default + genie::HEDISKinematicsGenerator/Default genie::HEDISGenerator/Default genie::UnstableParticleDecayer/BeforeHadronTransport genie::HEDISInteractionListGenerator/NC-Default @@ -706,47 +706,125 @@ XSecModel alg Yes Cross section model used at the thread - 5 - genie::InitialStateAppender/Default - genie::VertexGenerator/Default - genie::GLRESKinematicsGenerator/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/BeforeHadronTransport - genie::GLRESInteractionListGenerator/Mu-Default + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::GLRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/GLRES-Mu-Default - 5 - genie::InitialStateAppender/Default - genie::VertexGenerator/Default - genie::GLRESKinematicsGenerator/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/BeforeHadronTransport - genie::GLRESInteractionListGenerator/Tau-Default + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::GLRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/GLRES-Tau-Default - 5 - genie::InitialStateAppender/Default - genie::VertexGenerator/Default - genie::GLRESKinematicsGenerator/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/BeforeHadronTransport - genie::GLRESInteractionListGenerator/Ele-Default + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::GLRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/GLRES-Ele-Default - 5 - genie::InitialStateAppender/Default - genie::VertexGenerator/Default - genie::GLRESKinematicsGenerator/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/BeforeHadronTransport - genie::GLRESInteractionListGenerator/Had-Default + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::GLRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/GLRES-Had-Default + + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::HENuElGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/HENuEl-CC-Default + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::HENuElGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/HENuEl-NC-Default + + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::PhotonRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/PhotonRES-Mu-Default + + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::PhotonRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/PhotonRES-Ele-Default + + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::PhotonRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/PhotonRES-Tau-Default + + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::PhotonRESGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/PhotonRES-Had-Default + + + + + 5 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::HELeptonKinematicsGenerator/Default + genie::PhotonCOHGenerator/Default + genie::UnstableParticleDecayer/BeforeHadronTransport + genie::HELeptonInteractionListGenerator/PhotonCOH-Default + + diff --git a/config/EventGeneratorListAssembler.xml b/config/EventGeneratorListAssembler.xml index 25e39e350..b4f101090 100644 --- a/config/EventGeneratorListAssembler.xml +++ b/config/EventGeneratorListAssembler.xml @@ -57,45 +57,6 @@ Generator-%d alg No genie::EventGenerator/CEvNS - - 7 - genie::EventGenerator/GLRES-Mu - genie::EventGenerator/GLRES-Tau - genie::EventGenerator/GLRES-Ele - genie::EventGenerator/GLRES-Had - genie::EventGenerator/DIS-CC - genie::EventGenerator/DIS-NC - genie::EventGenerator/DIS-CC-CHARM - - - - 4 - genie::EventGenerator/GLRES-Mu - genie::EventGenerator/GLRES-Tau - genie::EventGenerator/GLRES-Ele - genie::EventGenerator/GLRES-Had - - - - 1 - genie::EventGenerator/GLRES-Mu - - - - 1 - genie::EventGenerator/GLRES-Tau - - - - 1 - genie::EventGenerator/GLRES-Ele - - - - 1 - genie::EventGenerator/GLRES-Had - - 2 genie::EventGenerator/DFR-CC @@ -455,7 +416,39 @@ Generator-%d alg No - + + + + 13 + genie::EventGenerator/HEDIS-CC + genie::EventGenerator/HEDIS-NC + genie::EventGenerator/GLRES-Mu + genie::EventGenerator/GLRES-Tau + genie::EventGenerator/GLRES-Ele + genie::EventGenerator/GLRES-Had + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Had + genie::EventGenerator/PhotonCOH + + + + 11 + genie::EventGenerator/GLRES-Mu + genie::EventGenerator/GLRES-Tau + genie::EventGenerator/GLRES-Ele + genie::EventGenerator/GLRES-Had + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Had + genie::EventGenerator/PhotonCOH + 6 @@ -492,5 +485,82 @@ Generator-%d alg No genie::EventGenerator/HEDIS-NC + + 4 + genie::EventGenerator/GLRES-Mu + genie::EventGenerator/GLRES-Tau + genie::EventGenerator/GLRES-Ele + genie::EventGenerator/GLRES-Had + + + + 1 + genie::EventGenerator/GLRES-Mu + + + + 1 + genie::EventGenerator/GLRES-Tau + + + + 1 + genie::EventGenerator/GLRES-Ele + + + + 1 + genie::EventGenerator/GLRES-Had + + + + 2 + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + + + + 1 + genie::EventGenerator/HENuEl-CC + + + + 1 + genie::EventGenerator/HENuEl-NC + + + + 4 + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Had + + + + 1 + genie::EventGenerator/PhotonRES-Mu + + + + 1 + genie::EventGenerator/PhotonRES-Tau + + + + 1 + genie::EventGenerator/PhotonRES-Ele + + + + 1 + genie::EventGenerator/PhotonRES-Had + + + + 1 + genie::EventGenerator/PhotonCOH + + diff --git a/config/G00_00a/EventGenerator.xml b/config/G00_00a/EventGenerator.xml index 6f0a5af21..2d873f1ad 100644 --- a/config/G00_00a/EventGenerator.xml +++ b/config/G00_00a/EventGenerator.xml @@ -345,15 +345,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/G00_00a/ModelConfiguration.xml b/config/G00_00a/ModelConfiguration.xml index 5aafbaeba..106bca2db 100644 --- a/config/G00_00a/ModelConfiguration.xml +++ b/config/G00_00a/ModelConfiguration.xml @@ -131,7 +131,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G00_00b/EventGenerator.xml b/config/G00_00b/EventGenerator.xml index ffdf7e6f5..4baf4d531 100644 --- a/config/G00_00b/EventGenerator.xml +++ b/config/G00_00b/EventGenerator.xml @@ -345,15 +345,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/G00_00b/ModelConfiguration.xml b/config/G00_00b/ModelConfiguration.xml index 74c015019..d11e03e9b 100644 --- a/config/G00_00b/ModelConfiguration.xml +++ b/config/G00_00b/ModelConfiguration.xml @@ -137,7 +137,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G00_00z/EventGenerator.xml b/config/G00_00z/EventGenerator.xml index 585823ae5..56f1d07f1 100644 --- a/config/G00_00z/EventGenerator.xml +++ b/config/G00_00z/EventGenerator.xml @@ -345,15 +345,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/G00_00z/ModelConfiguration.xml b/config/G00_00z/ModelConfiguration.xml index b7fbcd64e..6be7ec112 100644 --- a/config/G00_00z/ModelConfiguration.xml +++ b/config/G00_00z/ModelConfiguration.xml @@ -140,7 +140,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_01a/ModelConfiguration.xml b/config/G18_01a/ModelConfiguration.xml index 5cafa550e..732c724bf 100644 --- a/config/G18_01a/ModelConfiguration.xml +++ b/config/G18_01a/ModelConfiguration.xml @@ -131,7 +131,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_01b/ModelConfiguration.xml b/config/G18_01b/ModelConfiguration.xml index 0f6632dc9..6bbc38648 100644 --- a/config/G18_01b/ModelConfiguration.xml +++ b/config/G18_01b/ModelConfiguration.xml @@ -131,7 +131,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_02a/ModelConfiguration.xml b/config/G18_02a/ModelConfiguration.xml index 996865804..c77a4c844 100644 --- a/config/G18_02a/ModelConfiguration.xml +++ b/config/G18_02a/ModelConfiguration.xml @@ -134,7 +134,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_02b/ModelConfiguration.xml b/config/G18_02b/ModelConfiguration.xml index 1d98d4303..113ab225c 100644 --- a/config/G18_02b/ModelConfiguration.xml +++ b/config/G18_02b/ModelConfiguration.xml @@ -134,7 +134,6 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_10a/ModelConfiguration.xml b/config/G18_10a/ModelConfiguration.xml index bf310de04..95c99abbf 100644 --- a/config/G18_10a/ModelConfiguration.xml +++ b/config/G18_10a/ModelConfiguration.xml @@ -135,8 +135,6 @@ STFC, Rutherford Appleton Laboratory genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_10b/ModelConfiguration.xml b/config/G18_10b/ModelConfiguration.xml index bfddc32af..3c4d400a6 100644 --- a/config/G18_10b/ModelConfiguration.xml +++ b/config/G18_10b/ModelConfiguration.xml @@ -136,8 +136,6 @@ STFC, Rutherford Appleton Laboratory genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_10i/ModelConfiguration.xml b/config/G18_10i/ModelConfiguration.xml index 68697e587..5e91814f9 100644 --- a/config/G18_10i/ModelConfiguration.xml +++ b/config/G18_10i/ModelConfiguration.xml @@ -136,8 +136,6 @@ STFC, Rutherford Appleton Laboratory genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G18_10j/ModelConfiguration.xml b/config/G18_10j/ModelConfiguration.xml index afe5729bb..217c1e4da 100644 --- a/config/G18_10j/ModelConfiguration.xml +++ b/config/G18_10j/ModelConfiguration.xml @@ -137,7 +137,6 @@ STFC, Rutherford Appleton Laboratory genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G21_11a/EventGenerator.xml b/config/G21_11a/EventGenerator.xml index 3d8152dcc..ae179f857 100644 --- a/config/G21_11a/EventGenerator.xml +++ b/config/G21_11a/EventGenerator.xml @@ -426,15 +426,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/G21_11a/ModelConfiguration.xml b/config/G21_11a/ModelConfiguration.xml index 3a9df23d0..b92594217 100644 --- a/config/G21_11a/ModelConfiguration.xml +++ b/config/G21_11a/ModelConfiguration.xml @@ -95,7 +95,6 @@ STFC, Rutherford Appleton Laboratory genie::EmpiricalMECPXSec2015/Default genie::SuSAv2MECPXSec/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/G21_11b/EventGenerator.xml b/config/G21_11b/EventGenerator.xml index 3d8152dcc..ae179f857 100644 --- a/config/G21_11b/EventGenerator.xml +++ b/config/G21_11b/EventGenerator.xml @@ -426,15 +426,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/G21_11b/ModelConfiguration.xml b/config/G21_11b/ModelConfiguration.xml index 6e5656ddc..eba4f6b8d 100644 --- a/config/G21_11b/ModelConfiguration.xml +++ b/config/G21_11b/ModelConfiguration.xml @@ -95,7 +95,6 @@ STFC, Rutherford Appleton Laboratory genie::EmpiricalMECPXSec2015/Default genie::SuSAv2MECPXSec/Default - genie::GLRESPXSec/Default genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/GDM18_00a/ModelConfiguration.xml b/config/GDM18_00a/ModelConfiguration.xml index 34a6b45be..5e303b804 100644 --- a/config/GDM18_00a/ModelConfiguration.xml +++ b/config/GDM18_00a/ModelConfiguration.xml @@ -131,7 +131,7 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default + genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/GDM18_00b/ModelConfiguration.xml b/config/GDM18_00b/ModelConfiguration.xml index 947ea4060..cf54a261d 100644 --- a/config/GDM18_00b/ModelConfiguration.xml +++ b/config/GDM18_00b/ModelConfiguration.xml @@ -131,7 +131,7 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default + genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/GEM21_11a/EventGenerator.xml b/config/GEM21_11a/EventGenerator.xml index 3d8152dcc..ae179f857 100644 --- a/config/GEM21_11a/EventGenerator.xml +++ b/config/GEM21_11a/EventGenerator.xml @@ -426,15 +426,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/GEM21_11b/EventGenerator.xml b/config/GEM21_11b/EventGenerator.xml index 3d8152dcc..ae179f857 100644 --- a/config/GEM21_11b/EventGenerator.xml +++ b/config/GEM21_11b/EventGenerator.xml @@ -426,15 +426,6 @@ XSecModel alg Yes Cross section model used at the thread genie::AMNuGammaInteractionListGenerator/Default - - - 3 - genie::InitialStateAppender/Default - genie::GLRESGenerator/Default - genie::UnstableParticleDecayer/AfterHadronTransport - genie::GLRESInteractionListGenerator/Default - - 5 diff --git a/config/GHE19_00a/CommonDecay.xml b/config/GHE19_00a/CommonDecay.xml index 7969bfcc6..d76b12158 100644 --- a/config/GHE19_00a/CommonDecay.xml +++ b/config/GHE19_00a/CommonDecay.xml @@ -74,6 +74,7 @@ University of Liverpool false false false + false diff --git a/config/GHE19_00a/ModelConfiguration.xml b/config/GHE19_00a/ModelConfiguration.xml index 57c2f2859..337467736 100644 --- a/config/GHE19_00a/ModelConfiguration.xml +++ b/config/GHE19_00a/ModelConfiguration.xml @@ -26,6 +26,13 @@ STFC, Rutherford Appleton Laboratory genie::GLRESPXSec/Default genie::GLRESPXSec/Default genie::GLRESPXSec/Default + genie::HENuElPXSec/Default + genie::HENuElPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonCOHPXSec/Default genie::HEDISPXSec/Default genie::HEDISPXSec/Default diff --git a/config/GHE19_00a/TuneGeneratorList.xml b/config/GHE19_00a/TuneGeneratorList.xml index d3b0947d3..cd8e39f64 100644 --- a/config/GHE19_00a/TuneGeneratorList.xml +++ b/config/GHE19_00a/TuneGeneratorList.xml @@ -47,13 +47,20 @@ Generator-%d alg No --> - 6 + 13 genie::EventGenerator/HEDIS-CC genie::EventGenerator/HEDIS-NC genie::EventGenerator/GLRES-Mu genie::EventGenerator/GLRES-Tau genie::EventGenerator/GLRES-Ele genie::EventGenerator/GLRES-Had + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Had + genie::EventGenerator/PhotonCOH diff --git a/config/GHE19_00b/CommonDecay.xml b/config/GHE19_00b/CommonDecay.xml index 7969bfcc6..d76b12158 100644 --- a/config/GHE19_00b/CommonDecay.xml +++ b/config/GHE19_00b/CommonDecay.xml @@ -74,6 +74,7 @@ University of Liverpool false false false + false diff --git a/config/GHE19_00b/CommonParam.xml b/config/GHE19_00b/CommonParam.xml index 13d82841b..1815315f7 100644 --- a/config/GHE19_00b/CommonParam.xml +++ b/config/GHE19_00b/CommonParam.xml @@ -14,7 +14,7 @@ 1 400 400 - 1e-6 + 1e-8 1. 1e9 80.398 diff --git a/config/GHE19_00b/ModelConfiguration.xml b/config/GHE19_00b/ModelConfiguration.xml index 57c2f2859..337467736 100644 --- a/config/GHE19_00b/ModelConfiguration.xml +++ b/config/GHE19_00b/ModelConfiguration.xml @@ -26,6 +26,13 @@ STFC, Rutherford Appleton Laboratory genie::GLRESPXSec/Default genie::GLRESPXSec/Default genie::GLRESPXSec/Default + genie::HENuElPXSec/Default + genie::HENuElPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonCOHPXSec/Default genie::HEDISPXSec/Default genie::HEDISPXSec/Default diff --git a/config/GHE19_00b/TuneGeneratorList.xml b/config/GHE19_00b/TuneGeneratorList.xml index d3b0947d3..cd8e39f64 100644 --- a/config/GHE19_00b/TuneGeneratorList.xml +++ b/config/GHE19_00b/TuneGeneratorList.xml @@ -47,13 +47,20 @@ Generator-%d alg No --> - 6 + 13 genie::EventGenerator/HEDIS-CC genie::EventGenerator/HEDIS-NC genie::EventGenerator/GLRES-Mu genie::EventGenerator/GLRES-Tau genie::EventGenerator/GLRES-Ele genie::EventGenerator/GLRES-Had + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Had + genie::EventGenerator/PhotonCOH diff --git a/config/GHE19_00c/CommonDecay.xml b/config/GHE19_00c/CommonDecay.xml index 7969bfcc6..d76b12158 100644 --- a/config/GHE19_00c/CommonDecay.xml +++ b/config/GHE19_00c/CommonDecay.xml @@ -74,6 +74,7 @@ University of Liverpool false false false + false diff --git a/config/GHE19_00c/ModelConfiguration.xml b/config/GHE19_00c/ModelConfiguration.xml index 57c2f2859..337467736 100644 --- a/config/GHE19_00c/ModelConfiguration.xml +++ b/config/GHE19_00c/ModelConfiguration.xml @@ -26,6 +26,13 @@ STFC, Rutherford Appleton Laboratory genie::GLRESPXSec/Default genie::GLRESPXSec/Default genie::GLRESPXSec/Default + genie::HENuElPXSec/Default + genie::HENuElPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonCOHPXSec/Default genie::HEDISPXSec/Default genie::HEDISPXSec/Default diff --git a/config/GHE19_00c/TuneGeneratorList.xml b/config/GHE19_00c/TuneGeneratorList.xml index d3b0947d3..cd8e39f64 100644 --- a/config/GHE19_00c/TuneGeneratorList.xml +++ b/config/GHE19_00c/TuneGeneratorList.xml @@ -47,13 +47,20 @@ Generator-%d alg No --> - 6 + 13 genie::EventGenerator/HEDIS-CC genie::EventGenerator/HEDIS-NC genie::EventGenerator/GLRES-Mu genie::EventGenerator/GLRES-Tau genie::EventGenerator/GLRES-Ele genie::EventGenerator/GLRES-Had + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Had + genie::EventGenerator/PhotonCOH diff --git a/config/GHE19_00d/CommonDecay.xml b/config/GHE19_00d/CommonDecay.xml new file mode 100644 index 000000000..d76b12158 --- /dev/null +++ b/config/GHE19_00d/CommonDecay.xml @@ -0,0 +1,81 @@ + + + + + + + + + + + true + true + true + + + + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + false + + + + diff --git a/config/GHE19_00d/CommonParam.xml b/config/GHE19_00d/CommonParam.xml new file mode 100644 index 000000000..e1567d87b --- /dev/null +++ b/config/GHE19_00d/CommonParam.xml @@ -0,0 +1,51 @@ + + + + + + + + + + NNPDF31sx_nlo_as_0118_LHCb_nf_6 + 0 + true + GGHR + 2 + 400 + 400 + 1e-9 + 2.6896 + 1e10 + 80.385 + 91.1876 + 0.00940161 + 0. + 0.97427 + 0.22536 + 0.00355 + 0.22522 + 0.97343 + 0.04140 + 0.00886 + 0.04050 + 0.99914 + + + + 2. + 100000 + 100000 + 1 + + + + 1.4 + + + + 100.000 + 1000000000000.000 + + + diff --git a/config/GHE19_00d/ModelConfiguration.xml b/config/GHE19_00d/ModelConfiguration.xml new file mode 100644 index 000000000..337467736 --- /dev/null +++ b/config/GHE19_00d/ModelConfiguration.xml @@ -0,0 +1,42 @@ + + + + + + + + + + + genie::GLRESPXSec/Default + genie::GLRESPXSec/Default + genie::GLRESPXSec/Default + genie::GLRESPXSec/Default + genie::HENuElPXSec/Default + genie::HENuElPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonRESPXSec/Default + genie::PhotonCOHPXSec/Default + genie::HEDISPXSec/Default + genie::HEDISPXSec/Default + + + + + diff --git a/config/GHE19_00d/TuneGeneratorList.xml b/config/GHE19_00d/TuneGeneratorList.xml new file mode 100644 index 000000000..cd8e39f64 --- /dev/null +++ b/config/GHE19_00d/TuneGeneratorList.xml @@ -0,0 +1,67 @@ + + + + + + + + 13 + genie::EventGenerator/HEDIS-CC + genie::EventGenerator/HEDIS-NC + genie::EventGenerator/GLRES-Mu + genie::EventGenerator/GLRES-Tau + genie::EventGenerator/GLRES-Ele + genie::EventGenerator/GLRES-Had + genie::EventGenerator/HENuEl-CC + genie::EventGenerator/HENuEl-NC + genie::EventGenerator/PhotonRES-Mu + genie::EventGenerator/PhotonRES-Ele + genie::EventGenerator/PhotonRES-Tau + genie::EventGenerator/PhotonRES-Had + genie::EventGenerator/PhotonCOH + + + + diff --git a/config/GLRESGenerator.xml b/config/GLRESGenerator.xml index 3065b63f4..225e7a2d0 100644 --- a/config/GLRESGenerator.xml +++ b/config/GLRESGenerator.xml @@ -2,6 +2,10 @@ + + HEDIS-PYTHIA diff --git a/config/GLRESInteractionListGenerator.xml b/config/GLRESInteractionListGenerator.xml deleted file mode 100644 index def1d1ef1..000000000 --- a/config/GLRESInteractionListGenerator.xml +++ /dev/null @@ -1,31 +0,0 @@ - - - - - - - - true - - - - true - - - - true - - - - true - - - - diff --git a/config/GLRESPXSec.xml b/config/GLRESPXSec.xml index 6f70498c7..773990665 100644 --- a/config/GLRESPXSec.xml +++ b/config/GLRESPXSec.xml @@ -1,12 +1,14 @@ + + - - HEDIS-PYTHIA - genie::GLRESXSec/Default + HEDIS-PYTHIA + genie::HELeptonXSec/Default + diff --git a/config/GTEST18_02c/ModelConfiguration.xml b/config/GTEST18_02c/ModelConfiguration.xml index 63fc0024f..6ccae260f 100644 --- a/config/GTEST18_02c/ModelConfiguration.xml +++ b/config/GTEST18_02c/ModelConfiguration.xml @@ -134,7 +134,7 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default + genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/GTEST18_02d/ModelConfiguration.xml b/config/GTEST18_02d/ModelConfiguration.xml index ced0f5192..3d33cced7 100644 --- a/config/GTEST18_02d/ModelConfiguration.xml +++ b/config/GTEST18_02d/ModelConfiguration.xml @@ -134,7 +134,7 @@ STFC, Rutherford Appleton Laboratory --> genie::EmpiricalMECPXSec2015/Default genie::EmpiricalMECPXSec2015/Default - genie::GLRESPXSec/Default + genie::DummyPXSec/Default genie::NNBarOscDummyPXSec/Default diff --git a/config/HEDISXSec.xml b/config/HEDISXSec.xml index 4b38784dc..b251edcd6 100644 --- a/config/HEDISXSec.xml +++ b/config/HEDISXSec.xml @@ -9,10 +9,8 @@ Configuration for the HEDIS integration algorithm HEDIS-SF - adaptive - 5000000 - 100000 - 0.0001 + vegas + 200000 diff --git a/config/HELeptonInteractionListGenerator.xml b/config/HELeptonInteractionListGenerator.xml new file mode 100644 index 000000000..2a6a648eb --- /dev/null +++ b/config/HELeptonInteractionListGenerator.xml @@ -0,0 +1,59 @@ + + + + + + + + true + + + + true + + + + true + + + + true + + + + true + + + + true + + + + true + + + + true + + + + true + + + + true + + + + true + + + + diff --git a/config/GLRESKinematicsGenerator.xml b/config/HELeptonKinematicsGenerator.xml similarity index 64% rename from config/GLRESKinematicsGenerator.xml rename to config/HELeptonKinematicsGenerator.xml index 9c6275cf0..6d29106e4 100644 --- a/config/GLRESKinematicsGenerator.xml +++ b/config/HELeptonKinematicsGenerator.xml @@ -3,7 +3,7 @@ 1.05 + 0.000001 + 0.1 + 1.0 diff --git a/config/GLRESXSec.xml b/config/HELeptonXSec.xml similarity index 60% rename from config/GLRESXSec.xml rename to config/HELeptonXSec.xml index 5f7c2018a..5d07d0159 100644 --- a/config/GLRESXSec.xml +++ b/config/HELeptonXSec.xml @@ -3,20 +3,17 @@ - adaptive - 40000 - 0.001 + vegas + 500000 diff --git a/config/HENuElGenerator.xml b/config/HENuElGenerator.xml new file mode 100644 index 000000000..7801edb13 --- /dev/null +++ b/config/HENuElGenerator.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + diff --git a/config/HENuElPXSec.xml b/config/HENuElPXSec.xml new file mode 100644 index 000000000..e3a2c7afa --- /dev/null +++ b/config/HENuElPXSec.xml @@ -0,0 +1,14 @@ + + + + + + + + HEDIS-PYTHIA + genie::HELeptonXSec/Default + + + diff --git a/config/Messenger.xml b/config/Messenger.xml index 42b898f6b..95a3e472b 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -93,7 +93,7 @@ NOTICE WARN NOTICE - NOTICE + WARN NOTICE INFO INFO @@ -111,7 +111,7 @@ NOTICE NOTICE WARN - NOTICE + WARN NOTICE NOTICE WARN @@ -151,7 +151,7 @@ NOTICE WARN NOTICE - NOTICE + WARN INFO NOTICE INFO @@ -164,12 +164,12 @@ NOTICE NOTICE INFO - NOTICE - NOTICE + WARN + WARN WARN WARN WARN - INFO + WARN NOTICE WARN NOTICE @@ -204,19 +204,26 @@ WARN NOTICE NOTICE - NOTICE + WARN NOTICE - WARN + WARN + WARN WARN WARN - WARN + WARN + WARN + WARN + WARN + WARN + WARN WARN WARN WARN WARN + WARN ERROR - NOTICE + WARN NOTICE NOTICE NOTICE @@ -227,8 +234,11 @@ NOTICE INFO INFO - INFO + WARN INFO INFO + INFO + INFO + INFO diff --git a/config/Messenger_inuke_verbose.xml b/config/Messenger_inuke_verbose.xml index 4e741b006..c131b2c2d 100644 --- a/config/Messenger_inuke_verbose.xml +++ b/config/Messenger_inuke_verbose.xml @@ -180,6 +180,22 @@ NOTICE NOTICE NOTICE + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + ERROR NOTICE NOTICE @@ -193,5 +209,8 @@ INFO INFO INFO + INFO + INFO + INFO diff --git a/config/Messenger_laconic.xml b/config/Messenger_laconic.xml index a3845ccb6..ca600c5ea 100644 --- a/config/Messenger_laconic.xml +++ b/config/Messenger_laconic.xml @@ -199,6 +199,22 @@ WARN WARN WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN WARN WARN @@ -213,5 +229,8 @@ WARN WARN WARN + WARN + WARN + WARN diff --git a/config/Messenger_rambling.xml b/config/Messenger_rambling.xml index 40411a260..f2f5e4e92 100644 --- a/config/Messenger_rambling.xml +++ b/config/Messenger_rambling.xml @@ -197,6 +197,20 @@ NOTICE NOTICE NOTICE + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + WARN + ERROR NOTICE NOTICE @@ -212,5 +226,8 @@ INFO INFO INFO + INFO + INFO + INFO diff --git a/config/Messenger_whisper.xml b/config/Messenger_whisper.xml index 5e54805d0..a8450bb2c 100644 --- a/config/Messenger_whisper.xml +++ b/config/Messenger_whisper.xml @@ -196,6 +196,20 @@ FATAL FATAL FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL + FATAL FATAL FATAL @@ -211,5 +225,8 @@ FATAL FATAL FATAL + FATAL + FATAL + FATAL diff --git a/config/PhotonCOHGenerator.xml b/config/PhotonCOHGenerator.xml new file mode 100644 index 000000000..a5a8053b1 --- /dev/null +++ b/config/PhotonCOHGenerator.xml @@ -0,0 +1,14 @@ + + + + + + + + HEDIS-PYTHIA + + + + diff --git a/config/PhotonCOHPXSec.xml b/config/PhotonCOHPXSec.xml new file mode 100644 index 000000000..b49e48c0f --- /dev/null +++ b/config/PhotonCOHPXSec.xml @@ -0,0 +1,14 @@ + + + + + + + + HEDIS-PYTHIA + genie::HELeptonXSec/Default + + + diff --git a/config/PhotonRESGenerator.xml b/config/PhotonRESGenerator.xml new file mode 100644 index 000000000..8b2c2d092 --- /dev/null +++ b/config/PhotonRESGenerator.xml @@ -0,0 +1,14 @@ + + + + + + + + HEDIS-SF,HEDIS-PYTHIA + + + + diff --git a/config/PhotonRESPXSec.xml b/config/PhotonRESPXSec.xml new file mode 100644 index 000000000..d910fb7ae --- /dev/null +++ b/config/PhotonRESPXSec.xml @@ -0,0 +1,14 @@ + + + + + + + + HEDIS-SF,HEDIS-PYTHIA + genie::HELeptonXSec/Default + + + diff --git a/config/master_config.xml b/config/master_config.xml index 89001b065..65775c665 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -43,7 +43,7 @@ SKKinematicsGenerator.xml DMELKinematicsGenerator.xml DMDISKinematicsGenerator.xml - GLRESKinematicsGenerator.xml + HELeptonKinematicsGenerator.xml HEDISKinematicsGenerator.xml IBDHadronicSystemGenerator.xml QELHadronicSystemGenerator.xml @@ -59,6 +59,9 @@ AMNuGammaGenerator.xml MECGenerator.xml GLRESGenerator.xml + HENuElGenerator.xml + PhotonCOHGenerator.xml + PhotonRESGenerator.xml HEDISGenerator.xml CEvNSEventGenerator.xml COHDNuEventGenerator.xml @@ -79,7 +82,7 @@ COHInteractionListGenerator.xml AMNuGammaInteractionListGenerator.xml MECInteractionListGenerator.xml - GLRESInteractionListGenerator.xml + HELeptonInteractionListGenerator.xml DFRInteractionListGenerator.xml SKInteractionListGenerator.xml DMELInteractionListGenerator.xml @@ -160,7 +163,7 @@ DMElectronXSec.xml DMELXSec.xml DMDISXSec.xml - GLRESXSec.xml + HELeptonXSec.xml HEDISXSec.xml BostedChristyEMPXSec.xml @@ -193,6 +196,9 @@ NuElectronPXSec.xml DMElectronPXSec.xml GLRESPXSec.xml + HENuElPXSec.xml + PhotonCOHPXSec.xml + PhotonRESPXSec.xml BergerSehgalRESPXSec2014.xml KuzminLyubushkinNaumovRESPXSec2014.xml ReinDFRPXSec.xml diff --git a/src/Apps/Makefile b/src/Apps/Makefile index 32b839ba5..9b13d1645 100644 --- a/src/Apps/Makefile +++ b/src/Apps/Makefile @@ -19,20 +19,22 @@ GENIE_LIBS = $(shell $(GENIE)/src/scripts/setup/genie-config --libs) LIBRARIES := $(LIBRARIES) $(GENIE_LIBS) #LIBRARIES := $(LIBRARIES) $(CERN_LIBRARIES) $(GENIE_LIBS) -TGT_BASE = gevgen \ - gevgen_hadron \ - gevdump \ - gevpick \ - gevscan \ - gevcomp \ - gxscomp \ - gmkspl \ - gspladd \ - gspl2root \ - gntpc \ - gpdfcomp \ - gsfcomp \ - gmkhedissf \ +TGT_BASE = gevgen \ + gevgen_hadron \ + gevdump \ + gevpick \ + gevscan \ + gevcomp \ + gxscomp \ + gmkspl \ + gspladd \ + gspl2root \ + gntpc \ + gpdfcomp \ + gsfcomp \ + gmkhedissf \ + gcalchedisdiffxsec \ + gmkphotonsf \ gconfigdump ifeq ($(strip $(GOPT_ENABLE_FNAL)),YES) @@ -244,6 +246,18 @@ $(GENIE_BIN_PATH)/gsfcomp: gSFComp.o $(call find_libs,gsfcomp) @echo "** Building gsfcomp" $(LD) $(LDFLAGS) gSFComp.o $(LIBRARIES) -o $(GENIE_BIN_PATH)/gsfcomp +# App to calculate differential cross sections from HEDIS +# +$(GENIE_BIN_PATH)/gcalchedisdiffxsec: gCalcHEDISDiffXsec.o $(call find_libs,gcalchedisdiffxsec) + @echo "** Building gcalchedisdiffxsec" + $(LD) $(LDFLAGS) gCalcHEDISDiffXsec.o $(LIBRARIES) -o $(GENIE_BIN_PATH)/gcalchedisdiffxsec + +# App to create structure functions suitable for HELepton +# +$(GENIE_BIN_PATH)/gmkphotonsf: gMakePhotonStrucFunc.o $(call find_libs,gmkphotonsf) + @echo "** Building gmkphotonsf" + $(LD) $(LDFLAGS) gMakePhotonStrucFunc.o $(LIBRARIES) -o $(GENIE_BIN_PATH)/gmkphotonsf + # App to create structure functions suitable for HEDIS # $(GENIE_BIN_PATH)/gmkhedissf: gMakeHEDISStrucFunc.o $(call find_libs,gmkhedissf) diff --git a/src/Apps/gCalcHEDISDiffXsec.cxx b/src/Apps/gCalcHEDISDiffXsec.cxx new file mode 100644 index 000000000..be9cb69f2 --- /dev/null +++ b/src/Apps/gCalcHEDISDiffXsec.cxx @@ -0,0 +1,338 @@ +//____________________________________________________________________________ +/*! +\program gcalchedisdiffxsec +\brief GENIE utility program calculating differential cross sections from + HEDIS package. + Syntax : + gcalchedisdiffxsec -p nu -t tgt -o root_file + -x table_type + --tune genie_tune + --event-generator-list list_name + Note : + [] marks optional arguments. + <> marks a list of arguments out of which only one can be + selected at any given time. + Options : + -p + the neutrino pdg code + -t + the target pdg code (format: 10LZZZAAAI) + -x + path to ascii file with x values + -y + path to ascii file with y values + -e + path to ascii file with energy values + -o + output ROOT file name + --tune + Specifies a GENIE comprehensive neutrino interaction model tune. + --event-generator-list + List of event generators to load in event generation drivers. + [default: "Default"]. + --message-thresholds + Allows users to customize the message stream thresholds. + The thresholds are specified using an XML file. + See $GENIE/config/Messenger.xml for the XML schema. + *** See the User Manual for more details and examples. *** +\author Alfonso Garcia + NIKHEF (Amsterdam) +\created May 26, 2021 +\cpright Copyright (c) 2003-2021, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#include "TMath.h" +#include "TSystem.h" +#include "TTree.h" +#include "TFile.h" + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Framework/EventGen/InteractionList.h" +#include "Framework/EventGen/EventGeneratorI.h" +#include "Framework/EventGen/GEVGDriver.h" +#include "Framework/Utils/AppInit.h" +#include "Framework/Utils/RunOpt.h" +#include "Framework/Utils/KineUtils.h" +#include "Framework/Utils/CmdLnArgParser.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Interaction/InitialState.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/Conventions/Units.h" + +#include +#include + +using namespace genie; + +void PrintSyntax (void); +void DecodeCommandLine (int argc, char * argv[]); + +vector ReadListFromPath (string path); + +const double epsilon = 1e-5; + +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +//INPUT ARGUMENTS +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- + +int fNu = -1; +int fTgt = -1; +string fPathXlist = ""; +string fPathYlist = ""; +string fPathElist = ""; +string fOutFileName = ""; +bool fSaveAll = false; + +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +//MAIN PROGRAM +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- + +int main(int argc, char ** argv) { + + DecodeCommandLine(argc,argv); + + RunOpt::Instance()->BuildTune(); + utils::app_init::MesgThresholds(RunOpt::Instance()->MesgThresholdFiles()); + + string fChannel = RunOpt::Instance()->EventGeneratorList(); + + GEVGDriver evg_driver; + InitialState init_state(fTgt, fNu); + evg_driver.SetEventGeneratorList(fChannel); + evg_driver.Configure(init_state); + + vector ve = ReadListFromPath(fPathElist); + vector vx = ReadListFromPath(fPathXlist); + vector vy = ReadListFromPath(fPathYlist); + + double widthx = (TMath::Log10(vx[1])-TMath::Log10(vx[0])); + + LOG("gcalchedisdiffxsec", pINFO) << widthx; + + int quark; + double ei; + double bx; + double by; + double dxsec; + + TString treename = Form("diffxsec_nu%d_tgt%d_%s",fNu,fTgt,fChannel.c_str()); + + LOG("gcalchedisdiffxsec", pINFO) << treename; + + TTree * tree = new TTree(treename,treename); + tree->Branch( "Quark", &quark, "Quark/I" ); + tree->Branch( "Ei", &ei, "Ei/D" ); + tree->Branch( "Bx", &bx, "Bx/D" ); + tree->Branch( "By", &by, "By/D" ); + tree->Branch( "DiffXsec", &dxsec, "DiffXsec/D" ); + + const InteractionList * intlst = evg_driver.Interactions(); + + InteractionList::const_iterator intliter; + for(intliter = intlst->begin(); intliter != intlst->end(); ++intliter) { + + const Interaction * interaction = *intliter; + + quark = 10000 * interaction->InitState().Tgt().HitQrkPdg(); + if (!interaction->InitState().Tgt().HitSeaQrk()) { + if ( quark>0 ) quark += 100; + else quark -= 100; + } + quark += 1 * interaction->ExclTag().FinalQuarkPdg(); + + LOG("gcalchedisdiffxsec", pINFO) << "Current interaction: " << interaction->AsString(); + + const XSecAlgorithmI * xsec_alg = evg_driver.FindGenerator(interaction)->CrossSectionAlg(); + + for ( unsigned i=0; iInitStatePtr()->SetProbeP4(p4); + + for ( unsigned j=0; jby_prev) by = 1e-4; + } + + if ( by<0 || by>1 ) continue; + + interaction->KinePtr()->Sety(by); + + if (fSaveAll) { + for ( unsigned k=0; kKinePtr()->Setx(bx); + utils::kinematics::UpdateWQ2FromXY(interaction); + dxsec = xsec_alg->XSec(interaction, kPSxyfE) / units::cm2; + LOG("gcalchedisdiffxsec", pDEBUG) << "x: " << bx << " y: " << by << " -> d2sigmadxdy[E=" << ei << "GeV] = " << dxsec << " cm2"; + tree->Fill(); + } + } + else { + dxsec = 0.; + for ( unsigned k=0; kKinePtr()->Setx(bx); + utils::kinematics::UpdateWQ2FromXY(interaction); + dxsec += xsec_alg->XSec(interaction, kPSxyfE) * bx; + } + dxsec *= widthx*TMath::Log(10) / units::cm2; + LOG("gcalchedisdiffxsec", pDEBUG) << " y: " << by << " -> d2sigmady[E=" << ei << "GeV] = " << dxsec << " cm2"; + tree->Fill(); + } + + } + + } + + } + + TFile * outfile = new TFile(fOutFileName.c_str(),"RECREATE"); + tree->Write(tree->GetName()); + outfile->Close(); + +} + + +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +//READ LIST +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +vector ReadListFromPath(string path) { + + vector list; + + std::ifstream infile(path.c_str()); + + //check to see that the file was opened correctly: + if (!infile.is_open()) { + LOG("gcalchedisdiffxsec", pFATAL) << "There was a problem opening the input file!"; + exit(1);//exit or do additional error checking + } + + double val = 0.; + while (infile >> val) list.push_back(val); + + return list; + +} + + +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- +//INPUT PARSER +//---------------------------------------------------------------------------------------- +//---------------------------------------------------------------------------------------- + +void DecodeCommandLine(int argc, char * argv[]) { + + // Common run options. Set defaults and read. + RunOpt::Instance()->ReadFromCommandLine(argc,argv); + + // Parse run options for this app + CmdLnArgParser parser(argc,argv); + + if( parser.OptionExists('p') ){ + fNu = parser.ArgAsInt('p'); + LOG("gcalchedisdiffxsec", pINFO) << "Probe = " << fNu; + } + else { + LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input neutrino type!"; + PrintSyntax(); + exit(1); + } + + if( parser.OptionExists('t') ){ + fTgt = parser.ArgAsInt('t'); + LOG("gcalchedisdiffxsec", pINFO) << "Target = " << fTgt; + } + else { + LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input target type!"; + PrintSyntax(); + exit(1); + } + + if( parser.OptionExists('x') ){ + fPathXlist = parser.ArgAsString('x'); + LOG("gcalchedisdiffxsec", pINFO) << "PathXlist = " << fPathXlist; + } + else { + LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input path to Xlist!"; + PrintSyntax(); + exit(1); + } + + if( parser.OptionExists('y') ){ + fPathYlist = parser.ArgAsString('y'); + LOG("gcalchedisdiffxsec", pINFO) << "PathYlist = " << fPathYlist; + } + else { + LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input path to Ylist!"; + PrintSyntax(); + exit(1); + } + + if( parser.OptionExists('e') ){ + fPathElist = parser.ArgAsString('e'); + LOG("gcalchedisdiffxsec", pINFO) << "PathElist = " << fPathElist; + } + else { + LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified input path to Elist!"; + PrintSyntax(); + exit(1); + } + + if( parser.OptionExists('s') ) { + fSaveAll = true; + LOG("gcalchedisdiffxsec", pINFO) << "SaveAll = " << fSaveAll; + } + + if( parser.OptionExists('o') ){ + fOutFileName = parser.ArgAsString('o'); + LOG("gcalchedisdiffxsec", pINFO) << "OutFileName = " << fOutFileName; + } + else { + LOG("gcalchedisdiffxsec", pFATAL) << "Unspecified output file name!"; + PrintSyntax(); + exit(1); + } + +} + +//____________________________________________________________________________ +void PrintSyntax(void) +{ + LOG("gcalchedisdiffxsec", pNOTICE) + << "\n\n" << "Syntax:" << "\n" + << " gcalchedisdiffxsec -p nu -t tgt -o root_file\n" + << " -x pathxlist\n" + << " -y pathylist\n" + << " -e pathelist\n" + << " --tune genie_tune\n" + << " --event-generator-list list_name\n" + << " [--message-thresholds xml_file]\n"; +} \ No newline at end of file diff --git a/src/Apps/gMakePhotonStrucFunc.cxx b/src/Apps/gMakePhotonStrucFunc.cxx new file mode 100644 index 000000000..dd030d9f1 --- /dev/null +++ b/src/Apps/gMakePhotonStrucFunc.cxx @@ -0,0 +1,197 @@ +//____________________________________________________________________________ +/*! +\program gmkspl +\brief GENIE utility program building Structure Functions needed for HEDIS + package. + Syntax : + gmkphotonsf [-h] + Note : + [] marks optional arguments. + <> marks a list of arguments out of which only one can be + selected at any given time. + Options : + *** See the User Manual for more details and examples. *** +\author Alfonso Garcia + Harvard University & IFIC +\created Dev 8, 2020 +\cpright Copyright (c) 2003-2019, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#include "Framework/Messenger/Messenger.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" + +#include +#include +#include +#include + +#include +#include +#include + +#ifdef __GENIE_APFEL_ENABLED__ +#include "APFEL/APFEL.h" +#endif +#include "LHAPDF/LHAPDF.h" +#ifdef __GENIE_LHAPDF6_ENABLED__ +const LHAPDF::PDF* pdf_cache; +#endif + +using namespace std; + +using namespace genie; +using namespace genie::constants; + +int fNucPdg = 0; + +#ifdef __GENIE_APFEL_ENABLED__ +class PhotonConv: public ROOT::Math::IBaseFunctionOneDim +{ + public: + PhotonConv(double x, double q) : ROOT::Math::IBaseFunctionOneDim(),xmin(x),Qin(q){}; + ~PhotonConv(){}; + ROOT::Math::IBaseFunctionOneDim * Clone (void) const { return new PhotonConv(xmin,Qin); } + unsigned int NDim (void) const { return 1; } + double DoEval (double y) const { return 2. * ( TMath::Power(xmin/y,2)+TMath::Power( 1.-xmin/y, 2) ) * pdf_cache->xfxQ(22, y, Qin); } + void SetPar (double x, double q) { xmin=x; Qin=q; } + + private: + double xmin,Qin; +}; + +extern "C" void externalsetapfellept_(double* x, double* q, int* irep, double* xl, double* xf); + + +// Requires a global cache of pdf_cache LHAPDF::PDF +void externalsetapfellept_(double* x, double* q, int* irep, double* xl, double* xf){ + if (*x >= 1 || *x < 0) { + for ( int i=0; i<13; i++ ) xf[i] = 0.; + for ( int i=0; i<7; i++ ) xl[i] = 0.; + return; + } + else{ + for ( int i=0; i<13; i++ ) { + xf[i] = pdf_cache->xfxQ(i-6, *x, *q); + if ( pdg::IsNeutron(fNucPdg) ) { + if ( i==4 ) xf[i] = pdf_cache->xfxQ(-1, *x, *q); + else if( i==5 ) xf[i] = pdf_cache->xfxQ(-2, *x, *q); + else if( i==7 ) xf[i] = pdf_cache->xfxQ( 2, *x, *q); + else if( i==8 ) xf[i] = pdf_cache->xfxQ( 1, *x, *q); + } + } + for ( int i=0; i<7; i++ ) { + if ( pdg::IsProton (fNucPdg) ) { + if (i==0 || i==6) xl[i] = 0.; + else if (i==3 ) xl[i] = pdf_cache->xfxQ(22, *x, *q); + else{ + double mlep; + if (i==1 || i==5) mlep = kMuonMass; + else if (i==2 || i==4) mlep = kElectronMass; + ROOT::Math::IBaseFunctionOneDim * func = new PhotonConv(*x,*q); + ROOT::Math::IntegrationOneDim::Type ig_type = ROOT::Math::IntegrationOneDim::kADAPTIVE; + ROOT::Math::Integrator ig(*func,ig_type,1,0.01,100000); + double res = ig.Integral(*x,1.); + xl[i] = APFEL::AlphaQED(*q)/2./kPi*TMath::Log( *q/mlep ) * res; + } + } + else if ( pdg::IsNeutron(fNucPdg) ) xl[i] = 0.0; + } + } + + return; +} +#endif + +//____________________________________________________________________________ +int main(int argc, char ** argv) +{ + + const int nx = 1000.; + + int nucs[2] = { kPdgProton, kPdgNeutron }; + int pdgs[6] = { kPdgNuE, kPdgAntiNuE, kPdgNuMu, kPdgAntiNuMu, kPdgNuTau, kPdgAntiNuTau }; + +#ifdef __GENIE_APFEL_ENABLED__ + for (int k=0; k<2; k++) { + + fNucPdg = nucs[k]; + + // initialising APFEL framework + LOG("gmkphotonsf", pINFO) << "Initialising APFEL..." ; + string pdfset; + if (pdg::IsProton (fNucPdg) ) pdfset = "NNPDF31_nnlo_as_0118_luxqed"; + else if (pdg::IsNeutron(fNucPdg) ) pdfset = "NNPDF31_nnlo_as_0118"; + + const LHAPDF::PDFSet set(pdfset); + pdf_cache = LHAPDF::mkPDF(pdfset,0); + + double xPDFmin,QPDFmin,QPDFmax,mc,mb,mt; + stringstream( set.get_entry("MCharm") ) >> mc; + stringstream( set.get_entry("MBottom") ) >> mb; + stringstream( set.get_entry("MTop") ) >> mt; + stringstream( set.get_entry("QMin") ) >> QPDFmin; + stringstream( set.get_entry("QMax") ) >> QPDFmax; + stringstream( set.get_entry("XMin") ) >> xPDFmin; + LOG("gmkphotonsf", pINFO) << "xPDFmin = " << xPDFmin; + LOG("gmkphotonsf", pINFO) << "QPDFmin = " << QPDFmin; + LOG("gmkphotonsf", pINFO) << "QPDFmax = " << QPDFmax; + LOG("gmkphotonsf", pINFO) << "mc = " << mc; + LOG("gmkphotonsf", pINFO) << "mb = " << mb; + LOG("gmkphotonsf", pINFO) << "mt = " << mt; + + APFEL::CleanUp(); + + APFEL::SetPDFSet(pdfset); + APFEL::SetReplica(0); + APFEL::SetPerturbativeOrder(2); + APFEL::SetQLimits(QPDFmin,QPDFmax); + APFEL::SetMaxFlavourPDFs(6); + APFEL::SetMaxFlavourAlpha(6); + APFEL::SetNumberOfGrids(3); + APFEL::SetGridParameters(1,100,5,1e-9); + APFEL::SetGridParameters(2,40,3,1e-1); + APFEL::SetGridParameters(3,60,3,8e-1); + APFEL::SetGFermi(kGF); + APFEL::SetPoleMasses(mc,mb,mt); + APFEL::SetTheory("QUniD"); + APFEL::EnableLeptonEvolution(true); + APFEL::SetFastEvolution(true); + APFEL::SetPDFSet("leptexternal"); + APFEL::InitializeAPFEL(); + + LOG("gmkphotonsf", pWARN) << "Init EvolveAPFEL"; + APFEL::EvolveAPFEL(QPDFmin,kMw); + LOG("gmkphotonsf", pWARN) << "End EvolveAPFEL"; + + // open file in which SF will be stored + + double x[nx]; + for ( int i=0; iGetenv("GENIE")) + "/data/evgen/photon-sf/PhotonSF_hitnuc"+to_string(fNucPdg)+"_hitlep"+to_string(pdgs[j])+".dat"; + std::ofstream sf_stream(SFname); + for ( int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); } } if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) { for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); } } if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) { for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); } } if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) { for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); } } } + TGraph * gr_discc = new TGraph(kNSplineP, e, xsdiscc); + gr_discc->SetName("dis_cc"); + FormatXSecGraph(gr_discc); + topdir->Add(gr_discc); + TGraph * gr_disnc = new TGraph(kNSplineP, e, xsdisnc); + gr_disnc->SetName("dis_nc"); + FormatXSecGraph(gr_disnc); + topdir->Add(gr_disnc); TGraph * gr_disccp = new TGraph(kNSplineP, e, xsdisccp); gr_disccp->SetName("dis_cc_p"); FormatXSecGraph(gr_disccp); @@ -857,43 +875,57 @@ void SaveGraphsToRootFile(void) // for(int i=0; ibegin(); ilistiter != ilist->end(); ++ilistiter) { - const Interaction * interaction = *ilistiter; - const ProcessInfo & proc = interaction->ProcInfo(); - const XclsTag & xcls = interaction->ExclTag(); - const InitialState & init = interaction->InitState(); - const Target & tgt = init.Tgt(); + const Interaction * interaction = *ilistiter; + const ProcessInfo & proc = interaction->ProcInfo(); + const XclsTag & xcls = interaction->ExclTag(); + const InitialState & init = interaction->InitState(); + const Target & tgt = init.Tgt(); - const Spline * spl = evg_driver.XSecSpline(interaction); + const Spline * spl = evg_driver.XSecSpline(interaction); - if(!xcls.IsCharmEvent()) continue; + if(!xcls.IsCharmEvent()) continue; - if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) { - for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); - } - } - if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) { - for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); - } - } - if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) { - for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); - } - } - if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) { - for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); - } - } + if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsProton(tgt.HitNucPdg())) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + xsdisccp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); + } + } + if (proc.IsDeepInelastic() && proc.IsWeakCC() && pdg::IsNeutron(tgt.HitNucPdg())) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + xsdisccn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); + } + } + if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsProton(tgt.HitNucPdg())) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + xsdisncp[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); + } + } + if (proc.IsDeepInelastic() && proc.IsWeakNC() && pdg::IsNeutron(tgt.HitNucPdg())) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + xsdisncn[i] += (spl->Evaluate(e[i]) * (1E+38/units::cm2)); + } + } } + TGraph * gr_discc_charm = new TGraph(kNSplineP, e, xsdiscc); + gr_discc_charm->SetName("dis_cc_charm"); + FormatXSecGraph(gr_discc_charm); + topdir->Add(gr_discc_charm); + TGraph * gr_disnc_charm = new TGraph(kNSplineP, e, xsdisnc); + gr_disnc_charm->SetName("dis_nc_charm"); + FormatXSecGraph(gr_disnc_charm); + topdir->Add(gr_disnc_charm); TGraph * gr_disccp_charm = new TGraph(kNSplineP, e, xsdisccp); gr_disccp_charm->SetName("dis_cc_p_charm"); FormatXSecGraph(gr_disccp_charm); @@ -1002,6 +1034,55 @@ void SaveGraphsToRootFile(void) FormatXSecGraph(gr_cohtot); topdir->Add(gr_cohtot); + // + // add-up all glres and photon-res channels + // + + double * xsglrescc = new double[kNSplineP]; + double * xsglresnc = new double[kNSplineP]; + double * xsphrescc = new double[kNSplineP]; + for(int i=0; ibegin(); ilistiter != ilist->end(); ++ilistiter) { + const Interaction * interaction = *ilistiter; + const ProcessInfo & proc = interaction->ProcInfo(); + + const Spline * spl = evg_driver.XSecSpline(interaction); + + if (proc.IsGlashowResonance() && proc.IsWeakCC()) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + } + } + if (proc.IsGlashowResonance() && proc.IsWeakNC()) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + } + } + if (proc.IsPhotonResonance()) { + for(int i=0; iEvaluate(e[i]) * (1E+38/units::cm2)); + } + + } + } + TGraph * gr_glrescc = new TGraph(kNSplineP, e, xsglrescc); + gr_glrescc->SetName("glres_cc"); + FormatXSecGraph(gr_glrescc); + topdir->Add(gr_glrescc); + TGraph * gr_glresnc = new TGraph(kNSplineP, e, xsglresnc); + gr_glresnc->SetName("glres_nc"); + FormatXSecGraph(gr_glresnc); + topdir->Add(gr_glresnc); + TGraph * gr_phrescc = new TGraph(kNSplineP, e, xsphrescc); + gr_phrescc->SetName("phres_cc"); + FormatXSecGraph(gr_phrescc); + topdir->Add(gr_phrescc); + + // // total cross sections // @@ -1100,9 +1181,14 @@ void SaveGraphsToRootFile(void) delete [] xsdisccn; delete [] xsdisncp; delete [] xsdisncn; + delete [] xsdiscc; + delete [] xsdisnc; delete [] xscohcc; delete [] xscohnc; delete [] xscohtot; + delete [] xsglrescc; + delete [] xsglresnc; + delete [] xsphrescc; delete [] xstotcc; delete [] xstotccp; delete [] xstotccn; diff --git a/src/Framework/Conventions/KinePhaseSpace.h b/src/Framework/Conventions/KinePhaseSpace.h index 5c1a63e6e..9c3b714c8 100644 --- a/src/Framework/Conventions/KinePhaseSpace.h +++ b/src/Framework/Conventions/KinePhaseSpace.h @@ -70,7 +70,9 @@ typedef enum EKinePhaseSpace { kPSEgTlOgfE, kPSDMELEvGen, // Equivalent to kPSQELEvGen for Dark Matter scattering kPSlog10xlog10Q2fE, - kPSEDNufE // Used for Dark Neutrinos, two body final state + kPSEDNufE, // Used for Dark Neutrinos, two body final state + kPSn1n2fE, + kPSn1n2n3fE } KinePhaseSpace_t; class KinePhaseSpace @@ -128,6 +130,8 @@ class KinePhaseSpace case(kPSTAfE) : return "<{TA}|E>"; break; case(kPSlog10xlog10Q2fE) : return "<{log10x,log10Q2}|E>"; break; case(kPSEDNufE) : return "<{EDNu}|E>"; break; + case(kPSn1n2fE) : return "<{n1,n2}|E>"; break; + case(kPSn1n2n3fE) : return "<{n1,n2,n3}|E>"; break; } return "** Undefined kinematic phase space **"; } diff --git a/src/Framework/Conventions/KineVar.h b/src/Framework/Conventions/KineVar.h index 802e6719e..d79a43c3e 100644 --- a/src/Framework/Conventions/KineVar.h +++ b/src/Framework/Conventions/KineVar.h @@ -57,6 +57,9 @@ typedef enum EKineVar { kKVQ3, kKVSelQ0, kKVSelQ3, + kKVn1, + kKVn2, + kKVn3, // put all new enum names right before this line // do not change any previous ordering (neither insert nor delete) kNumOfKineVar @@ -101,6 +104,9 @@ class KineVar case(kKVQ3) : return " *Running* three momentum transfer" ; break; case(kKVSelQ0) : return "*Selected* energy transfer (Q0) " ; break; case(kKVSelQ3) : return "*Selected* three momentum transfer" ; break; + case(kKVn1) : return " *Running* Normalized variable n1" ; break; + case(kKVn2) : return " *Running* Normalized variable n2" ; break; + case(kKVn3) : return " *Running* Normalized variable n3" ; break; default : return "** Unknown kinematic variable **"; break; } diff --git a/src/Framework/Interaction/Interaction.cxx b/src/Framework/Interaction/Interaction.cxx index 16d8065b7..6b18f9f81 100644 --- a/src/Framework/Interaction/Interaction.cxx +++ b/src/Framework/Interaction/Interaction.cxx @@ -144,6 +144,9 @@ int Interaction::FSPrimLeptonPdg(void) const if (proc_info.IsNuElectronElastic()) return kPdgElectron; + if (proc_info.IsGlashowResonance() || proc_info.IsPhotonResonance()) + return xclstag.FinalLeptonPdg(); + // vN (Weak-NC) or eN (EM) if (proc_info.IsWeakNC() || proc_info.IsEM() || proc_info.IsWeakMix() || proc_info.IsDarkMatter()) return pdgc; // EDIT: DM does not change in FS @@ -152,12 +155,6 @@ int Interaction::FSPrimLeptonPdg(void) const int clpdgc; if (proc_info.IsIMDAnnihilation()) clpdgc = kPdgMuon; - else if (proc_info.IsGlashowResonance()) { - if ( pdg::IsMuon(xclstag.FinalLeptonPdg()) ) clpdgc = kPdgMuon; - else if ( pdg::IsTau(xclstag.FinalLeptonPdg()) ) clpdgc = kPdgTau; - else if ( pdg::IsElectron(xclstag.FinalLeptonPdg()) ) clpdgc = kPdgElectron; - else if ( pdg::IsPion(xclstag.FinalLeptonPdg()) ) clpdgc = kPdgPiP; - } else clpdgc = pdg::Neutrino2ChargedLepton(pdgc); return clpdgc; diff --git a/src/Framework/Interaction/KPhaseSpace.cxx b/src/Framework/Interaction/KPhaseSpace.cxx index 01e8756c6..04d76007c 100644 --- a/src/Framework/Interaction/KPhaseSpace.cxx +++ b/src/Framework/Interaction/KPhaseSpace.cxx @@ -202,6 +202,20 @@ double KPhaseSpace::Threshold(void) const double Ethr = 0.5 * (ml*ml-kElectronMass2)/kElectronMass; return TMath::Max(0.,Ethr); } + if(pi.IsPhotonResonance()) { + double Mn = tgt.HitNucP4Ptr()->M(); + double Ethr = 0.5 * (ml*ml-TMath::Power(Mn,2))/Mn; + return TMath::Max(0.,Ethr); + } + if(pi.IsPhotonCoherent()) { + double ml = 0; + if (pdg::IsNuE (TMath::Abs(init_state.ProbePdg()))) ml = kElectronMass; + else if (pdg::IsNuMu (TMath::Abs(init_state.ProbePdg()))) ml = kMuonMass; + else if (pdg::IsNuTau(TMath::Abs(init_state.ProbePdg()))) ml = kTauMass; + double MA = init_state.Tgt().Z()*kProtonMass + init_state.Tgt().N()*kNeutronMass; + double Ethr = 0.5 * (TMath::Power(kMw+ml,2)-TMath::Power(MA,2))/MA; + return TMath::Max(0.,Ethr); + } SLOG("KPhaseSpace", pERROR) @@ -261,6 +275,8 @@ bool KPhaseSpace::IsAboveThreshold(void) const pi.IsNuElectronElastic() || pi.IsDarkMatterElectronElastic() || pi.IsMEC() || + pi.IsPhotonCoherent() || + pi.IsPhotonResonance() || pi.IsGlashowResonance()) { E = init_state.ProbeE(kRfLab); @@ -330,7 +346,7 @@ bool KPhaseSpace::IsAllowed(void) const } //IMD - if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic() || pi.IsDarkMatterElectronElastic() || pi.IsGlashowResonance()) { + if(pi.IsInverseMuDecay() || pi.IsIMDAnnihilation() || pi.IsNuElectronElastic() || pi.IsDarkMatterElectronElastic()) { Range1D_t yl = this->YLim(); double y = kine.y(); bool in_phys = math::IsWithinLimits(y, yl); @@ -772,16 +788,6 @@ Range1D_t KPhaseSpace::YLim(void) const yl.max = 1. -ml/Ev - controls::kASmallNum; return yl; } - // GLRES - if(pi.IsGlashowResonance()) { - const InitialState & init_state = fInteraction->InitState(); - double Ev = init_state.ProbeE(kRfLab); - double ml = fInteraction->FSPrimLepton()->Mass(); - double me = kElectronMass; - yl.min = (ml*ml+me*me)/2/Ev/me + controls::kASmallNum; - yl.max = (4*Ev*(Ev+me) + (ml*ml+me*me))/2/Ev/(me+2*Ev) - controls::kASmallNum; - return yl; - } return yl; } //____________________________________________________________________________ diff --git a/src/Framework/Interaction/ProcessInfo.cxx b/src/Framework/Interaction/ProcessInfo.cxx index 9298237c1..2764e6d15 100644 --- a/src/Framework/Interaction/ProcessInfo.cxx +++ b/src/Framework/Interaction/ProcessInfo.cxx @@ -148,6 +148,16 @@ bool ProcessInfo::IsGlashowResonance(void) const return (fScatteringType == kScGlashowResonance); } //____________________________________________________________________________ +bool ProcessInfo::IsPhotonCoherent(void) const +{ + return (fScatteringType == kScPhotonCoherent); +} +//____________________________________________________________________________ +bool ProcessInfo::IsPhotonResonance(void) const +{ + return (fScatteringType == kScPhotonResonance); +} +//____________________________________________________________________________ bool ProcessInfo::IsAMNuGamma(void) const { return (fScatteringType == kScAMNuGamma); diff --git a/src/Framework/Interaction/ProcessInfo.h b/src/Framework/Interaction/ProcessInfo.h index 5af7040e9..5d4c50884 100644 --- a/src/Framework/Interaction/ProcessInfo.h +++ b/src/Framework/Interaction/ProcessInfo.h @@ -74,6 +74,8 @@ class ProcessInfo : public TObject { bool IsDarkMatterElectronElastic (void) const; bool IsInverseBetaDecay (void) const; bool IsGlashowResonance (void) const; + bool IsPhotonResonance (void) const; + bool IsPhotonCoherent (void) const; bool IsAMNuGamma (void) const; bool IsMEC (void) const; bool IsDiffractive (void) const; diff --git a/src/Framework/Interaction/ScatteringType.h b/src/Framework/Interaction/ScatteringType.h index 728407990..db4138ef2 100644 --- a/src/Framework/Interaction/ScatteringType.h +++ b/src/Framework/Interaction/ScatteringType.h @@ -50,6 +50,8 @@ typedef enum EScatteringType { kScInverseBetaDecay, kScGlashowResonance, kScIMDAnnihilation, + kScPhotonCoherent, + kScPhotonResonance, kScDarkMatterElastic = 101, kScDarkMatterDeepInelastic, kScDarkMatterElectron @@ -78,6 +80,8 @@ class ScatteringType case(kScInverseBetaDecay) : return "IBD"; break; case(kScGlashowResonance) : return "GLR"; break; case(kScIMDAnnihilation) : return "IMDAnh"; break; + case(kScPhotonCoherent) : return "PhotonCOH"; break; + case(kScPhotonResonance) : return "PhotonRES"; break; case(kScDarkMatterElastic) : return "DMEL"; break; case(kScDarkMatterDeepInelastic) : return "DMDIS"; break; case(kScDarkMatterElectron) : return "DME"; break; diff --git a/src/Framework/Numerical/MathUtils.h b/src/Framework/Numerical/MathUtils.h index cea200538..97e2454ea 100644 --- a/src/Framework/Numerical/MathUtils.h +++ b/src/Framework/Numerical/MathUtils.h @@ -80,16 +80,24 @@ namespace math } } - void Boost (long double bz) { - long double b2 = bz*bz; - long double gamma = 1.0 / sqrtl(1.0 - b2); - long double bp = bz*fPz; - long double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0; - fPz = fPz + gamma2*bp*bz + gamma*bz*fE; - fE = gamma*(fE + bp); - } - - + void BoostZ (long double bz) { + long double b2 = bz*bz; + long double gamma = 1.0 / sqrtl(1.0 - b2); + long double bp = bz*fPz; + long double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0; + fPz = fPz + gamma2*bp*bz + gamma*bz*fE; + fE = gamma*(fE + bp); + } + + void BoostY (long double by) { + long double b2 = by*by; + long double gamma = 1.0 / sqrtl(1.0 - b2); + long double bp = by*fPy; + long double gamma2 = b2 > 0 ? (gamma - 1.0)/b2 : 0.0; + fPy = fPy + gamma2*bp*by + gamma*by*fE; + fE = gamma*(fE + bp); + } + private : long double fPx; diff --git a/src/Framework/ParticleData/NaturalIsotopes.cxx b/src/Framework/ParticleData/NaturalIsotopes.cxx index fb1645fb2..b29309fcf 100644 --- a/src/Framework/ParticleData/NaturalIsotopes.cxx +++ b/src/Framework/ParticleData/NaturalIsotopes.cxx @@ -106,6 +106,29 @@ const NaturalIsotopeElementData * return vec[ielement]; } //____________________________________________________________________________ +const NaturalIsotopeElementData * + NaturalIsotopes::ElementDataPdg(int Z, int pdgcode) const +{ + map >::const_iterator miter; + + if( (miter=fNaturalIsotopesTable.find(Z)) == fNaturalIsotopesTable.end()) { + LOG("NatIsotop", pWARN) + << "Table has no elements for natural isotope Z = " << Z; + return 0; + } + + vector vec = miter->second; + for (int i; iPdgCode()==pdgcode) return vec[i]; + } + + LOG("NatIsotop", pWARN) + << "Natural isotope Z = " << Z << " has " << vec.size() << " elements" + << " (pdgcode = " << pdgcode << " was requested)"; + return 0; + +} +//____________________________________________________________________________ bool NaturalIsotopes::LoadTable(void) { // get the natural isotopes table filename @@ -157,7 +180,7 @@ bool NaturalIsotopes::LoadTable(void) LOG("NatIsotop", pDEBUG) << " - Element: " << n << ", pdg = " << pdgcode << ", A = " << atomicmass << ", abundance = " << abundance; - data = new NaturalIsotopeElementData(pdgcode, abundance); + data = new NaturalIsotopeElementData(pdgcode, abundance,atomicmass); vec.push_back(data); } fNaturalIsotopesTable.insert( diff --git a/src/Framework/ParticleData/NaturalIsotopes.h b/src/Framework/ParticleData/NaturalIsotopes.h index 4de483dfe..0fac6c948 100644 --- a/src/Framework/ParticleData/NaturalIsotopes.h +++ b/src/Framework/ParticleData/NaturalIsotopes.h @@ -37,6 +37,7 @@ class NaturalIsotopes int NElements(int Z) const; const NaturalIsotopeElementData * ElementData (int Z, int ielement) const; + const NaturalIsotopeElementData * ElementDataPdg (int Z, int pdgcode) const; private: NaturalIsotopes(); @@ -63,16 +64,18 @@ class NaturalIsotopes class NaturalIsotopeElementData { public: - NaturalIsotopeElementData() : fPdgCode(0), fAbundance(0) { } - NaturalIsotopeElementData(int code, double abundance) : fPdgCode(code), fAbundance(abundance) { } + NaturalIsotopeElementData() : fPdgCode(0), fAbundance(0), fAtomicMass(0) { } + NaturalIsotopeElementData(int code, double abundance, double atomicmass) : fPdgCode(code), fAbundance(abundance), fAtomicMass(atomicmass) { } ~NaturalIsotopeElementData() { } - int PdgCode (void) const { return fPdgCode; } - double Abundance (void) const { return fAbundance; } + int PdgCode (void) const { return fPdgCode; } + double Abundance (void) const { return fAbundance; } + double AtomicMass (void) const { return fAtomicMass; } private: int fPdgCode; double fAbundance; + double fAtomicMass; }; } // genie namespace diff --git a/src/Framework/Utils/KineUtils.cxx b/src/Framework/Utils/KineUtils.cxx index 2815354b5..9a1f8a45e 100644 --- a/src/Framework/Utils/KineUtils.cxx +++ b/src/Framework/Utils/KineUtils.cxx @@ -209,6 +209,19 @@ double genie::utils::kinematics::Jacobian( J = 1. / (kine.Q2() * kine.y()); } + // + // transformation: {x,y}|E -> {x,Q2}|E + // + else + if ( TransformMatched(fromps,tops,kPSxQ2fE,kPSxyfE,forward) ) + { + const InitialState & init_state = i->InitState(); + double Ev = init_state.ProbeE(kRfHitNucRest); + double M = init_state.Tgt().HitNucP4Ptr()->M(); + double x = kine.x(); + J = 2*x*Ev*M; + } + // // transformation: {Q2,y}|E -> {x,y}|E // diff --git a/src/Framework/Utils/XSecSplineList.cxx b/src/Framework/Utils/XSecSplineList.cxx index c7c8b30a9..d18082ab1 100644 --- a/src/Framework/Utils/XSecSplineList.cxx +++ b/src/Framework/Utils/XSecSplineList.cxx @@ -195,6 +195,13 @@ void XSecSplineList::CreateSpline(const XSecAlgorithmI * alg, SLOG("XSecSplLst", pNOTICE) << "Energy threshold for current interaction = " << Ethr << " GeV"; + if (Ethr>e_max) { + SLOG("XSecSplLst", pFATAL) << "Energy threshold higher than maximum."; + SLOG("XSecSplLst", pFATAL) << "Energy threshold = " << Ethr << " GeV"; + SLOG("XSecSplLst", pFATAL) << "Energy maximum = " << e_max << " GeV"; + return; + } + int nkb = (Ethr>e_min) ? 5 : 0; // number of knots < threshold int nka = nknots-nkb; // number of knots >= threshold diff --git a/src/Physics/Common/InitialStateAppender.cxx b/src/Physics/Common/InitialStateAppender.cxx index ee9da9dd5..6afbc8b32 100644 --- a/src/Physics/Common/InitialStateAppender.cxx +++ b/src/Physics/Common/InitialStateAppender.cxx @@ -87,7 +87,7 @@ void InitialStateAppender::AddNucleus(GHepRecord * evrec) const const ProcessInfo & proc_info = interaction->ProcInfo(); bool is_nucleus = init_state.Tgt().IsNucleus(); - if(!is_nucleus && !proc_info.IsGlashowResonance()) { + if(!is_nucleus && !proc_info.IsGlashowResonance() && !proc_info.IsPhotonCoherent()) { LOG("ISApp", pINFO) << "Not an interaction with a nuclear target - no nucleus to add"; return; diff --git a/src/Physics/Common/VertexGenerator.cxx b/src/Physics/Common/VertexGenerator.cxx index c2a179db1..0f102d6fc 100644 --- a/src/Physics/Common/VertexGenerator.cxx +++ b/src/Physics/Common/VertexGenerator.cxx @@ -117,7 +117,10 @@ TVector3 VertexGenerator::GenerateVertex(const Interaction * interaction, bool is_coh = proc_info.IsCoherentProduction() || proc_info.IsCoherentElastic(); bool is_ve = proc_info.IsInverseMuDecay() || proc_info.IsIMDAnnihilation() || - proc_info.IsNuElectronElastic(); + proc_info.IsNuElectronElastic() || + proc_info.IsGlashowResonance() || + proc_info.IsPhotonResonance() || + proc_info.IsPhotonCoherent(); if(is_coh||is_ve) { // ** For COH or ve- set a vertex positon on the nuclear boundary diff --git a/src/Physics/GlashowResonance/EventGen/GLRESGenerator.cxx b/src/Physics/GlashowResonance/EventGen/GLRESGenerator.cxx deleted file mode 100644 index 1627e9eab..000000000 --- a/src/Physics/GlashowResonance/EventGen/GLRESGenerator.cxx +++ /dev/null @@ -1,324 +0,0 @@ -//____________________________________________________________________________ -/* - Copyright (c) 2003-2020, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org - - Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory -*/ -//____________________________________________________________________________ - -#include - -#include -#include -#include - -#include "Framework/Conventions/Constants.h" -#include "Framework/GHEP/GHepStatus.h" -#include "Framework/GHEP/GHepParticle.h" -#include "Framework/GHEP/GHepRecord.h" -#include "Framework/GHEP/GHepFlags.h" -#include "Framework/EventGen/EVGThreadException.h" -#include "Framework/Messenger/Messenger.h" -#include "Framework/Numerical/RandomGen.h" -#include "Framework/Numerical/MathUtils.h" -#include "Framework/ParticleData/PDGCodes.h" -#include "Framework/ParticleData/PDGUtils.h" -#include "Framework/ParticleData/PDGLibrary.h" -#include "Physics/GlashowResonance/EventGen/GLRESGenerator.h" - -#ifdef __GENIE_PYTHIA6_ENABLED__ -#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) -#include -#else -#include -#endif -#endif // __GENIE_PYTHIA6_ENABLED__ - -using namespace genie; -using namespace genie::constants; -using namespace genie::utils::math; - -//___________________________________________________________________________ -GLRESGenerator::GLRESGenerator() : -EventRecordVisitorI("genie::GLRESGenerator") -{ - -} -//___________________________________________________________________________ -GLRESGenerator::GLRESGenerator(string config) : -EventRecordVisitorI("genie::GLRESGenerator", config) -{ - -} -//___________________________________________________________________________ -GLRESGenerator::~GLRESGenerator() -{ - -} -//___________________________________________________________________________ -void GLRESGenerator::ProcessEventRecord(GHepRecord * - #ifdef __GENIE_PYTHIA6_ENABLED__ - event // avoid unused variable warning if PYTHIA6 is not enabled - #endif -) const -{ - -#ifdef __GENIE_PYTHIA6_ENABLED__ - - Interaction * interaction = event->Summary(); - const InitialState & init_state = interaction->InitState(); - - //incoming v & struck particle & remnant nucleus - GHepParticle * nu = event -> Probe(); - GHepParticle * el = event -> HitElectron(); - GHepParticle * target = event -> TargetNucleus(); - assert(nu); - assert(el); - assert(target); - - if(target) event->AddParticle(target->Pdg(), kIStStableFinalState, 1,-1,-1,-1, *(target->P4()), *(target->X4()) ); - - LongLorentzVector p4_nu( * event->Probe()->P4() ); - LongLorentzVector p4_el( * event->HitElectron()->P4() ); - LongLorentzVector p4_Wlong( p4_nu.Px()+p4_el.Px(), p4_nu.Py()+p4_el.Py(), p4_nu.Pz()+p4_el.Pz(), p4_nu.E()+p4_el.E() ); - - long double Wmass = p4_Wlong.M(); - LOG("GLRESGenerator", pINFO) << "Wmass : " << Wmass; - - if(Wmass < fWmin) { - LOG("GLRESGenerator", pWARN) << "Low invariant mass, W = " << Wmass << " GeV!!"; - event->EventFlags()->SetBitNumber(kHadroSysGenErr, true); - genie::exceptions::EVGThreadException exception; - exception.SetReason("Could not simulate the hadronic system"); - exception.SwitchOnFastForward(); - throw exception; - return; - } - - // Final state primary lepton PDG code - int pdgl = interaction->FSPrimLeptonPdg(); - assert(pdgl!=0); - - LOG("GLRESGenerator", pINFO) << "Channel : " << interaction->FSPrimLeptonPdg(); - - if ( pdg::IsElectron(pdgl) || pdg::IsMuon(pdgl) || pdg::IsTau(pdgl) ) { - - // Get selected kinematics - long double y = interaction->Kine().y(true); - assert(y>0 && y<1); - - // Compute the neutrino and muon energy - long double Ev = init_state.ProbeE(kRfLab); - long double El = y*Ev; - - LOG("GLRESGenerator", pINFO) << "Ev = " << Ev << ", y = " << y << ", -> El = " << El; - - // Compute the momentum transfer and scattering angle - long double El2 = powl(El,2); - long double ml = interaction->FSPrimLepton()->Mass(); - long double ml2 = powl(ml,2); - long double pl = sqrtl(El2-ml2); - - long double Q2 = 2*(Ev-El)*kElectronMass + kElectronMass*kElectronMass; - long double costh = (El-0.5*(Q2+ml2)/Ev)/pl; - long double sinth = sqrtl( 1-powl(costh,2.) ); - // Randomize transverse components - RandomGen * rnd = RandomGen::Instance(); - long double phi = 2* constants::kPi * rnd->RndLep().Rndm(); - - LOG("GLRESGenerator", pINFO) << "Q2 = " << Q2 << ", cos(theta) = " << costh << ", phi = " << phi; - - long double plx = pl*sinth*cosl(phi); - long double ply = pl*sinth*sinl(phi); - long double plz = pl*costh; - - // Lepton 4-momentum in the LAB frame - LongLorentzVector p4lplong( plx, ply, plz, El ); - p4lplong.Rotate(p4_nu); - LongLorentzVector p4nulong( p4_Wlong.Px()-p4lplong.Px(), p4_Wlong.Py()-p4lplong.Py(), p4_Wlong.Pz()-p4lplong.Pz(), p4_Wlong.E()-p4lplong.E() ); - - TLorentzVector p4_W ( (double)p4_Wlong.Px(), (double)p4_Wlong.Py(), (double)p4_Wlong.Pz(), (double)p4_Wlong.E() ); - TLorentzVector p4_lpout( (double)p4lplong.Px(), (double)p4lplong.Py(), (double)p4lplong.Pz(), (double)p4lplong.E() ); - TLorentzVector p4_nuout( (double)p4nulong.Px(), (double)p4nulong.Py(), (double)p4nulong.Pz(), (double)p4nulong.E() ); - - int pdgvout = 0; - if ( pdg::IsElectron(pdgl) ) pdgvout = kPdgAntiNuE; - else if ( pdg::IsMuon(pdgl) ) pdgvout = kPdgAntiNuMu; - else if ( pdg::IsTau(pdgl) ) pdgvout = kPdgAntiNuTau; - - // Create a GHepParticle and add it to the event record - event->AddParticle( kPdgWM, kIStDecayedState, 0, -1, 5, 6, p4_W, *(nu->X4()) ); //W [mothers: nuebar_in,e_in][daugthers: nulbar_out,lp_out] - event->AddParticle( pdgvout, kIStStableFinalState, 4, -1, -1, -1, p4_nuout, *(nu->X4()) ); - event->AddParticle( pdgl, kIStStableFinalState, 4, -1, -1, -1, p4_lpout, *(nu->X4()) ); - event->Summary()->KinePtr()->SetFSLeptonP4(p4_lpout); - - } - else { - - char p6frame[10], p6nu[10], p6tgt[10]; - strcpy(p6frame, "CMS" ); - strcpy(p6nu, "nu_ebar"); - strcpy(p6tgt, "e-" ); - - - int mstp61 = fPythia->GetMSTP(61); - int mstp71 = fPythia->GetMSTP(71); - int mdme206 = fPythia->GetMDME(206,1); - int mdme207 = fPythia->GetMDME(207,1); - int mdme208 = fPythia->GetMDME(208,1); - double pmas1_W = fPythia->GetPMAS(24,1); - double pmas2_W = fPythia->GetPMAS(24,2); - double pmas2_t = fPythia->GetPMAS(6,2); - fPythia->SetMSTP(61,0); // (Default=2) master switch for initial-state QCD and QED radiation. - fPythia->SetMSTP(71,0); // (Default=2) master switch for initial-state QCD and QED radiation. - fPythia->SetMDME(206,1,0); //swicht off W decay leptonic modes - fPythia->SetMDME(207,1,0); - fPythia->SetMDME(208,1,0); - fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385) - fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation - fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation - - fPythia->Pyinit(p6frame, p6nu, p6tgt, (double)Wmass); - - fPythia->Pyevnt(); - - fPythia->SetMSTP(61,mstp61); - fPythia->SetMSTP(71,mstp71); - fPythia->SetMDME(206,1,mdme206); - fPythia->SetMDME(207,1,mdme207); - fPythia->SetMDME(208,1,mdme208); - fPythia->SetPMAS(24,1,pmas1_W); - fPythia->SetPMAS(24,2,pmas2_W); - fPythia->SetPMAS(6,2,pmas2_t); - - - // Use for debugging purposes - //fPythia->Pylist(3); - - // get LUJETS record - fPythia->GetPrimaries(); - TClonesArray * pythia_particles = (TClonesArray *) fPythia->ImportParticles("All"); - // copy PYTHIA container to a new TClonesArray so as to transfer ownership - // of the container and of its elements to the calling method - int np = pythia_particles->GetEntries(); - assert(np>0); - TClonesArray * particle_list = new TClonesArray("genie::GHepParticle", np); - particle_list->SetOwner(true); - - // Boost velocity LAB' -> Resonance rest frame - long double beta = p4_Wlong.P()/p4_Wlong.E(); - - TMCParticle * p = 0; - TIter particle_iter(pythia_particles); - while( (p = (TMCParticle *) particle_iter.Next()) ) { - - int pdgc = p->GetKF(); - int ks = p->GetKS(); - int parent = p->GetParent(); - - if ( ks==21 ) { continue; } //we dont want to save first particles from pythia (init states) - - LongLorentzVector p4long( p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy() ); - p4long.Boost(beta); - p4long.Rotate(p4_Wlong); - - TLorentzVector p4( (double)p4long.Px(), (double)p4long.Py(), (double)p4long.Pz(), (double)p4long.E() ); - - double massPDG = PDGLibrary::Instance()->Find(pdgc)->Mass(); - if ( (ks==1 || ks==4) && p4.E() < massPDG ) { - LOG("GLRESGenerator", pWARN) << "Putting at rest one stable particle generated by PYTHIA because E < m"; - LOG("GLRESGenerator", pWARN) << "PDG = " << pdgc << " // State = " << ks; - LOG("GLRESGenerator", pWARN) << "E = " << p4.E() << " // |p| = " << TMath::Sqrt(p4.P()); - LOG("GLRESGenerator", pWARN) << "p = [ " << p4.Px() << " , " << p4.Py() << " , " << p4.Pz() << " ]"; - LOG("GLRESGenerator", pWARN) << "m = " << p4.M() << " // mpdg = " << massPDG; - p4.SetXYZT(0,0,0,massPDG); - } - - // copy final state particles to the event record - GHepStatus_t ist = (ks==1 || ks==4) ? kIStStableFinalState : kIStDISPreFragmHadronicState; - - // fix numbering scheme used for mother/daughter assignments - int firstmother = -1; - int lastmother = -1; - int firstchild = -1; - int lastchild = -1; - - if ( parent < 10 ) { // I=10 is the position were W boson is stored - if ( TMath::Abs(pdgc)<7 ) { //outgoing quarks: mother will be the boson (saved in position 4) - firstmother = 4; - firstchild = p->GetFirstChild() - 6; //boson is stored in I=10 with pythia and I=4 for GENIE => shift 6 - lastchild = p->GetLastChild() - 6; - } - else if ( TMath::Abs(pdgc)==24 ) { //produced W boson: mother will be the incoming neutrino - firstmother = 0; - firstchild = p->GetFirstChild() - 6; - lastchild = p->GetLastChild() - 6; - } - else if ( pdgc==22 ) { //radiative photons: mother will be the incoming electron - firstmother = 2; - } - } - else { //rest - firstmother = parent - 6; //shift to match boson position - firstchild = (p->GetFirstChild()==0) ? p->GetFirstChild() - 1 : p->GetFirstChild() - 6; - lastchild = (p->GetLastChild()==0) ? p->GetLastChild() - 1 : p->GetLastChild() - 6; - } - - double vx = nu->X4()->X() + p->GetVx()*1e12; //pythia gives position in [mm] while genie uses [fm] - double vy = nu->X4()->Y() + p->GetVy()*1e12; - double vz = nu->X4()->Z() + p->GetVz()*1e12; - double vt = nu->X4()->T() + p->GetTime()*(units::millimeter/units::second); - TLorentzVector pos( vx, vy, vz, vt ); - - event->AddParticle(pdgc, ist, firstmother, lastmother, firstchild, lastchild, p4, pos ); - - } - - } - -#else - LOG("GLRESGenerator", pFATAL) - << "Calling GENIE/PYTHIA6 without enabling PYTHIA6"; - gAbortingInErr = true; - std::exit(1); -#endif - -} -//___________________________________________________________________________ -void GLRESGenerator::Configure(const Registry & config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void GLRESGenerator::Configure(string config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void GLRESGenerator::LoadConfig(void) -{ -#ifdef __GENIE_PYTHIA6_ENABLED__ - fPythia = TPythia6::Instance(); - - // sync GENIE/PYTHIA6 seed number - RandomGen::Instance(); - - // PYTHIA parameters only valid for HEDIS - GetParam( "Xsec-Wmin", fWmin ) ; - fPythia->SetPARP(2, fWmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes) - - int warnings; GetParam( "Warnings", warnings ) ; - int errors; GetParam( "Errors", errors ) ; - int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ; - fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed - fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed - fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term. - -#endif // __GENIE_PYTHIA6_ENABLED__ - -} -//____________________________________________________________________________ diff --git a/src/Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.cxx b/src/Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.cxx deleted file mode 100644 index 20525893e..000000000 --- a/src/Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.cxx +++ /dev/null @@ -1,121 +0,0 @@ -//____________________________________________________________________________ -/* - Copyright (c) 2003-2020, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org - - Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory -*/ -//____________________________________________________________________________ - -#include "Framework/EventGen/InteractionList.h" -#include "Framework/Interaction/Interaction.h" -#include "Framework/Messenger/Messenger.h" -#include "Framework/ParticleData/PDGCodes.h" -#include "Framework/ParticleData/PDGUtils.h" -#include "Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.h" - -using namespace genie; - -//___________________________________________________________________________ -GLRESInteractionListGenerator::GLRESInteractionListGenerator() : -InteractionListGeneratorI("genie::GLRESInteractionListGenerator") -{ - -} -//___________________________________________________________________________ -GLRESInteractionListGenerator::GLRESInteractionListGenerator(string config) : -InteractionListGeneratorI("genie::GLRESInteractionListGenerator", config) -{ - -} -//___________________________________________________________________________ -GLRESInteractionListGenerator::~GLRESInteractionListGenerator() -{ - -} -//___________________________________________________________________________ -InteractionList * - GLRESInteractionListGenerator::CreateInteractionList( - const InitialState & init_state) const -{ -// channels: -// nuebar + e- -> W- -> nuebar + e- -// nuebar + e- -> W- -> nuebar + mu- -// nuebar + e- -> W- -> nuebar + tau- -// nuebar + e- -> W- -> hadrons - - if(init_state.ProbePdg() != kPdgAntiNuE) { - LOG("IntLst", pDEBUG) - << "Return *null* interaction list"; - return 0; - } - - InitialState init(init_state); - init_state.TgtPtr()->SetHitNucPdg(0); - - ProcessInfo proc_info(kScGlashowResonance, kIntWeakCC); - - InteractionList * intlist = new InteractionList; - - if (fIsMu) { - Interaction * interaction = new Interaction(init_state, proc_info); - XclsTag exclusive_tag; - exclusive_tag.SetFinalLepton(kPdgMuon); - interaction->SetExclTag(exclusive_tag); - intlist->push_back(interaction); - } - else if (fIsTau) { - Interaction * interaction = new Interaction(init_state, proc_info); - XclsTag exclusive_tag; - exclusive_tag.SetFinalLepton(kPdgTau); - interaction->SetExclTag(exclusive_tag); - intlist->push_back(interaction); - } - else if (fIsEle) { - Interaction * interaction = new Interaction(init_state, proc_info); - XclsTag exclusive_tag; - exclusive_tag.SetFinalLepton(kPdgElectron); - interaction->SetExclTag(exclusive_tag); - intlist->push_back(interaction); - } - else if (fIsHad) { - Interaction * interaction = new Interaction(init_state, proc_info); - XclsTag exclusive_tag; - exclusive_tag.SetFinalLepton(kPdgPiP); - interaction->SetExclTag(exclusive_tag); - intlist->push_back(interaction); - } - - if(intlist->size() == 0) { - LOG("IntLst", pERROR) - << "Returning NULL InteractionList for init-state: " - << init_state.AsString(); - delete intlist; - return 0; - } - - return intlist; -} -//___________________________________________________________________________ -void GLRESInteractionListGenerator::Configure(const Registry & config) -{ - Algorithm::Configure(config); - this->LoadConfigData(); -} -//____________________________________________________________________________ -void GLRESInteractionListGenerator::Configure(string config) -{ - Algorithm::Configure(config); - this->LoadConfigData(); -} -//____________________________________________________________________________ -void GLRESInteractionListGenerator::LoadConfigData(void) -{ - - GetParamDef("is-Mu", fIsMu, false ) ; - GetParamDef("is-Tau", fIsTau, false ) ; - GetParamDef("is-Ele", fIsEle, false ) ; - GetParamDef("is-Had", fIsHad, false ) ; - -} \ No newline at end of file diff --git a/src/Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.h b/src/Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.h deleted file mode 100644 index 2b8945875..000000000 --- a/src/Physics/GlashowResonance/EventGen/GLRESInteractionListGenerator.h +++ /dev/null @@ -1,54 +0,0 @@ -//____________________________________________________________________________ -/*! - -\class genie::GLRESInteractionListGenerator - -\brief Generates a list of all the interactions that cab be generated by - the GLRES EventGenerator - -\author Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory - -\created December 14, 2009 - -\cpright Copyright (c) 2003-2020, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org -*/ -//____________________________________________________________________________ - -#ifndef _GLASHOW_RESONANCE_INTERACTION_GENERATOR_H_ -#define _GLASHOW_RESONANCE_INTERACTION_GENERATOR_H_ - -#include "Framework/EventGen/InteractionListGeneratorI.h" - -namespace genie { - -class GLRESInteractionListGenerator : public InteractionListGeneratorI { - -public : - GLRESInteractionListGenerator(); - GLRESInteractionListGenerator(string config); - ~GLRESInteractionListGenerator(); - - // implement the InteractionListGeneratorI interface - InteractionList * CreateInteractionList(const InitialState & init) const; - - // overload the Algorithm::Configure() methods to load private data - // members from configuration options - void Configure(const Registry & config); - void Configure(string config); - -private: - - void LoadConfigData(void); - - bool fIsMu; - bool fIsTau; - bool fIsEle; - bool fIsHad; - -}; - -} // genie namespace - -#endif // _GLASHOW_RESONANCE_INTERACTION_GENERATOR_H_ diff --git a/src/Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.cxx b/src/Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.cxx deleted file mode 100644 index 9f28c3eb8..000000000 --- a/src/Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.cxx +++ /dev/null @@ -1,211 +0,0 @@ -//____________________________________________________________________________ -/* - Copyright (c) 2003-2019, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org - or see $GENIE/LICENSE - - Author: Alfonso Garcia - NIKHEF (Amsterdam) - - For the class documentation see the corresponding header file. - - Important revisions after version 2.0.0 : - -*/ -//____________________________________________________________________________ - -#include "Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.h" -#include "Framework/Conventions/Controls.h" -#include "Framework/Conventions/KinePhaseSpace.h" -#include "Framework/EventGen/EVGThreadException.h" -#include "Framework/EventGen/EventGeneratorI.h" -#include "Framework/EventGen/RunningThreadInfo.h" -#include "Framework/GHEP/GHepRecord.h" -#include "Framework/GHEP/GHepFlags.h" -#include "Framework/Messenger/Messenger.h" -#include "Framework/Numerical/RandomGen.h" -#include "Framework/Utils/KineUtils.h" -#include "Framework/Utils/Range1.h" - -using namespace genie; -using namespace genie::controls; -using namespace genie::utils; - -//___________________________________________________________________________ -GLRESKinematicsGenerator::GLRESKinematicsGenerator() : -KineGeneratorWithCache("genie::GLRESKinematicsGenerator") -{ - -} -//___________________________________________________________________________ -GLRESKinematicsGenerator::GLRESKinematicsGenerator(string config) : -KineGeneratorWithCache("genie::GLRESKinematicsGenerator", config) -{ - -} -//___________________________________________________________________________ -GLRESKinematicsGenerator::~GLRESKinematicsGenerator() -{ - -} -//___________________________________________________________________________ -void GLRESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const -{ - if(fGenerateUniformly) { - LOG("GLRESKinematics", pNOTICE) - << "Generating kinematics uniformly over the allowed phase space"; - } - - //-- Get the random number generators - RandomGen * rnd = RandomGen::Instance(); - - //-- Access cross section algorithm for running thread - RunningThreadInfo * rtinfo = RunningThreadInfo::Instance(); - const EventGeneratorI * evg = rtinfo->RunningThread(); - fXSecModel = evg->CrossSectionAlg(); - - //-- For the subsequent kinematic selection with the rejection method: - // Calculate the max differential cross section or retrieve it from the - // cache. Throw an exception and quit the evg thread if a non-positive - // value is found. - // If the kinematics are generated uniformly over the allowed phase - // space the max xsec is irrelevant - double xsec_max = this->MaxXSec(evrec); - - //-- y range - const KPhaseSpace & kps = evrec->Summary()->PhaseSpace(); - Range1D_t yl = kps.Limits(kKVy); - double ymin = yl.min; - double ymax = yl.max; - double dy = ymax-ymin; - - double xsec = -1; - Interaction * interaction = evrec->Summary(); - - //-- Try to select a valid inelastisity y - unsigned int iter = 0; - bool accept = false; - while(1) { - iter++; - if(iter > kRjMaxIterations) { - LOG("GLRESKinematics", pWARN) - << "*** Could not select a valid y after " - << iter << " iterations"; - evrec->EventFlags()->SetBitNumber(kKineGenErr, true); - genie::exceptions::EVGThreadException exception; - exception.SetReason("Couldn't select kinematics"); - exception.SwitchOnFastForward(); - throw exception; - } - - double y = ymin + dy * rnd->RndKine().Rndm(); - interaction->KinePtr()->Sety(y); - - LOG("GLRESKinematics", pINFO) << "Trying: y = " << y; - - //-- computing cross section for the current kinematics - xsec = fXSecModel->XSec(interaction, kPSyfE); - - this->AssertXSecLimits(interaction, xsec, xsec_max); - - double t = xsec_max * rnd->RndKine().Rndm(); - LOG("GLRESKinematics", pDEBUG) << "xsec= "<< xsec<< ", J= 1, Rnd= "<< t; - - accept = (tAsString(); - SLOG("GLRESKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec; - SLOG("GLRESKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel; - - return max_xsec; -} -//___________________________________________________________________________ -double GLRESKinematicsGenerator::Energy(const Interaction * interaction) const -{ -// Override the base class Energy() method to cache the max xsec for the -// neutrino energy in the LAB rather than in the hit nucleon rest frame. - - const InitialState & init_state = interaction->InitState(); - double E = init_state.ProbeE(kRfLab); - return E; -} -//___________________________________________________________________________ -void GLRESKinematicsGenerator::Configure(const Registry & config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void GLRESKinematicsGenerator::Configure(string config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void GLRESKinematicsGenerator::LoadConfig(void) -{ -// Reads its configuration data from its configuration Registry and loads them -// in private data members to avoid looking up at the Registry all the time. - - //-- Safety factor for the maximum differential cross section - GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 2. ) ; - - //-- Maximum allowed fractional cross section deviation from maxim cross - // section used in rejection method - GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ; - assert(fMaxXSecDiffTolerance>=0); - -} -//____________________________________________________________________________ diff --git a/src/Physics/GlashowResonance/EventGen/LinkDef.h b/src/Physics/GlashowResonance/EventGen/LinkDef.h deleted file mode 100644 index aed3a665e..000000000 --- a/src/Physics/GlashowResonance/EventGen/LinkDef.h +++ /dev/null @@ -1,13 +0,0 @@ -#ifdef __CINT__ - -#pragma link off all globals; -#pragma link off all classes; -#pragma link off all functions; - -#pragma link C++ namespace genie; - -#pragma link C++ class genie::GLRESKinematicsGenerator; -#pragma link C++ class genie::GLRESGenerator; -#pragma link C++ class genie::GLRESInteractionListGenerator; - -#endif diff --git a/src/Physics/GlashowResonance/XSection/GLRESPXSec.cxx b/src/Physics/GlashowResonance/XSection/GLRESPXSec.cxx deleted file mode 100644 index 9ecd07729..000000000 --- a/src/Physics/GlashowResonance/XSection/GLRESPXSec.cxx +++ /dev/null @@ -1,181 +0,0 @@ -//____________________________________________________________________________ -/* - Copyright (c) 2003-2020, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org - - Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory -*/ -//____________________________________________________________________________ - -#include - -#include "Physics/XSectionIntegration/XSecIntegratorI.h" -#include "Framework/Conventions/Constants.h" -#include "Framework/Conventions/RefFrame.h" -#include "Framework/Messenger/Messenger.h" -#include "Framework/ParticleData/PDGCodes.h" -#include "Framework/ParticleData/PDGUtils.h" -#include "Framework/ParticleData/PDGLibrary.h" -#include "Framework/Utils/KineUtils.h" -#include "Physics/GlashowResonance/XSection/GLRESPXSec.h" - -using namespace genie; -using namespace genie::constants; - -//____________________________________________________________________________ -GLRESPXSec::GLRESPXSec() : -XSecAlgorithmI("genie::GLRESPXSec") -{ - -} -//____________________________________________________________________________ -GLRESPXSec::GLRESPXSec(string config) : -XSecAlgorithmI("genie::GLRESPXSec", config) -{ - -} -//____________________________________________________________________________ -GLRESPXSec::~GLRESPXSec() -{ - -} -//____________________________________________________________________________ -double GLRESPXSec::XSec( - const Interaction * interaction, KinePhaseSpace_t kps) const -{ - - if(! this -> ValidProcess (interaction) ) return 0.; - if(! this -> ValidKinematics (interaction) ) return 0.; - - const InitialState & init_state = interaction -> InitState(); - - double E = init_state.ProbeE(kRfLab); - double s = 2 * kElectronMass * E + kElectronMass2; - - // The W limit is because hadronization might have issues at low W (as in PYTHIA6). - // To be consistent the cross section must be computed in the W range where the events are generated. - if ( TMath::Sqrt(s) Kine(); - const XclsTag & xclstag = interaction -> ExclTag(); - int final_lepton = xclstag.FinalLeptonPdg(); - - double y = kinematics.y(); - double s0 = 2 * kElectronMass * E; - - double Gw = PDGLibrary::Instance()->Find(kPdgWM)->Width(); //PDGLibrary::Instance()->Find(kPdgWM)->Width() //genhen=2.124 - double gf = kGF2/kPi; //kGF2/kPi //genhen=1.16637e-05*1.16637e-05/3.141592653589793 - double m_w = kMw; //kMw //genhen=80.425 - double m_z = kMz; //kMz //genhen=91.1876 - double Sin2thw = 1 - kMw2 / kMz2; //1 - kMw2 / kMz2 //genhen=0.2312 - double branch = 64.41/10.63; //branching ratio hadrons/muons //genhen=67.96/10.57 - - double Gw2 = TMath::Power(Gw, 2); - double m_w2 = TMath::Power(m_w,2); - double m_z2 = TMath::Power(m_z,2); - double Sin4thw = TMath::Power(Sin2thw,2); - - double prop = TMath::Power(1-s/m_w2, 2) + Gw2/m_w2; - - double xsec = 0; - if ( pdg::IsMuon(final_lepton) ) { - double ml2 = kMuonMass2; - xsec = s0*TMath::Power(1-y, 2) + (3*kElectronMass2+ml2)*(1-y) + kElectronMass*(kElectronMass2+ml2)/E; - xsec *= gf/prop; - } - else if ( pdg::IsTau(final_lepton) ) { - double ml2 = kTauMass2; - xsec = s0*TMath::Power(1-y, 2) + (3*kElectronMass2+ml2)*(1-y) + kElectronMass*(kElectronMass2+ml2)/E; - xsec *= gf/prop; - } - else if ( pdg::IsPion(final_lepton) ) { - double ml2 = kMuonMass2; - xsec = ( s0*TMath::Power(1-y, 2) + (3*kElectronMass2+ml2)*(1-y) + kElectronMass*(kElectronMass2+ml2)/E ) * branch; - xsec *= gf/prop; - } - else if ( pdg::IsElectron(final_lepton) ) { - double u = 2 * kElectronMass2 - 2 * kElectronMass * E * y; - double t = 2 * kElectronMass2 - s - u; - double L = Sin2thw - 1./2.; - - double x1 = L*m_z2*(s-m_w2) + m_w2*(u-m_z2); - double y1 = L * m_z2 * Gw * m_w; - double x2 = (u - m_z2) * (s - m_w2); - double y2 = (u - m_z2) * Gw * m_w; - double y3 = ( x1*x2 + y1*y2 ) / ( TMath::Power(x2,2) + TMath::Power(y2,2) ); - double x3 = ( x2*y1 - x1*y2 ) / ( TMath::Power(x2,2) + TMath::Power(y2,2) ); - - xsec = gf * s0 * ( Sin4thw/TMath::Power(u/m_z2-1,2) + (TMath::Power(x3,2)+TMath::Power(y3,2))*TMath::Power(t-kElectronMass2,2)/TMath::Power(s0,2) ); - } - - - LOG("GLRESPXSec", pINFO) << "dxsec/dy (E= " << E << ", y= " << y << ") = " << xsec; - - //----- The algorithm computes dxsec/dy - // Check whether variable tranformation is needed - if(kps!=kPSyfE) { - double J = utils::kinematics::Jacobian(interaction,kPSyfE,kps); - xsec *= J; - } - - //----- If requested return the free electron xsec even for nuclear target - if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec; - - //----- Scale for the number of scattering centers at the target - int Ne = init_state.Tgt().Z(); // num of scattering centers - xsec *= Ne; - - return xsec; - -} -//____________________________________________________________________________ -double GLRESPXSec::Integral(const Interaction * interaction) const -{ - double xsec = fXSecIntegrator->Integrate(this,interaction); - - return xsec; -} -//____________________________________________________________________________ -bool GLRESPXSec::ValidProcess(const Interaction* interaction) const -{ - if(interaction->TestBit(kISkipProcessChk)) return true; - - const InitialState & init_state = interaction->InitState(); - const ProcessInfo & proc_info = interaction->ProcInfo(); - - bool nuok = pdg::IsAntiNuE(init_state.ProbePdg()); - bool nucok = !(init_state.Tgt().HitNucIsSet()); - bool ccprcok = proc_info.IsWeakCC(); - - if ( !nuok ) return false; - if ( !nucok ) return false; - if ( !ccprcok ) return false; - - return true; -} -//____________________________________________________________________________ -void GLRESPXSec::Configure(const Registry & config) -{ - - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void GLRESPXSec::Configure(string config) -{ - Algorithm::Configure(config); - this->LoadConfig(); -} -//____________________________________________________________________________ -void GLRESPXSec::LoadConfig(void) -{ - - //-- load the differential cross section integrator - fXSecIntegrator = dynamic_cast (this->SubAlg("XSec-Integrator")); - assert(fXSecIntegrator); - - GetParam( "Xsec-Wmin", fWmin ) ; - - -} \ No newline at end of file diff --git a/src/Physics/GlashowResonance/XSection/GLRESXSec.h b/src/Physics/GlashowResonance/XSection/GLRESXSec.h deleted file mode 100644 index ad588a592..000000000 --- a/src/Physics/GlashowResonance/XSection/GLRESXSec.h +++ /dev/null @@ -1,57 +0,0 @@ -//____________________________________________________________________________ -/*! - -\class genie::GLRESXSec - -\brief nubar + e- scattering glashow resonance. Integrates the loaded - differential cross section model. An analytical cross section - model also exists, so you cal also use that if you do not apply - any kinematical cuts. - - The cross section algorithm handles: - - nuebar + e- -> nuebar + e- [CC + NC + interference] - - nuebar + e- -> numubar + mu- [CC] - - nuebar + e- -> nutaubar + tau- [CC] - - nuebar + e- -> hadrons [CC] - - Is a concrete implementation of the XSecIntegratorI interface. \n - -\author Alfonso Garcia - NIKHEF (Amsterdam) - -\created November 8, 2019 - -\cpright Copyright (c) 2003-2019, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org - or see $GENIE/LICENSE -*/ -//____________________________________________________________________________ - -#ifndef _GLASHOW_RESONANCE_XSEC_H_ -#define _GLASHOW_RESONANCE_XSEC_H_ - -#include "Physics/XSectionIntegration/XSecIntegratorI.h" - -namespace genie { - -class GLRESXSec : public XSecIntegratorI { - -public: - GLRESXSec(); - GLRESXSec(string config); - virtual ~GLRESXSec(); - - //! XSecIntegratorI interface implementation - double Integrate(const XSecAlgorithmI * model, const Interaction * i) const; - - //! Overload the Algorithm::Configure() methods to load private data - //! members from configuration options - void Configure(const Registry & config); - void Configure(string config); - -private: - void LoadConfig (void); -}; - -} // genie namespace -#endif // _GLASHOW_RESONANCE_XSEC_H_ diff --git a/src/Physics/HEDIS/EventGen/HEDISGenerator.cxx b/src/Physics/HEDIS/EventGen/HEDISGenerator.cxx index 2aa527c88..c033fdfd9 100644 --- a/src/Physics/HEDIS/EventGen/HEDISGenerator.cxx +++ b/src/Physics/HEDIS/EventGen/HEDISGenerator.cxx @@ -60,6 +60,8 @@ void HEDISGenerator::ProcessEventRecord(GHepRecord * evrec) const //-- Add the target remnant this->AddTargetNucleusRemnant(evrec); + GHepParticle * target = evrec -> TargetNucleus(); + if(target) evrec->Particle(evrec->RemnantNucleusPosition())->SetStatus(kIStFinalStateNuclearRemnant); //-- Add the primary lepton this->AddPrimaryLepton(evrec); diff --git a/src/Physics/HEDIS/EventGen/HEDISInteractionListGenerator.cxx b/src/Physics/HEDIS/EventGen/HEDISInteractionListGenerator.cxx index b09efea90..44756ff1e 100644 --- a/src/Physics/HEDIS/EventGen/HEDISInteractionListGenerator.cxx +++ b/src/Physics/HEDIS/EventGen/HEDISInteractionListGenerator.cxx @@ -227,3 +227,5 @@ void HEDISInteractionListGenerator::LoadConfigData(void) GetParamDef("is-NC", fIsNC, false ) ; } + + diff --git a/src/Physics/HEDIS/EventGen/LinkDef.h b/src/Physics/HEDIS/EventGen/LinkDef.h index a46231df0..dca6d1f6c 100644 --- a/src/Physics/HEDIS/EventGen/LinkDef.h +++ b/src/Physics/HEDIS/EventGen/LinkDef.h @@ -6,7 +6,6 @@ #pragma link C++ namespace genie; -#pragma link C++ class genie::LongLorentzVector; #pragma link C++ class genie::HEDISKinematicsGenerator; #pragma link C++ class genie::HEDISGenerator; #pragma link C++ class genie::HEDISInteractionListGenerator; diff --git a/src/Physics/HEDIS/XSection/HEDISStrucFunc.cxx b/src/Physics/HEDIS/XSection/HEDISStrucFunc.cxx index 7a1f2e0f9..f258c88ae 100644 --- a/src/Physics/HEDIS/XSection/HEDISStrucFunc.cxx +++ b/src/Physics/HEDIS/XSection/HEDISStrucFunc.cxx @@ -225,6 +225,12 @@ HEDISStrucFunc::HEDISStrucFunc(string basedir, SF_info sfinfo) APFEL::SetMassScheme("ZM-VFNS"); APFEL::SetPoleMasses(mPDFQrk[4],mPDFQrk[5],mPDFQrk[5]+0.1); } + else if (fSF.Scheme=="GGHR") { + APFEL::SetMassScheme("FFNS5"); + APFEL::SetPoleMasses(mPDFQrk[4],mPDFQrk[5],mPDFQrk[6]); + APFEL::SetMaxFlavourPDFs(5); + APFEL::SetMaxFlavourAlpha(5); + } else { LOG("HEDISStrucFunc", pERROR) << "Mass Scheme is not set properly"; assert(0); diff --git a/src/Physics/HEDIS/XSection/HEDISXSec.cxx b/src/Physics/HEDIS/XSection/HEDISXSec.cxx index c229410e3..45238f909 100644 --- a/src/Physics/HEDIS/XSection/HEDISXSec.cxx +++ b/src/Physics/HEDIS/XSection/HEDISXSec.cxx @@ -127,8 +127,9 @@ double HEDISXSec::Integrate( // If a GSL option has been chosen, then the total xsec is recomptued ROOT::Math::IBaseFunctionMultiDim * func = new utils::gsl::d2XSec_dlog10xdlog10Q2_E(model, interaction); ROOT::Math::IntegrationMultiDim::Type ig_type = utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); - double abstol = 1; //We mostly care about relative tolerance. - ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); + double abstol = 0; //In vegas we care about number of interactions + double reltol = 0; //In vegas we care about number of interactions + ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, reltol, fGSLMaxEval); double kine_min[2] = { TMath::Log10(xl.min), TMath::Log10(Q2l.min) }; double kine_max[2] = {TMath::Log10(xl.max), TMath::Log10(Q2l.max) }; xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); @@ -158,14 +159,11 @@ void HEDISXSec::LoadConfig(void) { // Get GSL integration type & relative tolerance - GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ; - GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ; + GetParamDef( "gsl-integration-type", fGSLIntgType, string("vegas") ) ; - int max_eval, min_eval ; + int max_eval ; GetParamDef( "gsl-max-eval", max_eval, 500000 ) ; - GetParamDef( "gsl-min-eval", min_eval, 10000 ) ; fGSLMaxEval = (unsigned int) max_eval ; - fGSLMinEval = (unsigned int) min_eval ; // Limits from the SF tables that are useful to reduce computation // time of the total cross section diff --git a/src/Physics/HEDIS/XSection/Makefile b/src/Physics/HEDIS/XSection/Makefile index 37d88e9b0..ccbd7270e 100644 --- a/src/Physics/HEDIS/XSection/Makefile +++ b/src/Physics/HEDIS/XSection/Makefile @@ -25,14 +25,4 @@ install : install-inc install-lib # include $(GENIE)/src/make/Make.std-package-targets -ifeq ($(strip $(GOPT_ENABLE_APFEL)),YES) - # extra flags, include paths, and libraries to link to - CPP_INCLUDES += $(APFEL_INCLUDES) - ROOT_DICT_GEN_INCLUDES += $(APFEL_INCLUDES) - EXTRA_EXT_LIBS += $(APFEL_LIBRARIES) -else - $(info $(PACKAGE) not build against APFEL++) -endif - - FORCE: diff --git a/src/Physics/HELepton/EventGen/GLRESGenerator.cxx b/src/Physics/HELepton/EventGen/GLRESGenerator.cxx new file mode 100644 index 000000000..1c42eb60a --- /dev/null +++ b/src/Physics/HELepton/EventGen/GLRESGenerator.cxx @@ -0,0 +1,297 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include + +#include +#include +#include + +#include "Framework/Conventions/Constants.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Physics/HELepton/EventGen/GLRESGenerator.h" + +#ifdef __GENIE_PYTHIA6_ENABLED__ +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) +#include +#else +#include +#endif +#endif // __GENIE_PYTHIA6_ENABLED__ + +using namespace genie; +using namespace genie::constants; +using namespace genie::utils::math; + +//___________________________________________________________________________ +GLRESGenerator::GLRESGenerator() : +EventRecordVisitorI("genie::GLRESGenerator") +{ + born = new Born(); +} +//___________________________________________________________________________ +GLRESGenerator::GLRESGenerator(string config) : +EventRecordVisitorI("genie::GLRESGenerator", config) +{ + +} +//___________________________________________________________________________ +GLRESGenerator::~GLRESGenerator() +{ + +} +//___________________________________________________________________________ +void GLRESGenerator::ProcessEventRecord(GHepRecord * + #ifdef __GENIE_PYTHIA6_ENABLED__ + event // avoid unused variable warning if PYTHIA6 is not enabled + #endif +) const +{ + +#ifdef __GENIE_PYTHIA6_ENABLED__ + + Interaction * interaction = event->Summary(); + const InitialState & init_state = interaction->InitState(); + + //incoming v & struck particle & remnant nucleus + GHepParticle * nu = event->Probe(); + GHepParticle * el = event->HitElectron(); + + GHepParticle * target = event -> TargetNucleus(); + if(target) event->AddParticle(target->Pdg(), kIStFinalStateNuclearRemnant, 1,-1,-1,-1, *(target->P4()), *(target->X4()) ); + + TVector3 unit_nu = nu->P4()->Vect().Unit(); + + long double mlout = interaction->FSPrimLepton()->Mass(); //mass of charged lepton + long double mlin = kElectronMass; //mass of incoming charged lepton + + long double Enuin = init_state.ProbeE(kRfLab); + long double s = born->GetS(mlin,Enuin); + + long double n1 = interaction->Kine().GetKV(kKVn1); + long double n2 = interaction->Kine().GetKV(kKVn2); + + long double costhCM = n1; + long double sinthCM = sqrtl(1-costhCM*costhCM); + + long double t = born->GetT(mlin,mlout,s,n1); + long double zeta = born->GetReAlpha()/kPi*(2.0*logl(sqrtl(-t)/kElectronMass)-1.0); + long double omx = powl(n2, 1.0/zeta ); + long double s_r = s*( 1.-omx ); + + // Boost velocity CM -> LAB + long double EnuinCM = (s_r-mlin*mlin)/sqrtl(s_r)/2.; + long double beta = (powl(Enuin,2)-powl(EnuinCM,2))/(powl(Enuin,2)+powl(EnuinCM,2)); + + // Final state primary lepton PDG code + int pdgl = interaction->FSPrimLeptonPdg(); + assert(pdgl!=0); + + if ( pdg::IsElectron(TMath::Abs(pdgl)) || pdg::IsMuon(TMath::Abs(pdgl)) || pdg::IsTau(TMath::Abs(pdgl)) ) { + + long double ElpoutCM = (s_r+mlout*mlout)/sqrtl(s_r)/2.; + long double EnuoutCM = (s_r-mlout*mlout)/sqrtl(s_r)/2.; + LongLorentzVector p4_lpout( 0., EnuoutCM*sinthCM, EnuoutCM*costhCM, ElpoutCM ); + LongLorentzVector p4_nuout( 0., -EnuoutCM*sinthCM, -EnuoutCM*costhCM, EnuoutCM ); + + p4_lpout.BoostZ(beta); + p4_nuout.BoostZ(beta); + + TLorentzVector p4lp_o( (double)p4_lpout.Px(), (double)p4_lpout.Py(), (double)p4_lpout.Pz(), (double)p4_lpout.E() ); + TLorentzVector p4nu_o( (double)p4_nuout.Px(), (double)p4_nuout.Py(), (double)p4_nuout.Pz(), (double)p4_nuout.E() ); + + // Randomize transverse components + RandomGen * rnd = RandomGen::Instance(); + double phi = 2* kPi * rnd->RndLep().Rndm(); + p4lp_o.RotateZ(phi); + p4nu_o.RotateZ(phi); + + //rotate from LAB=[0,0,Ev,Ev]->[px,py,pz,E] + p4lp_o.RotateUz(unit_nu); + p4nu_o.RotateUz(unit_nu); + + int pdgvout = 0; + if ( pdg::IsElectron(pdgl) ) pdgvout = kPdgAntiNuE; + else if ( pdg::IsPositron(pdgl) ) pdgvout = kPdgNuE; + else if ( pdg::IsMuon(pdgl) ) pdgvout = kPdgAntiNuMu; + else if ( pdg::IsAntiMuon(pdgl) ) pdgvout = kPdgNuMu; + else if ( pdg::IsTau(pdgl) ) pdgvout = kPdgAntiNuTau; + else if ( pdg::IsAntiTau(pdgl) ) pdgvout = kPdgNuTau; + + int pdgboson = pdg::IsNeutrino(init_state.ProbePdg()) ? kPdgWP : kPdgWM; + + // Create a GHepParticle and add it to the event record + event->AddParticle( pdgboson, kIStDecayedState, 0, -1, 5, 6, *nu->P4()+*el->P4(), *(nu->X4()) ); //W [mothers: nuebar_in,e_in][daugthers: nulbar_out,lp_out] + event->AddParticle( pdgl, kIStStableFinalState, 4, -1, -1, -1, p4lp_o, *(nu->X4()) ); + event->AddParticle( pdgvout, kIStStableFinalState, 4, -1, -1, -1, p4nu_o, *(nu->X4()) ); + event->Summary()->KinePtr()->SetFSLeptonP4(p4lp_o); + + } + else { + + char p6frame[10],p6nu[10],p6tgt[10]; + strcpy(p6frame, "CMS" ); + strcpy(p6nu, "nu_ebar" ); + strcpy(p6tgt, "e-" ); + + int def61 = fPythia->GetMSTP(61); + int def71 = fPythia->GetMSTP(71); + int def206 = fPythia->GetMDME(206,1); + int def207 = fPythia->GetMDME(207,1); + int def208 = fPythia->GetMDME(208,1); + fPythia->SetMSTP(61,0); // (Default=2) master switch for initial-state QCD and QED radiation. + fPythia->SetMSTP(71,0); // (Default=2) master switch for initial-state QCD and QED radiation. + fPythia->SetMDME(206,1,0); //swicht off W decay leptonic modes + fPythia->SetMDME(207,1,0); + fPythia->SetMDME(208,1,0); + + fPythia->Pyinit(p6frame, p6nu, p6tgt, sqrtl(s_r)); + fPythia->Pyevnt(); + + fPythia->SetMSTP(61,def61); + fPythia->SetMSTP(71,def71); + fPythia->SetMDME(206,1,def206); + fPythia->SetMDME(207,1,def207); + fPythia->SetMDME(208,1,def208); + + // get LUJETS record + fPythia->GetPrimaries(); + TClonesArray * pythia_particles = (TClonesArray *) fPythia->ImportParticles("All"); + int np = pythia_particles->GetEntries(); + assert(np>0); + + TMCParticle * particle = 0; + TIter piter(pythia_particles); + while( (particle = (TMCParticle *) piter.Next()) ) { + + int pdgc = particle->GetKF(); + int ks = particle->GetKS(); + + if ( ks==21 ) { continue; } //we dont want to save first particles from pythia (init states) + + LongLorentzVector p4longo(particle->GetPx(), particle->GetPy(), particle->GetPz(), particle->GetEnergy()); + p4longo.BoostZ(beta); + + TLorentzVector p4o( (double)p4longo.Px(), (double)p4longo.Py(), (double)p4longo.Pz(), (double)p4longo.E() ); + p4o.RotateUz(unit_nu); + + TParticlePDG * part = PDGLibrary::Instance()->Find(pdgc); + if ( (ks==1 || ks==4) && p4o.E() < part->Mass() ) { + LOG("GLRESGenerator", pWARN) << "Putting at rest one stable particle generated by PYTHIA because E < m"; + LOG("GLRESGenerator", pWARN) << "PDG = " << pdgc << " // State = " << ks; + LOG("GLRESGenerator", pWARN) << "E = " << p4o.E() << " // |p| = " << TMath::Sqrt(p4o.P()); + LOG("GLRESGenerator", pWARN) << "p = [ " << p4o.Px() << " , " << p4o.Py() << " , " << p4o.Pz() << " ]"; + LOG("GLRESGenerator", pWARN) << "m = " << p4o.M() << " // mpdg = " << part->Mass(); + p4o.SetXYZT(0,0,0,part->Mass()); + } + + // copy final state particles to the event record + GHepStatus_t ist = (ks==1 || ks==4) ? kIStStableFinalState : kIStDISPreFragmHadronicState; + + // fix numbering scheme used for mother/daughter assignments + int firstmother = -1; + int lastmother = -1; + int firstchild = -1; + int lastchild = -1; + + if ( particle->GetParent() < 10 ) { + if ( TMath::Abs(pdgc)<7 ) { //outgoing quarks: mother will be the boson (saved in position 4) + firstmother = 4; + firstchild = particle->GetFirstChild() - 6; + lastchild = particle->GetLastChild() - 6; + } + else if ( TMath::Abs(pdgc)==24 ) { //produced W boson: mother will be the incoming neutrino + firstmother = 0; + firstchild = particle->GetFirstChild() - 6; + lastchild = particle->GetLastChild() - 6; + } + else if ( pdgc==22 ) { //radiative photons: mother will be the incoming electron + firstmother = 2; + } + } + else { //rest + firstmother = particle->GetParent() - 6; //shift to match boson position + firstchild = (particle->GetFirstChild()==0) ? particle->GetFirstChild() - 1 : particle->GetFirstChild() - 6; + lastchild = (particle->GetLastChild()==0) ? particle->GetLastChild() - 1 : particle->GetLastChild() - 6; + } + + double lightspeed = 299792458e3; //c in mm/s. Used for time in PYTHIA t[s]=t_pythia[mm]/c[mm/s] + double vx = nu->X4()->X() + particle->GetVx()*1e12; //pythia gives position in [mm] while genie uses [fm] + double vy = nu->X4()->Y() + particle->GetVy()*1e12; + double vz = nu->X4()->Z() + particle->GetVz()*1e12; + double vt = nu->X4()->T() + particle->GetTime()/lightspeed; + TLorentzVector pos( vx, vy, vz, vt ); + + event->AddParticle(pdgc, ist, firstmother, lastmother, firstchild, lastchild, p4o, pos ); + + } + + delete particle; + pythia_particles->Clear("C"); + + } + +#else + LOG("GLRESGenerator", pFATAL) + << "Calling GENIE/PYTHIA6 without enabling PYTHIA6"; + gAbortingInErr = true; + std::exit(1); +#endif + +} +//___________________________________________________________________________ +void GLRESGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void GLRESGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void GLRESGenerator::LoadConfig(void) +{ + +#ifdef __GENIE_PYTHIA6_ENABLED__ + fPythia = TPythia6::Instance(); + + // sync GENIE/PYTHIA6 seed number + RandomGen::Instance(); + + // PYTHIA parameters only valid for HEDIS + double wmin; GetParam( "Xsec-Wmin", wmin ) ; + int warnings; GetParam( "Warnings", warnings ) ; + int errors; GetParam( "Errors", errors ) ; + int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ; + fPythia->SetPARP(2, wmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes) + fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed + fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed + fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term. + + fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385) + fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation + fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation +#endif // __GENIE_PYTHIA6_ENABLED__ + +} +//____________________________________________________________________________ diff --git a/src/Physics/GlashowResonance/EventGen/GLRESGenerator.h b/src/Physics/HELepton/EventGen/GLRESGenerator.h similarity index 80% rename from src/Physics/GlashowResonance/EventGen/GLRESGenerator.h rename to src/Physics/HELepton/EventGen/GLRESGenerator.h index b46fe5459..ee28e3901 100644 --- a/src/Physics/GlashowResonance/EventGen/GLRESGenerator.h +++ b/src/Physics/HELepton/EventGen/GLRESGenerator.h @@ -3,15 +3,16 @@ \class genie::GLRESGenerator -\brief Glashow resonance event generator +\brief Generator for glashow resonance. -\author Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory +\author Alfonso Garcia + IFIC & Harvard University -\created Feb 15, 2008 +\created Dec 8, 2021 \cpright Copyright (c) 2003-2020, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE */ //____________________________________________________________________________ @@ -21,6 +22,7 @@ #define __GENIE_PYTHIA6_ENABLED__ #include "Framework/EventGen/EventRecordVisitorI.h" +#include "Physics/HELepton/XSection/Born.h" #ifdef __GENIE_PYTHIA6_ENABLED__ #include @@ -51,7 +53,7 @@ public : mutable TPythia6 * fPythia; ///< PYTHIA6 wrapper class #endif - double fWmin; // Minimum value of W + Born * born; }; diff --git a/src/Physics/HELepton/EventGen/HELeptonInteractionListGenerator.cxx b/src/Physics/HELepton/EventGen/HELeptonInteractionListGenerator.cxx new file mode 100644 index 000000000..d20761797 --- /dev/null +++ b/src/Physics/HELepton/EventGen/HELeptonInteractionListGenerator.cxx @@ -0,0 +1,229 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include "Framework/EventGen/InteractionList.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Physics/HELepton/EventGen/HELeptonInteractionListGenerator.h" + +using namespace genie; + +//___________________________________________________________________________ +HELeptonInteractionListGenerator::HELeptonInteractionListGenerator() : +InteractionListGeneratorI("genie::HELeptonInteractionListGenerator") +{ + +} +//___________________________________________________________________________ +HELeptonInteractionListGenerator::HELeptonInteractionListGenerator(string config) : +InteractionListGeneratorI("genie::HELeptonInteractionListGenerator", config) +{ + +} +//___________________________________________________________________________ +HELeptonInteractionListGenerator::~HELeptonInteractionListGenerator() +{ + +} +//___________________________________________________________________________ +InteractionList * + HELeptonInteractionListGenerator::GLRESInteraction( + const InitialState & init_state) const +{ + + + InteractionList * intlist = new InteractionList; + + ProcessInfo proc_info(kScGlashowResonance, kIntWeakCC); + + int probepdg = init_state.ProbePdg(); + + if(probepdg == kPdgAntiNuE) { + InitialState init(init_state); + init_state.TgtPtr()->SetHitNucPdg(0); + Interaction * interaction = new Interaction(init_state, proc_info); + XclsTag exclusive_tag; + if (fIsGLRESMu) exclusive_tag.SetFinalLepton(kPdgMuon); + else if (fIsGLRESTau) exclusive_tag.SetFinalLepton(kPdgTau); + else if (fIsGLRESEle) exclusive_tag.SetFinalLepton(kPdgElectron); + else if (fIsGLRESHad) exclusive_tag.SetFinalLepton(kPdgPiP); + interaction->SetExclTag(exclusive_tag); + intlist->push_back(interaction); + } + + return intlist; + +} +//___________________________________________________________________________ +InteractionList * + HELeptonInteractionListGenerator::HENuElectronInteraction( + const InitialState & init_state) const +{ + + + InteractionList * intlist = new InteractionList; + + int probepdg = init_state.ProbePdg(); + + if (fIsHENuElCC) { + ProcessInfo proc_info(kScGlashowResonance, kIntWeakCC); + InitialState init(init_state); + init_state.TgtPtr()->SetHitNucPdg(0); + Interaction * interaction = new Interaction(init_state, proc_info); + XclsTag exclusive_tag; //charged lepton + if ( pdg::IsNuMu(probepdg) ) exclusive_tag.SetFinalLepton(kPdgMuon); + else if ( pdg::IsNuTau(probepdg) ) exclusive_tag.SetFinalLepton(kPdgTau); + else if ( pdg::IsNuE(probepdg) ) exclusive_tag.SetFinalLepton(kPdgElectron); + else return intlist; + interaction->SetExclTag(exclusive_tag); + intlist->push_back(interaction); + } + else if (fIsHENuElNC) { + ProcessInfo proc_info(kScGlashowResonance, kIntWeakNC); + InitialState init(init_state); + init_state.TgtPtr()->SetHitNucPdg(0); + Interaction * interaction = new Interaction(init_state, proc_info); + XclsTag exclusive_tag; //charged lepton + if ( pdg::IsNuMu(probepdg) ) exclusive_tag.SetFinalLepton(kPdgElectron); + else if ( pdg::IsNuTau(probepdg) ) exclusive_tag.SetFinalLepton(kPdgElectron); + else if ( pdg::IsAntiNuMu(probepdg) ) exclusive_tag.SetFinalLepton(kPdgElectron); + else if ( pdg::IsAntiNuTau(probepdg) ) exclusive_tag.SetFinalLepton(kPdgElectron); + else return intlist; + interaction->SetExclTag(exclusive_tag); + intlist->push_back(interaction); + } + + return intlist; + +} +//___________________________________________________________________________ +InteractionList * + HELeptonInteractionListGenerator::PhotonRESInteraction( + const InitialState & init_state) const +{ + + InteractionList * intlist = new InteractionList; + + ProcessInfo proc_info(kScPhotonResonance, kIntWeakCC); + + int probepdg = init_state.ProbePdg(); + bool hasP = (init_state.Tgt().Z() > 0); + bool hasN = (init_state.Tgt().N() > 0); + + int nuclpdg[2] = { kPdgProton, kPdgNeutron }; + for(int inucl=0; inucl<2; inucl++) { + int struck_nucleon = nuclpdg[inucl]; + if( (struck_nucleon == kPdgProton && hasP) || (struck_nucleon == kPdgNeutron && hasN) ) { + Interaction * interaction = new Interaction(init_state, proc_info); + Target * target = interaction->InitStatePtr()->TgtPtr(); + target->SetHitNucPdg(struck_nucleon); + XclsTag exclusive_tag; + if (fIsPhotonRESMu) exclusive_tag.SetFinalLepton( (probepdg>0) ? kPdgAntiMuon : kPdgMuon ); + else if (fIsPhotonRESTau) exclusive_tag.SetFinalLepton( (probepdg>0) ? kPdgAntiTau : kPdgTau ); + else if (fIsPhotonRESEle) exclusive_tag.SetFinalLepton( (probepdg>0) ? kPdgPositron : kPdgElectron ); + else if (fIsPhotonRESHad) exclusive_tag.SetFinalLepton( (probepdg>0) ? kPdgPiP : kPdgPiM ); + interaction->SetExclTag(exclusive_tag); + intlist->push_back(interaction); + } + } + + return intlist; + +} +//___________________________________________________________________________ +InteractionList * + HELeptonInteractionListGenerator::PhotonCOHInteraction( + const InitialState & init_state) const +{ + + InteractionList * intlist = new InteractionList; + ProcessInfo proc_info(kScPhotonCoherent, kIntWeakCC); + Interaction * interaction = new Interaction(init_state, proc_info); + intlist->push_back(interaction); + return intlist; + +} +//___________________________________________________________________________ +InteractionList * + HELeptonInteractionListGenerator::CreateInteractionList( + const InitialState & init_state) const +{ +// channels: +// nuebar + e- -> W- -> nuebar + e- [CC+NC] +// nuebar + e- -> W- -> nuebar + mu- [CC] +// nuebar + e- -> W- -> nuebar + tau- [CC] +// nuebar + e- -> W- -> hadrons [CC] +// nue + e- -> e + nue [CC+NC] +// numu + e- -> mu + nue [CC] +// nutau + e- -> tau + nue [CC] +// numu + e- -> numu + e [NC] +// nutau + e- -> nutau + e [NC] +// numubar + e- -> numubar + e [NC] +// nutaubar + e- -> nutaubar + e [NC] +// nu + gamma* -> l- + W+ (coherent & resonant) +// nubar + gamma* -> l+ + W- (coherent & resonant) + + int ppdg = init_state.ProbePdg(); + if( !pdg::IsNeutralLepton(ppdg) ) { + LOG("IntLst", pWARN) + << "Can not handle probe! Returning NULL InteractionList " + << "for init-state: " << init_state.AsString(); + return 0; + } + + if (fIsGLRESMu) return GLRESInteraction(init_state); + else if (fIsGLRESTau) return GLRESInteraction(init_state); + else if (fIsGLRESEle) return GLRESInteraction(init_state); + else if (fIsGLRESHad) return GLRESInteraction(init_state); + else if (fIsHENuElCC) return HENuElectronInteraction(init_state); + else if (fIsHENuElNC) return HENuElectronInteraction(init_state); + else if (fIsPhotonRESMu) return PhotonRESInteraction(init_state); + else if (fIsPhotonRESTau) return PhotonRESInteraction(init_state); + else if (fIsPhotonRESEle) return PhotonRESInteraction(init_state); + else if (fIsPhotonRESHad) return PhotonRESInteraction(init_state); + else if (fIsPhotonCOH) return PhotonCOHInteraction(init_state); + else { + LOG("IntLst", pERROR) + << "Returning NULL InteractionList for init-state: " << init_state.AsString(); + return 0; + } + +} +//___________________________________________________________________________ +void HELeptonInteractionListGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfigData(); +} +//____________________________________________________________________________ +void HELeptonInteractionListGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfigData(); +} +//____________________________________________________________________________ +void HELeptonInteractionListGenerator::LoadConfigData(void) +{ + + GetParamDef("is-GLRES-Mu", fIsGLRESMu, false ) ; + GetParamDef("is-GLRES-Tau", fIsGLRESTau, false ) ; + GetParamDef("is-GLRES-Ele", fIsGLRESEle, false ) ; + GetParamDef("is-GLRES-Had", fIsGLRESHad, false ) ; + GetParamDef("is-HENuEl-CC", fIsHENuElCC, false ) ; + GetParamDef("is-HENuEl-NC", fIsHENuElNC, false ) ; + GetParamDef("is-PhotonRES-Mu", fIsPhotonRESMu, false ) ; + GetParamDef("is-PhotonRES-Tau", fIsPhotonRESTau, false ) ; + GetParamDef("is-PhotonRES-Ele", fIsPhotonRESEle, false ) ; + GetParamDef("is-PhotonRES-Had", fIsPhotonRESHad, false ) ; + GetParamDef("is-PhotonCOH", fIsPhotonCOH, false ) ; + +} \ No newline at end of file diff --git a/src/Physics/HELepton/EventGen/HELeptonInteractionListGenerator.h b/src/Physics/HELepton/EventGen/HELeptonInteractionListGenerator.h new file mode 100644 index 000000000..85e143c41 --- /dev/null +++ b/src/Physics/HELepton/EventGen/HELeptonInteractionListGenerator.h @@ -0,0 +1,66 @@ +//____________________________________________________________________________ +/*! + +\class genie::HELeptonInteractionListGenerator + +\brief Interaction list generator in HELepton + +\author Alfonso Garcia + IFIC & Harvard University + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _HE_LEPTON_INTERACTION_GENERATOR_H_ +#define _HE_LEPTON_INTERACTION_GENERATOR_H_ + +#include "Framework/EventGen/InteractionListGeneratorI.h" + +namespace genie { + +class HELeptonInteractionListGenerator : public InteractionListGeneratorI { + +public : + HELeptonInteractionListGenerator(); + HELeptonInteractionListGenerator(string config); + ~HELeptonInteractionListGenerator(); + + // implement the InteractionListGeneratorI interface + InteractionList * CreateInteractionList(const InitialState & init) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + + InteractionList * GLRESInteraction (const InitialState & init_state) const; + InteractionList * HENuElectronInteraction (const InitialState & init_state) const; + InteractionList * PhotonRESInteraction (const InitialState & init_state) const; + InteractionList * PhotonCOHInteraction (const InitialState & init_state) const; + + void LoadConfigData(void); + + bool fIsGLRESMu; + bool fIsGLRESTau; + bool fIsGLRESEle; + bool fIsGLRESHad; + bool fIsHENuElCC; + bool fIsHENuElNC; + bool fIsPhotonRESMu; + bool fIsPhotonRESEle; + bool fIsPhotonRESTau; + bool fIsPhotonRESHad; + bool fIsPhotonCOH; + +}; + +} // genie namespace + +#endif // _HE_LEPTON_INTERACTION_GENERATOR_H_ diff --git a/src/Physics/HELepton/EventGen/HELeptonKinematicsGenerator.cxx b/src/Physics/HELepton/EventGen/HELeptonKinematicsGenerator.cxx new file mode 100644 index 000000000..78fdc2849 --- /dev/null +++ b/src/Physics/HELepton/EventGen/HELeptonKinematicsGenerator.cxx @@ -0,0 +1,356 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include "Physics/HELepton/EventGen/HELeptonKinematicsGenerator.h" +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/KinePhaseSpace.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/EventGen/EventGeneratorI.h" +#include "Framework/EventGen/RunningThreadInfo.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/Utils/KineUtils.h" +#include "Framework/Utils/Range1.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Physics/XSectionIntegration/GSLXSecFunc.h" + +#include "Math/Minimizer.h" +#include "Math/Factory.h" + +using namespace genie; +using namespace genie::controls; +using namespace genie::utils; + +//___________________________________________________________________________ +HELeptonKinematicsGenerator::HELeptonKinematicsGenerator() : +KineGeneratorWithCache("genie::HELeptonKinematicsGenerator") +{ + +} +//___________________________________________________________________________ +HELeptonKinematicsGenerator::HELeptonKinematicsGenerator(string config) : +KineGeneratorWithCache("genie::HELeptonKinematicsGenerator", config) +{ + +} +//___________________________________________________________________________ +HELeptonKinematicsGenerator::~HELeptonKinematicsGenerator() +{ + +} +//___________________________________________________________________________ +void HELeptonKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ + if(fGenerateUniformly) { + LOG("HELeptonKinematics", pNOTICE) + << "Generating kinematics uniformly over the allowed phase space"; + } + + //-- Get the random number generators + RandomGen * rnd = RandomGen::Instance(); + + //-- Access cross section algorithm for running thread + RunningThreadInfo * rtinfo = RunningThreadInfo::Instance(); + const EventGeneratorI * evg = rtinfo->RunningThread(); + fXSecModel = evg->CrossSectionAlg(); + + Interaction * interaction = evrec->Summary(); + + //-- For the subsequent kinematic selection with the rejection method: + // Calculate the max differential cross section or retrieve it from the + // cache. Throw an exception and quit the evg thread if a non-positive + // value is found. + // If the kinematics are generated uniformly over the allowed phase + // space the max xsec is irrelevant + double xsec_max = this->MaxXSec(evrec); + + const ProcessInfo & proc_info = interaction->ProcInfo(); + if(proc_info.IsPhotonCoherent()) { + + double nupdg = interaction->InitState().ProbePdg(); + + double n2min = 0.; + double n2max = 1.; + double n3min = 0.; + double n3max = 1.; + double dn2 = n2max-n2min; + double dn3 = n3max-n3min; + + double n1max = 0.; + double n1min = 0.; + if (pdg::IsNuE (TMath::Abs(nupdg))) { n1min = 1.-fDeltaN1NuE; n1max = 1.+fDeltaN1NuE; } + else if (pdg::IsNuMu (TMath::Abs(nupdg))) { n1min = 1.-fDeltaN1NuMu; n1max = 1.+fDeltaN1NuMu; } + else if (pdg::IsNuTau(TMath::Abs(nupdg))) { n1min = 1.-fDeltaN1NuTau; n1max = 1.+fDeltaN1NuTau; } + double dn1 = n1max-n1min; + + + //-- Try to select a valid inelastisity y + double xsec = -1; + unsigned int iter = 0; + bool accept = false; + + while(1) { + iter++; + if(iter > 1000000) { + LOG("HELeptonKinematics", pWARN) + << "*** Could not select a valid y after " + << iter << " iterations"; + evrec->EventFlags()->SetBitNumber(kKineGenErr, true); + genie::exceptions::EVGThreadException exception; + exception.SetReason("Couldn't select kinematics"); + exception.SwitchOnFastForward(); + throw exception; + } + + double n2 = n2min + dn2 * rnd->RndKine().Rndm(); + double n3 = n3min + dn3 * rnd->RndKine().Rndm(); + double n1 = n1min + dn1 * rnd->RndKine().Rndm(); + n1 = (n1>1.) ? n1-2. : n1; + + interaction->KinePtr()->SetKV(kKVn1,n1); + interaction->KinePtr()->SetKV(kKVn2,n2); + interaction->KinePtr()->SetKV(kKVn3,n3); + + LOG("HELeptonKinematics", pDEBUG) << "Trying: n1 = " << n1 << ", n2 = " << n2 << ", n3 = " << n3; + + //-- computing cross section for the current kinematics + xsec = fXSecModel->XSec(interaction, kPSn1n2n3fE); + + this->AssertXSecLimits(interaction, xsec, xsec_max); + + double t = xsec_max * rnd->RndKine().Rndm(); + LOG("HELeptonKinematics", pDEBUG) << "xsec= "<< xsec<< ", J= 1, Rnd= "<< t; + + accept = (tEventFlags()->SetBitNumber(kKineGenErr, true); + genie::exceptions::EVGThreadException exception; + exception.SetReason("Couldn't select kinematics"); + exception.SwitchOnFastForward(); + throw exception; + } + + double n1 = n1min + dn1 * rnd->RndKine().Rndm(); + double n2 = n2min + dn2 * rnd->RndKine().Rndm(); + interaction->KinePtr()->SetKV(kKVn1,n1); + interaction->KinePtr()->SetKV(kKVn2,n2); + + LOG("HELeptonKinematics", pDEBUG) << "Trying: n1 = " << n1 << ", n2 = " << n2; + + //-- computing cross section for the current kinematics + xsec = fXSecModel->XSec(interaction, kPSn1n2fE); + + this->AssertXSecLimits(interaction, xsec, xsec_max); + + double t = xsec_max * rnd->RndKine().Rndm(); + LOG("HELeptonKinematics", pDEBUG) << "xsec= "<< xsec<< ", J= 1, Rnd= "<< t; + + accept = (tProcInfo(); + if(proc_info.IsPhotonCoherent()) { + + utils::gsl::d2Xsec_dn1dn2dn3_E f(fXSecModel,interaction,-1); + min->SetFunction( f ); + min->SetMaxFunctionCalls(10000); + min->SetTolerance(0.05); + + min->SetFixedVariable ( 0, "n1", 1.); + min->SetLimitedVariable ( 1, "n2", 0., 0.01, 0., 1.); + min->SetLimitedVariable ( 2, "n3", 0., 0.01, 0., 1.); + min->Minimize(); + min->SetFixedVariable ( 0, "n1", 1.); + min->SetLimitedVariable ( 1, "n2", min->X()[1], 0.01, TMath::Max(0.,min->X()[1]-0.1), TMath::Min(1.,min->X()[1]+0.1)); + min->SetLimitedVariable ( 2, "n3", min->X()[2], 0.01, TMath::Max(0.,min->X()[2]-0.1), TMath::Min(1.,min->X()[2]+0.1)); + min->Minimize(); + interaction->KinePtr()->SetKV(kKVn1,min->X()[0]); + interaction->KinePtr()->SetKV(kKVn2,min->X()[1]); + interaction->KinePtr()->SetKV(kKVn3,min->X()[2]); + SLOG("HELeptonKinematics", pDEBUG) << "Minimum found -> n1: 1, n2: " << min->X()[1] << ", n3: " << min->X()[2] << ", xsec: " << fXSecModel->XSec(interaction, kPSn1n2n3fE); + max_xsec = TMath::Max(fXSecModel->XSec(interaction, kPSn1n2n3fE),max_xsec); + + min->SetFixedVariable ( 0, "n1", -1.); + min->SetLimitedVariable ( 1, "n2", 0., 0.01, 0., 1.); + min->SetLimitedVariable ( 2, "n3", 0., 0.01, 0., 1.); + min->Minimize(); + min->SetFixedVariable ( 0, "n1", -1.); + min->SetLimitedVariable ( 1, "n2", min->X()[1], 0.01, TMath::Max(0.,min->X()[1]-0.1), TMath::Min(1.,min->X()[1]+0.1)); + min->SetLimitedVariable ( 2, "n3", min->X()[2], 0.01, TMath::Max(0.,min->X()[2]-0.1), TMath::Min(1.,min->X()[2]+0.1)); + interaction->KinePtr()->SetKV(kKVn1,min->X()[0]); + interaction->KinePtr()->SetKV(kKVn2,min->X()[1]); + interaction->KinePtr()->SetKV(kKVn3,min->X()[2]); + SLOG("HELeptonKinematics", pDEBUG) << "Minimum found -> n1: -1, n2: " << min->X()[1] << ", n3: " << min->X()[2] << ", xsec: " << fXSecModel->XSec(interaction, kPSn1n2n3fE); + max_xsec = TMath::Max(fXSecModel->XSec(interaction, kPSn1n2n3fE),max_xsec); + + } + else { + + const int Nscan = 100; + const int n1min = -1.; + const int n1max = 1.; + const int n2min = 0.; + const int n2max = 1.; + const double dn1 = (n1max-n1min)/(double)Nscan; + const double dn2 = (n2max-n2min)/(double) + Nscan; + + double scan_n1 = 0.; + double scan_n2 = 0.; + for (int i=0; iKinePtr()->SetKV(kKVn1,n1); + interaction->KinePtr()->SetKV(kKVn2,n2); + double dxsec = fXSecModel->XSec(interaction, kPSn1n2fE); + if ( dxsec > max_xsec ) { + scan_n1 = n1; + scan_n2 = n2; + max_xsec = dxsec; + } + } + } + + utils::gsl::d2Xsec_dn1dn2_E f(fXSecModel,interaction,-1); + min->SetFunction( f ); + min->SetMaxFunctionCalls(10000); + min->SetTolerance(0.05); + min->SetLimitedVariable ( 0, "n1", scan_n1, 0.001, TMath::Max(-1.,scan_n1-0.1), TMath::Min(1.,scan_n1+0.1)); + min->SetLimitedVariable ( 1, "n2", scan_n2, 0.1, TMath::Max(-0.,scan_n2-0.1), TMath::Min(1.,scan_n2+0.1)); + min->Minimize(); + interaction->KinePtr()->SetKV(kKVn1,min->X()[0]); + interaction->KinePtr()->SetKV(kKVn2,min->X()[1]); + max_xsec = fXSecModel->XSec(interaction, kPSn1n2fE); + SLOG("HELeptonKinematics", pDEBUG) << "Minimum found -> n1: " << min->X()[0] << ", n2: " << min->X()[1]; + + } + + // Apply safety factor, since value retrieved from the cache might + // correspond to a slightly different energy. + max_xsec *= fSafetyFactor; + + SLOG("HELeptonKinematics", pDEBUG) << interaction->AsString(); + SLOG("HELeptonKinematics", pDEBUG) << "Max xsec in phase space = " << max_xsec; + SLOG("HELeptonKinematics", pDEBUG) << "Computed using alg = " << *fXSecModel; + + return max_xsec; +} +//___________________________________________________________________________ +double HELeptonKinematicsGenerator::Energy(const Interaction * interaction) const +{ +// Override the base class Energy() method to cache the max xsec for the +// neutrino energy in the LAB rather than in the hit nucleon rest frame. + + const InitialState & init_state = interaction->InitState(); + double E = init_state.ProbeE(kRfLab); + return E; +} +//___________________________________________________________________________ +void HELeptonKinematicsGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void HELeptonKinematicsGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void HELeptonKinematicsGenerator::LoadConfig(void) +{ +// Reads its configuration data from its configuration Registry and loads them +// in private data members to avoid looking up at the Registry all the time. + + //-- Safety factor for the maximum differential cross section + GetParamDef( "MaxXSec-SafetyFactor", fSafetyFactor, 2. ) ; + + //-- Maximum allowed fractional cross section deviation from maxim cross + // section used in rejection method + GetParamDef( "MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999. ) ; + assert(fMaxXSecDiffTolerance>=0); + + //-- Sampling range for n1 variable + GetParamDef( "Delta-N1-NuE", fDeltaN1NuE, 0. ) ; + GetParamDef( "Delta-N1-NuMu", fDeltaN1NuMu, 0. ) ; + GetParamDef( "Delta-N1-NuTau", fDeltaN1NuTau, 0. ) ; + + +} +//____________________________________________________________________________ diff --git a/src/Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.h b/src/Physics/HELepton/EventGen/HELeptonKinematicsGenerator.h similarity index 51% rename from src/Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.h rename to src/Physics/HELepton/EventGen/HELeptonKinematicsGenerator.h index 720d1a353..66c3fac77 100644 --- a/src/Physics/GlashowResonance/EventGen/GLRESKinematicsGenerator.h +++ b/src/Physics/HELepton/EventGen/HELeptonKinematicsGenerator.h @@ -1,37 +1,34 @@ //____________________________________________________________________________ /*! -\class genie::GLRESKinematicsGenerator -\brief Generates values for the kinematic variables describing Glashow resonance. - Is a concrete implementation of the EventRecordVisitorI interface. - Part of its implementation, related with the caching and retrieval of - previously computed values, is inherited from the KineGeneratorWithCache - abstract class. +\class genie::HELeptonKinematicsGenerator -\author Alfonso Garcia - NIKHEF (Amsterdam) +\brief Kinematics generator for HELepton. -\created November 8, 2019 +\author Alfonso Garcia + IFIC & Harvard University -\cpright Copyright (c) 2003-2019, The GENIE Collaboration +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org or see $GENIE/LICENSE */ //____________________________________________________________________________ -#ifndef _GLASHOW_RESONANCE_KINEMATICS_GENERATOR_H_ -#define _GLASHOW_RESONANCE_KINEMATICS_GENERATOR_H_ +#ifndef _HE_LEPTON_KINEMATICS_GENERATOR_H_ +#define _HE_LEPTON_KINEMATICS_GENERATOR_H_ #include "Physics/Common/KineGeneratorWithCache.h" namespace genie { -class GLRESKinematicsGenerator : public KineGeneratorWithCache { +class HELeptonKinematicsGenerator : public KineGeneratorWithCache { public : - GLRESKinematicsGenerator(); - GLRESKinematicsGenerator(string config); - ~GLRESKinematicsGenerator(); + HELeptonKinematicsGenerator(); + HELeptonKinematicsGenerator(string config); + ~HELeptonKinematicsGenerator(); //-- implement the EventRecordVisitorI interface void ProcessEventRecord(GHepRecord * event_rec) const; @@ -41,7 +38,7 @@ public : void Configure(const Registry & config); void Configure(string config); -public: +private: //-- methods to load sub-algorithms and config data from the Registry void LoadConfig (void); @@ -49,7 +46,12 @@ public : //-- overload KineGeneratorWithCache methods double ComputeMaxXSec (const Interaction * in) const; double Energy (const Interaction * in) const; + + double fDeltaN1NuE; + double fDeltaN1NuMu; + double fDeltaN1NuTau; + }; } // genie namespace -#endif // _GLASHOW_RESONANCE_KINEMATICS_GENERATOR_H_ +#endif // _HE_LEPTON_KINEMATICS_GENERATOR_H_ diff --git a/src/Physics/HELepton/EventGen/HENuElGenerator.cxx b/src/Physics/HELepton/EventGen/HENuElGenerator.cxx new file mode 100644 index 000000000..5a94a4223 --- /dev/null +++ b/src/Physics/HELepton/EventGen/HENuElGenerator.cxx @@ -0,0 +1,152 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include + +#include +#include +#include + +#include "Framework/Conventions/Constants.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Physics/HELepton/EventGen/HENuElGenerator.h" + +#ifdef __GENIE_PYTHIA6_ENABLED__ +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) +#include +#else +#include +#endif +#endif // __GENIE_PYTHIA6_ENABLED__ + +using namespace genie; +using namespace genie::constants; +using namespace genie::utils::math; + +//___________________________________________________________________________ +HENuElGenerator::HENuElGenerator() : +EventRecordVisitorI("genie::HENuElGenerator") +{ + born = new Born(); +} +//___________________________________________________________________________ +HENuElGenerator::HENuElGenerator(string config) : +EventRecordVisitorI("genie::HENuElGenerator", config) +{ + +} +//___________________________________________________________________________ +HENuElGenerator::~HENuElGenerator() +{ + +} +//___________________________________________________________________________ +void HENuElGenerator::ProcessEventRecord(GHepRecord * event) const +{ + + Interaction * interaction = event->Summary(); + const InitialState & init_state = interaction->InitState(); + const ProcessInfo & proc_info = interaction->ProcInfo(); + + //incoming v & struck particle & remnant nucleus + GHepParticle * nu = event->Probe(); + + GHepParticle * target = event -> TargetNucleus(); + if(target) event->AddParticle(target->Pdg(), kIStFinalStateNuclearRemnant, 1,-1,-1,-1, *(target->P4()), *(target->X4()) ); + + TVector3 unit_nu = nu->P4()->Vect().Unit(); + + bool isCC = proc_info.IsWeakCC(); + + long double mlout = interaction->FSPrimLepton()->Mass(); //mass of charged lepton + long double mlin = kElectronMass; //mass of incoming charged lepton + + long double Enuin = init_state.ProbeE(kRfLab); + long double s = born->GetS(mlin,Enuin); + + long double n1 = interaction->Kine().GetKV(kKVn1); + long double n2 = interaction->Kine().GetKV(kKVn2); + + long double costhCM = n1; + long double sinthCM = sqrtl(1-costhCM*costhCM); + + long double t = born->GetT( mlin, mlout, s, n1 ); + long double zeta = born->GetReAlpha()/kPi*(2.0*logl(sqrtl(-t)/kElectronMass)-1.0); + long double omx = powl(n2, 1.0/zeta ); + long double s_r = s*( 1.-omx ); + + // Boost velocity CM -> LAB + long double EnuinCM = (s_r-mlin*mlin)/sqrtl(s_r)/2.; + long double beta = (powl(Enuin,2)-powl(EnuinCM,2))/(powl(Enuin,2)+powl(EnuinCM,2)); + + // Final state primary lepton PDG code + int pdgl = interaction->FSPrimLeptonPdg(); + assert(pdgl!=0); + + long double ElpoutCM = (s_r+mlout*mlout)/sqrtl(s_r)/2.; + long double EnuoutCM = (s_r-mlout*mlout)/sqrtl(s_r)/2.; + LongLorentzVector p4_lpout( 0., EnuoutCM*sinthCM, EnuoutCM*costhCM, ElpoutCM ); + LongLorentzVector p4_nuout( 0., -EnuoutCM*sinthCM, -EnuoutCM*costhCM, EnuoutCM ); + + p4_lpout.BoostZ(beta); + p4_nuout.BoostZ(beta); + + TLorentzVector p4lp_o( (double)p4_lpout.Px(), (double)p4_lpout.Py(), (double)p4_lpout.Pz(), (double)p4_lpout.E() ); + TLorentzVector p4nu_o( (double)p4_nuout.Px(), (double)p4_nuout.Py(), (double)p4_nuout.Pz(), (double)p4_nuout.E() ); + + // Randomize transverse components + RandomGen * rnd = RandomGen::Instance(); + double phi = 2* kPi * rnd->RndLep().Rndm(); + p4lp_o.RotateZ(phi); + p4nu_o.RotateZ(phi); + + //rotate from LAB=[0,0,Ev,Ev]->[px,py,pz,E] + p4lp_o.RotateUz(unit_nu); + p4nu_o.RotateUz(unit_nu); + + int pdgvout = 0; + if (isCC) pdgvout = kPdgNuE; + else pdgvout = nu->Pdg(); + + // Create a GHepParticle and add it to the event record + event->AddParticle( pdgl, kIStStableFinalState, 4, -1, -1, -1, p4lp_o, *(nu->X4()) ); + event->AddParticle( pdgvout, kIStStableFinalState, 4, -1, -1, -1, p4nu_o, *(nu->X4()) ); + event->Summary()->KinePtr()->SetFSLeptonP4(p4lp_o); + +} +//___________________________________________________________________________ +void HENuElGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void HENuElGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void HENuElGenerator::LoadConfig(void) +{ + + +} +//____________________________________________________________________________ diff --git a/src/Physics/HELepton/EventGen/HENuElGenerator.h b/src/Physics/HELepton/EventGen/HENuElGenerator.h new file mode 100644 index 000000000..523f12792 --- /dev/null +++ b/src/Physics/HELepton/EventGen/HENuElGenerator.h @@ -0,0 +1,53 @@ +//____________________________________________________________________________ +/*! + +\class genie::HENuElGenerator + +\brief Generator for high energy neutrino-electron scattering. + +\author Alfonso Garcia + IFIC & Harvard University + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _HE_NUEL_GENERATOR_H_ +#define _HE_NUEL_GENERATOR_H_ + +#define __GENIE_PYTHIA6_ENABLED__ + +#include "Framework/EventGen/EventRecordVisitorI.h" +#include "Physics/HELepton/XSection/Born.h" + +namespace genie { + +class HENuElGenerator : public EventRecordVisitorI { + +public : + HENuElGenerator(); + HENuElGenerator(string config); + ~HENuElGenerator(); + + // implement the EventRecordVisitorI interface + void ProcessEventRecord (GHepRecord * event) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + + void LoadConfig(void); + + Born * born; + +}; + +} // genie namespace +#endif // _GLASHOW_RESONANCE_GENERATOR_H_ diff --git a/src/Physics/HELepton/EventGen/LinkDef.h b/src/Physics/HELepton/EventGen/LinkDef.h new file mode 100644 index 000000000..85fbf142d --- /dev/null +++ b/src/Physics/HELepton/EventGen/LinkDef.h @@ -0,0 +1,17 @@ +#ifdef __CINT__ + +#pragma link off all globals; +#pragma link off all classes; +#pragma link off all functions; + +#pragma link C++ namespace genie; + +#pragma link C++ class genie::HELeptonInteractionListGenerator; +#pragma link C++ class genie::HELeptonKinematicsGenerator; + +#pragma link C++ class genie::GLRESGenerator; +#pragma link C++ class genie::HENuElGenerator; +#pragma link C++ class genie::PhotonRESGenerator; +#pragma link C++ class genie::PhotonCOHGenerator; + +#endif diff --git a/src/Physics/GlashowResonance/XSection/Makefile b/src/Physics/HELepton/EventGen/Makefile similarity index 86% rename from src/Physics/GlashowResonance/XSection/Makefile rename to src/Physics/HELepton/EventGen/Makefile index 4b8c35f8b..555200ac8 100644 --- a/src/Physics/GlashowResonance/XSection/Makefile +++ b/src/Physics/HELepton/EventGen/Makefile @@ -12,8 +12,8 @@ MAKEFILE = Makefile # include $(GENIE)/src/make/Make.include -PACKAGE = Physics/GlashowResonance/XSection -PACKAGE_ABBREV = PhGlwResXS +PACKAGE = Physics/HELepton/EventGen +PACKAGE_ABBREV = PhHELptnEG DICTIONARY = _ROOT_DICT_$(PACKAGE_ABBREV) LIBNAME = libG$(PACKAGE_ABBREV) EXTRA_EXT_LIBS = diff --git a/src/Physics/HELepton/EventGen/PhotonCOHGenerator.cxx b/src/Physics/HELepton/EventGen/PhotonCOHGenerator.cxx new file mode 100644 index 000000000..85171b623 --- /dev/null +++ b/src/Physics/HELepton/EventGen/PhotonCOHGenerator.cxx @@ -0,0 +1,251 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include + +#include +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) +#include +#else +#include +#endif +#include +#include + +#include "Physics/HELepton/EventGen/PhotonCOHGenerator.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" + +using namespace genie; +using namespace genie::constants; +using namespace genie::utils::math; + +//___________________________________________________________________________ +PhotonCOHGenerator::PhotonCOHGenerator() : +EventRecordVisitorI("genie::PhotonCOHGenerator") +{ + this->Initialize(); +} +//___________________________________________________________________________ +PhotonCOHGenerator::PhotonCOHGenerator(string config) : +EventRecordVisitorI("genie::PhotonCOHGenerator", config) +{ + this->Initialize(); +} +//___________________________________________________________________________ +PhotonCOHGenerator::~PhotonCOHGenerator() +{ + +} +//____________________________________________________________________________ +void PhotonCOHGenerator::Initialize(void) const +{ + fPythia = TPythia6::Instance(); + + // sync GENIE/PYTHIA6 seed number + RandomGen::Instance(); +} +//___________________________________________________________________________ +void PhotonCOHGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ + + Interaction * interaction = evrec->Summary(); + const InitialState & init_state = interaction->InitState(); + + //incoming v & struck particle & remnant nucleus + GHepParticle * nu = evrec->Probe(); + + GHepParticle * target = evrec -> TargetNucleus(); + if(target) evrec->AddParticle(target->Pdg(), kIStFinalStateNuclearRemnant, 1,-1,-1,-1, *(target->P4()), *(target->X4()) ); + + TVector3 unit_nu = nu->P4()->Vect().Unit(); + + long double Ev = init_state.ProbeE(kRfLab); + + long double Mtgt = init_state.Tgt().Z()*kProtonMass + init_state.Tgt().N()*kNeutronMass; + + long double n1 = interaction->Kine().GetKV(kKVn1); + long double n2 = interaction->Kine().GetKV(kKVn2); + long double n3 = interaction->Kine().GetKV(kKVn3); + + long double costh = n1; + long double sinth = sqrtl(1-costh*costh); + + long double mlout = 0; + if (pdg::IsNuE (TMath::Abs(nu->Pdg()))) mlout = kElectronMass; + else if (pdg::IsNuMu (TMath::Abs(nu->Pdg()))) mlout = kMuonMass; + else if (pdg::IsNuTau(TMath::Abs(nu->Pdg()))) mlout = kTauMass; + long double mlout2 = mlout*mlout; + + long double mL = mlout+kMw; + long double Delta = sqrtl( powl(2.*Ev*Mtgt-mL*mL,2)-4.*powl(Mtgt*mL,2) ); + long double s12_min = Ev/(2.*Ev+Mtgt)*(mL*mL+2.*Ev*Mtgt-Delta); + long double s12_max = Ev/(2.*Ev+Mtgt)*(mL*mL+2.*Ev*Mtgt+Delta); + long double s12 = expl( logl(s12_min)+(logl(s12_max)-logl(s12_min))*n2); + long double Q2_min = powl(s12,2)*Mtgt/(2.*Ev*(2.*Ev*Mtgt-s12)); + long double Q2_max = s12 - mL*mL; + double Q2 = expl( logl(Q2_min) + (logl(Q2_max)-logl(Q2_min))*n3 ); + double s_r = s12 - Q2; + + long double EW = (s_r+kMw2-mlout2)/sqrtl(s_r)/2.; + long double El = (s_r-kMw2+mlout2)/sqrtl(s_r)/2.; + long double p = sqrtl( EW*EW - kMw2 ); + LongLorentzVector p4_lout( 0., -p*sinth, -p*costh, El ); + + long double bz = 4.*Ev*Mtgt*Q2/(Q2+s_r)/(2.*Ev*Mtgt-Q2) - (2.*Ev*Mtgt+Q2)/(2.*Ev*Mtgt-Q2); + long double by = sqrtl(Mtgt*powl(Q2+s_r,2)/(2.*Ev*Q2*(s_r+Q2-2.*Ev*Mtgt))+1.); + + p4_lout.BoostZ(-bz); + p4_lout.BoostY(-by); + + TLorentzVector p4l_o(p4_lout.Px(),p4_lout.Py(),p4_lout.Pz(),p4_lout.E()); + p4l_o.RotateX((double)acosl(by)-kPi/2.); + + + // Randomize transverse components + RandomGen * rnd = RandomGen::Instance(); + double phi = 2 * kPi * rnd->RndLep().Rndm(); + + p4l_o.RotateZ(phi); + + //rotate from LAB=[0,0,Ev,Ev]->[px,py,pz,E] + p4l_o.RotateUz(unit_nu); + + int pdglout = 0; + if (pdg::IsAntiNuE (nu->Pdg())) pdglout = kPdgPositron; + else if (pdg::IsNuE (nu->Pdg())) pdglout = kPdgElectron; + else if (pdg::IsAntiNuMu (nu->Pdg())) pdglout = kPdgAntiMuon; + else if (pdg::IsNuMu (nu->Pdg())) pdglout = kPdgMuon; + else if (pdg::IsAntiNuTau(nu->Pdg())) pdglout = kPdgAntiTau; + else if (pdg::IsNuTau (nu->Pdg())) pdglout = kPdgTau; + + // Create a GHepParticle and add it to the event record + evrec->AddParticle( pdglout, kIStStableFinalState, 0, -1, -1, -1, p4l_o, *(nu->X4()) ); + + int pdgboson = pdg::IsNeutrino(init_state.ProbePdg()) ? kPdgWP : kPdgWM; + + int def61 = fPythia->GetMSTP(61); + int def71 = fPythia->GetMSTP(71); + fPythia->SetMSTP(61,0); // (Default=2) master switch for initial-state QCD and QED radiation. + fPythia->SetMSTP(71,0); // (Default=2) master switch for initial-state QCD and QED radiation. + + fPythia->Py1ent( -1, pdgboson, EW, acosl(costh), kPi/2. ); //k(1,2) = 2 + fPythia->Pyexec(); + + fPythia->SetMSTP(61,def61); + fPythia->SetMSTP(71,def71); + + fPythia->GetPrimaries(); + TClonesArray * pythia_particles = (TClonesArray *) fPythia->ImportParticles("All"); + int np = pythia_particles->GetEntries(); + assert(np>0); + + TMCParticle * particle = 0; + TIter piter(pythia_particles); + while( (particle = (TMCParticle *) piter.Next()) ) { + + int pdgc = particle->GetKF(); + int ks = particle->GetKS(); + + LongLorentzVector p4longo(particle->GetPx(), particle->GetPy(), particle->GetPz(), particle->GetEnergy()); + p4longo.BoostZ(-bz); + p4longo.BoostY(-by); + + TLorentzVector p4o(p4longo.Px(),p4longo.Py(),p4longo.Pz(),p4longo.E()); + p4o.RotateX((double)acosl(by)-kPi/2.); + p4o.RotateZ(phi); + p4o.RotateUz(unit_nu); + + TParticlePDG * part = PDGLibrary::Instance()->Find(pdgc); + if ( (ks==1 || ks==4) && p4o.E() < part->Mass() ) { + LOG("PhotonCOHGenerator", pWARN) << "Putting at rest one stable particle generated by PYTHIA because E < m"; + LOG("PhotonCOHGenerator", pWARN) << "PDG = " << pdgc << " // State = " << ks; + LOG("PhotonCOHGenerator", pWARN) << "E = " << p4o.E() << " // |p| = " << TMath::Sqrt(p4o.P()); + LOG("PhotonCOHGenerator", pWARN) << "p = [ " << p4o.Px() << " , " << p4o.Py() << " , " << p4o.Pz() << " ]"; + LOG("PhotonCOHGenerator", pWARN) << "m = " << p4o.M() << " // mpdg = " << part->Mass(); + p4o.SetXYZT(0,0,0,part->Mass()); + } + + // copy final state particles to the event record + GHepStatus_t ist = (ks==1 || ks==4) ? kIStStableFinalState : kIStDISPreFragmHadronicState; + + // fix numbering scheme used for mother/daughter assignments + int firstmother = -1; + int lastmother = -1; + int firstchild = -1; + int lastchild = -1; + + if (particle->GetParent()==0) { + firstmother = 0; + } + else { + firstmother = particle->GetParent() + 3; + if (particle->GetFirstChild()!=0) firstchild = particle->GetFirstChild() + 3; + if (particle->GetLastChild() !=0) lastchild = particle->GetLastChild() + 3; + + } + + double lightspeed = 299792458e3; //c in mm/s. Used for time in PYTHIA t[s]=t_pythia[mm]/c[mm/s] + double vx = nu->X4()->X() + particle->GetVx()*1e12; //pythia gives position in [mm] while genie uses [fm] + double vy = nu->X4()->Y() + particle->GetVy()*1e12; + double vz = nu->X4()->Z() + particle->GetVz()*1e12; + double vt = nu->X4()->T() + particle->GetTime()/lightspeed; + TLorentzVector pos( vx, vy, vz, vt ); + + evrec->AddParticle(pdgc, ist, firstmother, lastmother, firstchild, lastchild, p4o, pos ); + + } + + delete particle; + pythia_particles->Clear("C"); + +} +//___________________________________________________________________________ +void PhotonCOHGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonCOHGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonCOHGenerator::LoadConfig(void) +{ + + // PYTHIA parameters only valid for HEDIS & HELepton + double wmin; GetParam( "Xsec-Wmin", wmin ) ; + int warnings; GetParam( "Warnings", warnings ) ; + int errors; GetParam( "Errors", errors ) ; + int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ; + fPythia->SetPARP(2, wmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes) + fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed + fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed + fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term. + + fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385) + fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation + fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation + +} +//____________________________________________________________________________ \ No newline at end of file diff --git a/src/Physics/HELepton/EventGen/PhotonCOHGenerator.h b/src/Physics/HELepton/EventGen/PhotonCOHGenerator.h new file mode 100644 index 000000000..dc68593d7 --- /dev/null +++ b/src/Physics/HELepton/EventGen/PhotonCOHGenerator.h @@ -0,0 +1,55 @@ +//____________________________________________________________________________ +/*! + +\class genie::PhotonCOHGenerator + +\brief Generator for W boson production. + +\author Alfonso Garcia + IFIC & Harvard University + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _PHOTON_COH_GENERATOR_H_ +#define _PHOTON_COH_GENERATOR_H_ + +#include +#include + +#include "Framework/EventGen/EventRecordVisitorI.h" + +namespace genie { + +class PhotonCOHGenerator : public EventRecordVisitorI { + +public : + PhotonCOHGenerator(); + PhotonCOHGenerator(string config); + ~PhotonCOHGenerator(); + + // implement the EventRecordVisitorI interface + void Initialize (void) const; + void ProcessEventRecord (GHepRecord * evrec) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + + void Configure(string config); + +private: + + void LoadConfig(void); + + mutable TPythia6 * fPythia; ///< PYTHIA6 wrapper class + +}; + +} // genie namespace +#endif // _PHOTON_COH_GENERATOR_H_ \ No newline at end of file diff --git a/src/Physics/HELepton/EventGen/PhotonRESGenerator.cxx b/src/Physics/HELepton/EventGen/PhotonRESGenerator.cxx new file mode 100644 index 000000000..c10abee30 --- /dev/null +++ b/src/Physics/HELepton/EventGen/PhotonRESGenerator.cxx @@ -0,0 +1,289 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include + +#include +#if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) +#include +#else +#include +#endif +#include +#include + +#include "Physics/HELepton/EventGen/PhotonRESGenerator.h" +#include "Physics/HEDIS/EventGen/HEDISGenerator.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/RandomGen.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" + +using namespace genie; +using namespace genie::constants; +using namespace genie::utils::math; + +//___________________________________________________________________________ +PhotonRESGenerator::PhotonRESGenerator() : +EventRecordVisitorI("genie::PhotonRESGenerator") +{ + this->Initialize(); + born = new Born(); +} +//___________________________________________________________________________ +PhotonRESGenerator::PhotonRESGenerator(string config) : +EventRecordVisitorI("genie::PhotonRESGenerator", config) +{ + this->Initialize(); +} +//___________________________________________________________________________ +PhotonRESGenerator::~PhotonRESGenerator() +{ + +} +//____________________________________________________________________________ +void PhotonRESGenerator::Initialize(void) const +{ + fPythia = TPythia6::Instance(); + + // sync GENIE/PYTHIA6 seed number + RandomGen::Instance(); +} +//___________________________________________________________________________ +void PhotonRESGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ + + Interaction * interaction = evrec->Summary(); + const InitialState & init_state = interaction->InitState(); + + //incoming v & struck particle & remnant nucleus + GHepParticle * nu = evrec->Probe(); + GHepParticle * el = evrec->HitNucleon(); + + GHepParticle * target = evrec -> TargetNucleus(); + if(target) evrec->AddParticle(target->Pdg(), kIStFinalStateNuclearRemnant, 1,-1,-1,-1, *(target->P4()), *(target->X4()) ); + + TVector3 unit_nu = nu->P4()->Vect().Unit(); + + int probepdg = init_state.ProbePdg(); + + long double Mtarget = init_state.Tgt().HitNucMass(); + long double mlout = interaction->FSPrimLepton()->Mass(); //mass of charged lepton + + long double Enuin = init_state.ProbeE(kRfLab); + long double s = born->GetS(Mtarget,Enuin); + + long double n1 = interaction->Kine().GetKV(kKVn1); + long double n2 = interaction->Kine().GetKV(kKVn2); + + long double costhCM = n1; + long double sinthCM = sqrtl(1-costhCM*costhCM); + + long double xmin = fQ2PDFmin/2./Enuin/Mtarget; + long double x = expl( logl(xmin) + (logl(1.0)-logl(xmin))*n2 ); + long double s_r = s*x; + + // Boost velocity CM -> LAB + long double EnuinCM = sqrtl(s_r)/2.; + long double beta = (powl(Enuin,2)-powl(EnuinCM,2))/(powl(Enuin,2)+powl(EnuinCM,2)); + + // Final state primary lepton PDG code + int pdgl = interaction->FSPrimLeptonPdg(); + assert(pdgl!=0); + + if ( pdg::IsElectron(TMath::Abs(pdgl)) || pdg::IsMuon(TMath::Abs(pdgl)) || pdg::IsTau(TMath::Abs(pdgl)) ) { + + long double ElpoutCM = (s_r+mlout*mlout)/sqrtl(s_r)/2.; + long double EnuoutCM = (s_r-mlout*mlout)/sqrtl(s_r)/2.; + LongLorentzVector p4_lpout( 0., EnuoutCM*sinthCM, EnuoutCM*costhCM, ElpoutCM ); + LongLorentzVector p4_nuout( 0., -EnuoutCM*sinthCM, -EnuoutCM*costhCM, EnuoutCM ); + + p4_lpout.BoostZ(beta); + p4_nuout.BoostZ(beta); + + TLorentzVector p4lp_o( (double)p4_lpout.Px(), (double)p4_lpout.Py(), (double)p4_lpout.Pz(), (double)p4_lpout.E() ); + TLorentzVector p4nu_o( (double)p4_nuout.Px(), (double)p4_nuout.Py(), (double)p4_nuout.Pz(), (double)p4_nuout.E() ); + + // Randomize transverse components + RandomGen * rnd = RandomGen::Instance(); + double phi = 2* kPi * rnd->RndLep().Rndm(); + p4lp_o.RotateZ(phi); + p4nu_o.RotateZ(phi); + + //rotate from LAB=[0,0,Ev,Ev]->[px,py,pz,E] + p4lp_o.RotateUz(unit_nu); + p4nu_o.RotateUz(unit_nu); + + int pdgvout = 0; + if ( pdg::IsElectron(pdgl) ) pdgvout = kPdgAntiNuE; + else if ( pdg::IsPositron(pdgl) ) pdgvout = kPdgNuE; + else if ( pdg::IsMuon(pdgl) ) pdgvout = kPdgAntiNuMu; + else if ( pdg::IsAntiMuon(pdgl) ) pdgvout = kPdgNuMu; + else if ( pdg::IsTau(pdgl) ) pdgvout = kPdgAntiNuTau; + else if ( pdg::IsAntiTau(pdgl) ) pdgvout = kPdgNuTau; + + int pdgboson = pdg::IsNeutrino(probepdg) ? kPdgWP : kPdgWM; + + // Create a GHepParticle and add it to the event record + evrec->AddParticle( pdgboson, kIStDecayedState, 0, -1, 5, 6, *nu->P4()+*el->P4(), *(nu->X4()) ); //W [mothers: nuebar_in,e_in][daugthers: nulbar_out,lp_out] + evrec->AddParticle( pdgl, kIStStableFinalState, 4, -1, -1, -1, p4lp_o, *(nu->X4()) ); + evrec->AddParticle( pdgvout, kIStStableFinalState, 4, -1, -1, -1, p4nu_o, *(nu->X4()) ); + evrec->Summary()->KinePtr()->SetFSLeptonP4(p4lp_o); + + } + else { + + char p6frame[10]; + strcpy(p6frame, "CMS" ); + + char p6nu[10], p6tgt[10]; + if (pdg::IsAntiNeutrino(nu->Pdg())) { strcpy(p6nu, "nu_ebar"); strcpy(p6tgt, "e-"); } + else if (pdg::IsNeutrino (nu->Pdg())) { strcpy(p6nu, "nu_e"); strcpy(p6tgt, "e+"); } + + int def61 = fPythia->GetMSTP(61); + int def71 = fPythia->GetMSTP(71); + int def206 = fPythia->GetMDME(206,1); + int def207 = fPythia->GetMDME(207,1); + int def208 = fPythia->GetMDME(208,1); + fPythia->SetMSTP(61,0); // (Default=2) master switch for initial-state QCD and QED radiation. + fPythia->SetMSTP(71,0); // (Default=2) master switch for initial-state QCD and QED radiation. + fPythia->SetMDME(206,1,0); //swicht off W decay leptonic modes + fPythia->SetMDME(207,1,0); + fPythia->SetMDME(208,1,0); + + fPythia->Pyinit(p6frame, p6nu, p6tgt, sqrtl(s_r)); + fPythia->Pyevnt(); + + fPythia->SetMSTP(61,def61); + fPythia->SetMSTP(71,def71); + fPythia->SetMDME(206,1,def206); + fPythia->SetMDME(207,1,def207); + fPythia->SetMDME(208,1,def208); + + // get LUJETS record + fPythia->GetPrimaries(); + TClonesArray * pythia_particles = (TClonesArray *) fPythia->ImportParticles("All"); + int np = pythia_particles->GetEntries(); + assert(np>0); + + TMCParticle * particle = 0; + TIter piter(pythia_particles); + while( (particle = (TMCParticle *) piter.Next()) ) { + + int pdgc = particle->GetKF(); + int ks = particle->GetKS(); + + if ( ks==21 ) { continue; } //we dont want to save first particles from pythia (init states) + + LongLorentzVector p4longo(particle->GetPx(), particle->GetPy(), particle->GetPz(), particle->GetEnergy()); + p4longo.BoostZ(beta); + + TLorentzVector p4o( (double)p4longo.Px(), (double)p4longo.Py(), (double)p4longo.Pz(), (double)p4longo.E() ); + p4o.RotateUz(unit_nu); + + TParticlePDG * part = PDGLibrary::Instance()->Find(pdgc); + if ( (ks==1 || ks==4) && p4o.E() < part->Mass() ) { + LOG("PhotonRESGenerator", pWARN) << "Putting at rest one stable particle generated by PYTHIA because E < m"; + LOG("PhotonRESGenerator", pWARN) << "PDG = " << pdgc << " // State = " << ks; + LOG("PhotonRESGenerator", pWARN) << "E = " << p4o.E() << " // |p| = " << TMath::Sqrt(p4o.P()); + LOG("PhotonRESGenerator", pWARN) << "p = [ " << p4o.Px() << " , " << p4o.Py() << " , " << p4o.Pz() << " ]"; + LOG("PhotonRESGenerator", pWARN) << "m = " << p4o.M() << " // mpdg = " << part->Mass(); + p4o.SetXYZT(0,0,0,part->Mass()); + } + + // copy final state particles to the event record + GHepStatus_t ist = (ks==1 || ks==4) ? kIStStableFinalState : kIStDISPreFragmHadronicState; + + // fix numbering scheme used for mother/daughter assignments + int firstmother = -1; + int lastmother = -1; + int firstchild = -1; + int lastchild = -1; + + if ( particle->GetParent() < 10 ) { + if ( TMath::Abs(pdgc)<7 ) { //outgoing quarks: mother will be the boson (saved in position 4) + firstmother = 4; + firstchild = particle->GetFirstChild() - 6; + lastchild = particle->GetLastChild() - 6; + } + else if ( TMath::Abs(pdgc)==24 ) { //produced W boson: mother will be the incoming neutrino + firstmother = 0; + firstchild = particle->GetFirstChild() - 6; + lastchild = particle->GetLastChild() - 6; + } + else if ( pdgc==22 ) { //radiative photons: mother will be the incoming electron + firstmother = 2; + } + } + else { //rest + firstmother = particle->GetParent() - 6; //shift to match boson position + firstchild = (particle->GetFirstChild()==0) ? particle->GetFirstChild() - 1 : particle->GetFirstChild() - 6; + lastchild = (particle->GetLastChild()==0) ? particle->GetLastChild() - 1 : particle->GetLastChild() - 6; + } + + double lightspeed = 299792458e3; //c in mm/s. Used for time in PYTHIA t[s]=t_pythia[mm]/c[mm/s] + double vx = nu->X4()->X() + particle->GetVx()*1e12; //pythia gives position in [mm] while genie uses [fm] + double vy = nu->X4()->Y() + particle->GetVy()*1e12; + double vz = nu->X4()->Z() + particle->GetVz()*1e12; + double vt = nu->X4()->T() + particle->GetTime()/lightspeed; + TLorentzVector pos( vx, vy, vz, vt ); + + evrec->AddParticle(pdgc, ist, firstmother, lastmother, firstchild, lastchild, p4o, pos ); + + } + + delete particle; + pythia_particles->Clear("C"); + + } + +} +//___________________________________________________________________________ +void PhotonRESGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonRESGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonRESGenerator::LoadConfig(void) +{ + + // PYTHIA parameters only valid for HEDIS + double wmin; GetParam( "Xsec-Wmin", wmin ) ; + int warnings; GetParam( "Warnings", warnings ) ; + int errors; GetParam( "Errors", errors ) ; + int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ; + fPythia->SetPARP(2, wmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes) + fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed + fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed + fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term. + + fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385) + fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation + fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation + + GetParam( "Q2Grid-Min", fQ2PDFmin ); + +} +//____________________________________________________________________________ diff --git a/src/Physics/HELepton/EventGen/PhotonRESGenerator.h b/src/Physics/HELepton/EventGen/PhotonRESGenerator.h new file mode 100644 index 000000000..da90ef083 --- /dev/null +++ b/src/Physics/HELepton/EventGen/PhotonRESGenerator.h @@ -0,0 +1,59 @@ +//____________________________________________________________________________ +/*! + +\class genie::PhotonRESGenerator + +\brief Generator for trident production. + +\author Alfonso Garcia + IFIC & Harvard University + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _PHOTON_RES_GENERATOR_H_ +#define _PHOTON_RES_GENERATOR_H_ + +#include + +#include "Framework/EventGen/EventRecordVisitorI.h" +#include "Physics/HELepton/XSection/Born.h" + +namespace genie { + +class PhotonRESGenerator : public EventRecordVisitorI { + +public : + PhotonRESGenerator(); + PhotonRESGenerator(string config); + ~PhotonRESGenerator(); + + // implement the EventRecordVisitorI interface + void Initialize (void) const; + void ProcessEventRecord (GHepRecord * evrec) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + + void Configure(string config); + +private: + + void LoadConfig(void); + + mutable TPythia6 * fPythia; ///< PYTHIA6 wrapper class + + double fQ2PDFmin; + + Born * born; + +}; + +} // genie namespace +#endif // _PHOTON_RES_GENERATOR_H_ \ No newline at end of file diff --git a/src/Physics/HELepton/XSection/Born.cxx b/src/Physics/HELepton/XSection/Born.cxx new file mode 100644 index 000000000..7dfc07122 --- /dev/null +++ b/src/Physics/HELepton/XSection/Born.cxx @@ -0,0 +1,226 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include "Physics/HELepton/XSection/Born.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/Conventions/Constants.h" + +#include + +#include + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +Born::Born() +{ + + fGw = PDGLibrary::Instance()->Find(kPdgWM)->Width(); + fGz = PDGLibrary::Instance()->Find(kPdgZ0)->Width(); + fmw2c = TComplex(kMw2,fGw*kMw); + fmz2c = TComplex(kMz2,fGz*kMz); + TComplex rat = fmw2c/fmz2c; + fsw2 = TComplex(1.-rat.Re(),-rat.Im()); + fcw2 = 1.-fsw2; + falpha = TMath::Sqrt(2.)*kGF/kPi * fmw2c * fsw2; + + fgae = -1./2. + 2.*fsw2; + fgbe = -1./2.; + fgav = 1./2.; + +} +//____________________________________________________________________________ +Born::~Born() +{ + +/* +Make sure the p3 is always the charged lepton: + +nu (p1) + lp (p2) -> lp (p3) + nu (p4) + + s-channel t-channel u-channel +1 \ / 3 1 --------- 3 1 ------\ / 3 + ------ | | X +2 / \ 4 2 --------- 4 2 ------/ \ 4 +*/ + +} +//____________________________________________________________________________ +double Born::PXSecCCR(double s, double t, double mlin, double mlout) +{ + + TComplex prop = falpha/fsw2/(s-fmw2c); + + return (t-mlout*mlout)*(t-mlin*mlin) * prop.Rho2(); + +} +//____________________________________________________________________________ +double Born::PXSecCCV(double s, double t, double mlin, double mlout) +{ + + TComplex prop = falpha/fsw2/(t-fmw2c); + + return (s-mlout*mlout)*(s-mlin*mlin) * prop.Rho2(); + +} +//____________________________________________________________________________ +double Born::PXSecCCRNC(double s, double t, double mlin, double mlout) +{ + + double u = GetU(mlin,mlout,s,t); + + TComplex a = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2; + TComplex b = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2/fsw2 + 1./(s-fmw2c)/fsw2; + return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() ); + +} +//____________________________________________________________________________ +double Born::PXSecCCVNC(double s, double t, double mlin, double mlout) +{ + + double u = GetU(mlin,mlout,s,t); + + TComplex a = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2/fsw2 + 1./(t-fmw2c)/fsw2; + TComplex b = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2; + return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() ); + +} +//____________________________________________________________________________ +double Born::PXSecNCVnu(double s, double t, double mlin, double mlout) +{ + + double u = GetU(mlin,mlout,s,t); + + TComplex a = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2/fsw2; + TComplex b = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2; + return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() ); + +} +//____________________________________________________________________________ +double Born::PXSecNCVnubar(double s, double t, double mlin, double mlout) +{ + + double u = GetU(mlin,mlout,s,t); + + TComplex a = fgav*(fgae-fgbe)/(u-fmz2c)/fcw2/fsw2; + TComplex b = fgav*(fgae+fgbe)/(u-fmz2c)/fcw2; + return falpha.Rho2() * ( (t-mlout*mlout)*(t-mlin*mlin)*b.Rho2()/fsw2.Rho2() + (s-mlout*mlout)*(s-mlin*mlin)*a.Rho2() ); + +} +//____________________________________________________________________________ +double Born::PXSecPhoton(double s, double t, double mlout2) +{ + + double u = kMw2 + mlout2 - s - t; + + double ME = 8*kPi2*falpha.Rho2()*s/kMw2/fsw2.Re()/TMath::Power(mlout2 - t,2)/TMath::Power(kMw2 - u,2)* + ( -2*(kMw2 - u)*TMath::Power(mlout2,3) - 2*TMath::Power(mlout2,2)*(-2*kMw2*u + u*(s + u) + TMath::Power(kMw2,2)) + + mlout2*(-(kMw2*u*(4*s + 5*u)) + (s + u)*TMath::Power(kMw2,2) + 3*TMath::Power(kMw2,3) + (s + u)*TMath::Power(u,2)) + + 2*kMw2*((3*s + u)*TMath::Power(kMw2,2) - TMath::Power(kMw2,3) + 4*u*TMath::Power(s,2) + 2*TMath::Power(s,3) + 3*s*TMath::Power(u,2) + + TMath::Power(u,3) - kMw2*TMath::Power(2*s + u,2)) ); + + return TMath::Max(0.,ME); + +} +//____________________________________________________________________________ +double Born::PXSecPhoton_T(double s12, double s13, double Q2, double ml2) +{ + double ME2 = 0.0; + ME2 = (4*falpha.Rho2()*kPi2*(TMath::Power(ml2,4)*s12*(2*TMath::Power(kMw2,2)*TMath::Power(s12,2) - 2*kMw2*Q2*s12*s13 + TMath::Power(Q2,2)*s13*(-s12 + s13)) + + TMath::Power(ml2,3)*(-2*TMath::Power(kMw2,3)*TMath::Power(s12,3) + TMath::Power(Q2,2)*(s12 - s13)*s13*(-(Q2*s12) + TMath::Power(s12,2) + Q2*s13 - 3*s12*s13) + + 2*TMath::Power(kMw2,2)*TMath::Power(s12,2)*(3*Q2*s12 - 2*TMath::Power(s12,2) + 2*Q2*s13 + s12*s13) + + kMw2*Q2*s12*(Q2*TMath::Power(s12,2) - Q2*TMath::Power(s13,2) - 2*s12*TMath::Power(s13,2))) + + TMath::Power(ml2,2)*(-6*TMath::Power(kMw2,4)*TMath::Power(s12,3) + 2*kMw2*Q2* + (-(Q2*s12*(s12 - 3*s13)) + TMath::Power(Q2,2)*(s12 - s13) + TMath::Power(s12,2)*(s12 - s13))*(s12 - s13)*s13 + + TMath::Power(Q2,2)*(s12 - s13)*TMath::Power(s13,2)*(-2*Q2*s12 + 2*TMath::Power(s12,2) + 2*Q2*s13 - 3*s12*s13) + + 2*TMath::Power(kMw2,3)*TMath::Power(s12,2)*(2*TMath::Power(s12,2) - 3*Q2*(2*s12 + s13)) + + TMath::Power(kMw2,2)*s12*(2*TMath::Power(s12,2)*TMath::Power(s12 - s13,2) + TMath::Power(Q2,2)*(s12 - s13)*s13 - + 2*Q2*s12*(TMath::Power(s12,2) - 8*s12*s13 - 2*TMath::Power(s13,2)))) - + 2*TMath::Power(kMw2,2)*TMath::Power(s12,2)*(2*TMath::Power(kMw2,4)*s12 + TMath::Power(Q2,2)*TMath::Power(s13,2)*(-s12 + s13) - + 2*TMath::Power(kMw2,3)*(2*TMath::Power(s12,2) - Q2*s13 + s12*s13) + + TMath::Power(kMw2,2)*(4*Q2*s12*(s12 - 2*s13) + TMath::Power(Q2,2)*(-s12 + s13) + 2*s12*TMath::Power(s12 + s13,2)) + + 2*kMw2*s13*(TMath::Power(Q2,2)*(s12 - s13) - s12*(TMath::Power(s12,2) + TMath::Power(s13,2)) + Q2*(-TMath::Power(s12,2) + 2*s12*s13 + TMath::Power(s13,2)))) + + ml2*(10*TMath::Power(kMw2,5)*TMath::Power(s12,3) - TMath::Power(Q2,2)*(Q2 - s12)*TMath::Power(s12 - s13,2)*TMath::Power(s13,3) + + 2*TMath::Power(kMw2,4)*TMath::Power(s12,2)*(3*Q2*s12 - 4*TMath::Power(s12,2) + 4*Q2*s13 - 3*s12*s13) + + kMw2*Q2*(s12 - s13)*TMath::Power(s13,2)*(2*TMath::Power(Q2,2)*(s12 - s13) + 2*TMath::Power(s12,2)*(s12 - s13) + Q2*s12*(-3*s12 + 5*s13)) + + TMath::Power(kMw2,3)*s12*(2*Q2*s12*(5*TMath::Power(s12,2) - 16*s12*s13 - TMath::Power(s13,2)) + + TMath::Power(Q2,2)*(-3*TMath::Power(s12,2) + 2*s12*s13 + TMath::Power(s13,2)) + 2*TMath::Power(s12,2)*(TMath::Power(s12,2) + 6*s12*s13 + TMath::Power(s13,2))) - + TMath::Power(kMw2,2)*s13*(TMath::Power(Q2,3)*TMath::Power(s12 - s13,2) - 2*TMath::Power(s12,3)*TMath::Power(s12 - s13,2) + + TMath::Power(Q2,2)*s12*(-5*TMath::Power(s12,2) + 8*s12*s13 - 3*TMath::Power(s13,2)) + + 2*Q2*TMath::Power(s12,2)*(3*TMath::Power(s12,2) - 9*s12*s13 + 2*TMath::Power(s13,2))))))/(TMath::Power(kMw2,3)*TMath::Power(s12,2)*TMath::Power(ml2 - kMw2 + s13,2)*TMath::Power(ml2 - kMw2 - s12 + s13,2)*fsw2.Re()); + return TMath::Max(0.,ME2); +} +//____________________________________________________________________________ +double Born::PXSecPhoton_L(double s12, double s13, double Q2, double ml2) +{ + double ME2 = 0.0; + ME2 = 2*falpha.Rho2()*Q2*TMath::Power(kMw2,-3)*kPi2*TMath::Power(s12,-2)*TMath::Power(ml2 - kMw2 + s13,-2)*TMath::Power(ml2 - kMw2 - s12 + s13,-2)* + ((s12 - s13)*TMath::Power(ml2,5)*TMath::Power(s12,2) - 2*s12*TMath::Power(ml2,4)*(2*kMw2*TMath::Power(s12,2) - (Q2 - s12)*(-3*s12*s13 + TMath::Power(s12,2) + 2*TMath::Power(s13,2))) + + TMath::Power(ml2,3)*(-2*(s12 - s13)*TMath::Power(kMw2,2)*TMath::Power(s12,2) + + 2*kMw2*s12*(Q2*(9*s12*s13 - 5*TMath::Power(s12,2) - 2*TMath::Power(s13,2)) + s12*(-11*s12*s13 + 3*TMath::Power(s12,2) + 4*TMath::Power(s13,2))) + + (s12 - s13)*(TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) + TMath::Power(s12,2)*(-6*s12*s13 + TMath::Power(s12,2) + 6*TMath::Power(s13,2)) - + 2*Q2*s12*(-5*s12*s13 + TMath::Power(s12,2) + 6*TMath::Power(s13,2)))) + + 2*TMath::Power(ml2,2)*(4*(2*s12 - s13)*TMath::Power(kMw2,3)*TMath::Power(s12,2) + + (Q2 - s12)*s13*(Q2*(s12 - 2*s13) + s12*(-s12 + s13))*(-3*s12*s13 + TMath::Power(s12,2) + 2*TMath::Power(s13,2)) + + s12*TMath::Power(kMw2,2)*(Q2*(-11*s12*s13 + 7*TMath::Power(s12,2) - 2*TMath::Power(s13,2)) + s12*(15*s12*s13 - 11*TMath::Power(s12,2) + 2*TMath::Power(s13,2))) - + kMw2*((s12 - s13)*TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) + TMath::Power(s12,2)*(-11*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 20*s12*TMath::Power(s13,2) - 8*TMath::Power(s13,3)) - + 2*Q2*s12*(-10*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 17*s12*TMath::Power(s13,2) - 6*TMath::Power(s13,3)))) + + 4*TMath::Power(kMw2,2)*TMath::Power(s12,2)*(-(s13*(Q2 - 7*s12 + 4*s13)*TMath::Power(kMw2,2)) + (s12 - 2*s13)*TMath::Power(kMw2,3) + kMw2*(2*Q2 - s12 - 2*s13)*TMath::Power(s13,2) + + (-Q2 + s12)*TMath::Power(s13,3)) + ml2*(-15*(s12 - s13)*TMath::Power(kMw2,4)*TMath::Power(s12,2) + + 2*s12*TMath::Power(kMw2,3)*(Q2*(7*s12*s13 - 3*TMath::Power(s12,2) + 2*TMath::Power(s13,2)) + s12*(-21*s12*s13 + 9*TMath::Power(s12,2) + 4*TMath::Power(s13,2))) + + TMath::Power(kMw2,2)*((s12 - s13)*TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) + + TMath::Power(s12,2)*(-31*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 60*s12*TMath::Power(s13,2) - 14*TMath::Power(s13,3)) - + 2*Q2*s12*(-14*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 27*s12*TMath::Power(s13,2) - 6*TMath::Power(s13,3))) - + 2*kMw2*s13*((s12 - s13)*TMath::Power(Q2,2)*TMath::Power(s12 - 2*s13,2) + + TMath::Power(s12,2)*(-8*s13*TMath::Power(s12,2) + TMath::Power(s12,3) + 11*s12*TMath::Power(s13,2) - 4*TMath::Power(s13,3)) + + Q2*s12*(15*s13*TMath::Power(s12,2) - 2*TMath::Power(s12,3) - 25*s12*TMath::Power(s13,2) + 10*TMath::Power(s13,3))) + + (s12 - s13)*TMath::Power(s13,2)*TMath::Power(Q2*(s12 - 2*s13) + s12*(-s12 + s13),2)))/fsw2.Re(); + return TMath::Max(0.,ME2); +} +//____________________________________________________________________________ +double Born::Lambda(double a, double b, double c) +{ + return a*a + b*b + c*c - 2*a*b - 2*a*c - 2*b*c; +} +//____________________________________________________________________________ +double Born::GetS(double mlin, double Enuin) +{ + return 2. * mlin * Enuin + mlin*mlin; +} +//____________________________________________________________________________ +double Born::GetT(double mlin, double mlout, double s, double costhCM) +{ + //http://edu.itp.phys.ethz.ch/hs10/ppp1/PPP1_2.pdf [Sec 2.2.1] + double sum = mlin*mlin+mlout*mlout; + return ( (TMath::Sqrt(Lambda(s,0.,mlin*mlin)*Lambda(s,mlout*mlout,0.))*costhCM+mlin*mlin*mlout*mlout)/s + sum - s ) /2.; +} +//____________________________________________________________________________ +double Born::GetU(double mlin, double mlout, double s, double t) +{ + return mlin*mlin+mlout*mlout-s-t; +} +//____________________________________________________________________________ +bool Born::IsInPhaseSpace(double mlin, double mlout, double Enuin, double Enuout) +{ + + //https://arxiv.org/pdf/2007.14426.pdf [section 2.2] + double frac = Enuout/Enuin; + if ( frac < mlin/(mlin+2.*Enuin)+(mlout*mlout-mlin*mlin)/2./Enuin/(mlin+2.*Enuin) ) return false; + else if ( frac > 1.-(mlout*mlout-mlin*mlin)/2./Enuin/mlin ) return false; + + return true; + +} + + + + diff --git a/src/Physics/HELepton/XSection/Born.h b/src/Physics/HELepton/XSection/Born.h new file mode 100644 index 000000000..0f320bf10 --- /dev/null +++ b/src/Physics/HELepton/XSection/Born.h @@ -0,0 +1,66 @@ +//____________________________________________________________________________ +/*! + +\class genie::Born + +\brief Born level nu-electron cross section. + +\author Alfonso Garcia + IFIC & Harvard University + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _BORN_PXSEC_H_ +#define _BORN_PXSEC_H_ + +#include + +namespace genie { + +class Born { + +public: + Born (); + virtual ~Born (); + + double GetReAlpha (void) { return falpha.Re(); } + double PXSecCCR (double s, double t, double mlin, double mlout); + double PXSecCCV (double s, double t, double mlin, double mlout); + double PXSecCCRNC (double s, double t, double mlin, double mlout); + double PXSecCCVNC (double s, double t, double mlin, double mlout); + double PXSecNCVnu (double s, double t, double mlin, double mlout); + double PXSecNCVnubar (double s, double t, double mlin, double mlout); + double PXSecPhoton (double s, double t, double mlout2); + double PXSecPhoton_T (double s12, double s13, double Q2, double ml2); + double PXSecPhoton_L (double s12, double s13, double Q2, double ml2); + double GetS (double mlin, double Enuin); + double GetT (double mlin, double mlout, double s, double costhCM); + double GetU (double mlin, double mlout, double s, double t); + bool IsInPhaseSpace (double mlin, double mlout, double Enuin, double Enuout); + double Lambda (double a, double b, double c); + +private: + + double fGw; + double fGz; + + TComplex falpha; + TComplex fsw2; + TComplex fcw2; + TComplex fmw2c; + TComplex fmz2c; + TComplex fgae; + TComplex fgbe; + TComplex fgav; + +}; + +} // genie namespace + +#endif // _BORN_H_ diff --git a/src/Physics/HELepton/XSection/GLRESPXSec.cxx b/src/Physics/HELepton/XSection/GLRESPXSec.cxx new file mode 100644 index 000000000..e672670ca --- /dev/null +++ b/src/Physics/HELepton/XSection/GLRESPXSec.cxx @@ -0,0 +1,160 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include + +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/RefFrame.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/Utils/KineUtils.h" +#include "Physics/HELepton/XSection/GLRESPXSec.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +GLRESPXSec::GLRESPXSec() : +XSecAlgorithmI("genie::GLRESPXSec") +{ + born = new Born(); + +} +//____________________________________________________________________________ +GLRESPXSec::GLRESPXSec(string config) : +XSecAlgorithmI("genie::GLRESPXSec", config) +{ + +} +//____________________________________________________________________________ +GLRESPXSec::~GLRESPXSec() +{ + +} +//____________________________________________________________________________ +double GLRESPXSec::XSec( + const Interaction * interaction, KinePhaseSpace_t kps) const +{ + + if(! this -> ValidProcess (interaction) ) return 0.; + + const InitialState & init_state = interaction -> InitState(); + const Kinematics & kinematics = interaction -> Kine(); + const XclsTag & xclstag = interaction -> ExclTag(); + + int loutpdg = xclstag.FinalLeptonPdg(); + + double mlout = interaction->FSPrimLepton()->Mass(); //mass of charged lepton + double mlin = kElectronMass; //mass of incoming charged lepton + + double Enuin = init_state.ProbeE(kRfLab); + double s = born->GetS(mlin,Enuin); + + double n1 = kinematics.GetKV(kKVn1); + double n2 = kinematics.GetKV(kKVn2); + double t = born->GetT(mlin,mlout,s,n1); + if (t>0) return 0.; + + //nlo correction + double zeta = born->GetReAlpha()/kPi*(2.*TMath::Log(TMath::Sqrt(-t)/kElectronMass)-1.); + double omx = TMath::Power(n2, 1./zeta ); + double pdf_soft = TMath::Exp(zeta*(3./4.-TMath::EulerGamma()))/TMath::Gamma(1.+zeta) + omx*(omx-2.)/2./n2; + if ( omx<0. || omx>1. ) return 0.; + double s_r = s*(1. - omx); + double t_r = t*(1. - omx); + + //http://users.jyu.fi/~tulappi/fysh300sl11/l2.pdf [slide 22] + //remember we always define nuout as p4 + double Enuout = (mlin*mlin-t_r)/2./mlin; + if ( !born->IsInPhaseSpace(mlin,mlout,Enuin,Enuout) ) return 0.; + + double xsec = kPi/4./(s-mlin*mlin) * pdf_soft ; + + if ( pdg::IsPion(loutpdg) ) { + if ( TMath::Sqrt(s_r)PXSecCCRNC(s_r,t_r,mlin,mlout); + else ME = born->PXSecCCR (s_r,t_r,mlin,mlout); + xsec *= TMath::Max(0.,ME); + + //----- If requested return the free electron xsec even for nuclear target + if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec; + + //----- Scale for the number of scattering centers at the target + int Ne = init_state.Tgt().Z(); // num of scattering centers + xsec *= Ne; + + if(kps!=kPSn1n2fE) { + LOG("GLRESPXSec", pWARN) + << "Doesn't support transformation from " + << KinePhaseSpace::AsString(kPSn1n2fE) << " to " + << KinePhaseSpace::AsString(kps); + xsec = 0; + } + + LOG("GLRESPXSec", pINFO) << "dxsec/dn1dn2 (E= " << Enuin << ", n1= " << n1 << ", n2=" << n2 << ") = " << xsec; + + return xsec; + +} +//____________________________________________________________________________ +double GLRESPXSec::Integral(const Interaction * interaction) const +{ + double xsec = fXSecIntegrator->Integrate(this,interaction); + + return xsec; +} +//____________________________________________________________________________ +bool GLRESPXSec::ValidProcess(const Interaction* interaction) const +{ + if(interaction->TestBit(kISkipProcessChk)) return true; + + const ProcessInfo & proc_info = interaction->ProcInfo(); + if(!proc_info.IsGlashowResonance()) return false; + if(!proc_info.IsWeakCC()) return false; + + const InitialState & init_state = interaction -> InitState(); + if(!pdg::IsAntiNuE(init_state.ProbePdg())) return false; + + if(init_state.Tgt().HitNucIsSet()) return false; + + return true; +} +//____________________________________________________________________________ +void GLRESPXSec::Configure(const Registry & config) +{ + + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void GLRESPXSec::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void GLRESPXSec::LoadConfig(void) +{ + + //-- load the differential cross section integrator + fXSecIntegrator = dynamic_cast (this->SubAlg("XSec-Integrator")); + assert(fXSecIntegrator); + + GetParam( "Xsec-Wmin", fWmin ) ; + + +} \ No newline at end of file diff --git a/src/Physics/GlashowResonance/XSection/GLRESPXSec.h b/src/Physics/HELepton/XSection/GLRESPXSec.h similarity index 75% rename from src/Physics/GlashowResonance/XSection/GLRESPXSec.h rename to src/Physics/HELepton/XSection/GLRESPXSec.h index 6042f6bb7..67ac5c2f7 100644 --- a/src/Physics/GlashowResonance/XSection/GLRESPXSec.h +++ b/src/Physics/HELepton/XSection/GLRESPXSec.h @@ -3,18 +3,18 @@ \class genie::GLRESPXSec -\brief Nuebar cross section at the Glashow resonance (nuebar + e- -> W-). - Is a concrete implementation of the XSecAlgorithmI interface. +\brief Differential cross section for glashow resonance -\ref T.K.Gaisser, F.Halzen and T.Stanev, Physics Reports 258:173 (1995) +\author Alfonso Garcia + IFIC & Harvard University -\author Costas Andreopoulos - University of Liverpool & STFC Rutherford Appleton Laboratory +\ref Phys. Rev. D 100, 091301 (2019) -\created May 04, 2005 +\created Dec 8, 2021 \cpright Copyright (c) 2003-2020, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE */ //____________________________________________________________________________ @@ -22,6 +22,7 @@ #define _GLASHOW_RESONANCE_PXSEC_H_ #include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/HELepton/XSection/Born.h" namespace genie { @@ -51,6 +52,8 @@ class GLRESPXSec : public XSecAlgorithmI { double fWmin; ///< Minimum value of W + Born * born; + }; } // genie namespace diff --git a/src/Physics/GlashowResonance/XSection/GLRESXSec.cxx b/src/Physics/HELepton/XSection/HELeptonXSec.cxx similarity index 52% rename from src/Physics/GlashowResonance/XSection/GLRESXSec.cxx rename to src/Physics/HELepton/XSection/HELeptonXSec.cxx index 5f897aa3f..9cad6b56b 100644 --- a/src/Physics/GlashowResonance/XSection/GLRESXSec.cxx +++ b/src/Physics/HELepton/XSection/HELeptonXSec.cxx @@ -1,16 +1,10 @@ //____________________________________________________________________________ /* - Copyright (c) 2003-2019, The GENIE Collaboration + Copyright (c) 2003-2020, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org - or see $GENIE/LICENSE - - Author: Alfonso Garcia - NIKHEF (Amsterdam) - - For the class documentation see the corresponding header file. - - Important revisions after version 2.0.0 : + Alfonso Garcia + IFIC & Harvard University */ //____________________________________________________________________________ @@ -18,7 +12,7 @@ #include #include -#include "Physics/GlashowResonance/XSection/GLRESXSec.h" +#include "Physics/HELepton/XSection/HELeptonXSec.h" #include "Framework/Algorithm/AlgConfigPool.h" #include "Framework/Conventions/GBuild.h" #include "Framework/Conventions/Controls.h" @@ -40,39 +34,38 @@ using namespace genie::controls; using namespace genie::constants; //____________________________________________________________________________ -GLRESXSec::GLRESXSec() : -XSecIntegratorI("genie::GLRESXSec") +HELeptonXSec::HELeptonXSec() : +XSecIntegratorI("genie::HELeptonXSec") { } //____________________________________________________________________________ -GLRESXSec::GLRESXSec(string config) : -XSecIntegratorI("genie::GLRESXSec", config) +HELeptonXSec::HELeptonXSec(string config) : +XSecIntegratorI("genie::HELeptonXSec", config) { } //____________________________________________________________________________ -GLRESXSec::~GLRESXSec() +HELeptonXSec::~HELeptonXSec() { } //____________________________________________________________________________ -double GLRESXSec::Integrate( +double HELeptonXSec::Integrate( const XSecAlgorithmI * model, const Interaction * in) const { if(! model->ValidProcess(in) ) return 0.; const KPhaseSpace & kps = in->PhaseSpace(); if(!kps.IsAboveThreshold()) { - LOG("GLRESXSec", pDEBUG) << "*** Below energy threshold"; + LOG("HELeptonXSec", pDEBUG) << "*** Below energy threshold"; return 0; } + const ProcessInfo & proc_info = in->ProcInfo(); const InitialState & init_state = in->InitState(); double Ev = init_state.ProbeE(kRfLab); - - int NNucl = init_state.Tgt().Z(); // If the input interaction is off a nuclear target, then chek whether // the corresponding free nucleon cross section already exists at the @@ -80,64 +73,93 @@ double GLRESXSec::Integrate( // If yes, calculate the nuclear cross section based on that value. // XSecSplineList * xsl = XSecSplineList::Instance(); - if( !xsl->IsEmpty() ) { + if( !xsl->IsEmpty() && !proc_info.IsPhotonCoherent() ) { + Interaction * interaction = new Interaction(*in); Target * target = interaction->InitStatePtr()->TgtPtr(); - target->SetId(kPdgTgtFreeP); - if(xsl->SplineExists(model,interaction)) { + + int NNucl = 0; + bool skip = false; + if ( proc_info.IsGlashowResonance() ) { + NNucl = init_state.Tgt().Z(); + target->SetId(kPdgTgtFreeP); + } + else if ( proc_info.IsPhotonResonance() ) { + int nucpdgc = init_state.Tgt().HitNucPdg(); + if (pdg::IsProton(nucpdgc)) { + NNucl = init_state.Tgt().Z(); + target->SetId(kPdgTgtFreeP); + } + else { + NNucl = init_state.Tgt().N(); + target->SetId(kPdgTgtFreeN); + } + } + + if( xsl->SplineExists(model,interaction) ) { const Spline * spl = xsl->GetSpline(model, interaction); double xsec = spl->Evaluate(Ev); - LOG("GLRESXSec", pINFO) << "From XSecSplineList: XSec[ve-,free nucleon] (E = " << Ev << " GeV) = " << xsec; + LOG("HELeptonXSec", pINFO) << "From XSecSplineList: XSec["<ExclTag().FinalLeptonPdg()<< ",free nucleon] (E = " << Ev << " GeV) = " << xsec; if( !interaction->TestBit(kIAssumeFreeNucleon) ) { - xsec *= NNucl; - LOG("GLRESXSec", pINFO) << "XSec[ve-] (E = " << Ev << " GeV) = " << xsec; + xsec *= NNucl; + LOG("HELeptonXSec", pINFO) << "XSec["<ExclTag().FinalLeptonPdg()<< "] (E = " << Ev << " GeV) = " << xsec; } delete interaction; return xsec; } delete interaction; + } + + Interaction * interaction = new Interaction(*in); + interaction->SetBit(kISkipProcessChk); + ROOT::Math::IntegrationMultiDim::Type ig_type = utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); - Range1D_t yl = kps.Limits(kKVy); - - LOG("GLRESXSec", pDEBUG) << "y = (" << yl.min << ", " << yl.max << ")"; + double xsec = 0.; + ROOT::Math::IBaseFunctionMultiDim * func; + if ( proc_info.IsPhotonCoherent() ) { + double kine_min[3] = { -1., 0., 0. }; + double kine_max[3] = { 1., 1., 1. }; + func = new utils::gsl::d2Xsec_dn1dn2dn3_E(model, interaction); + ROOT::Math::IntegratorMultiDim ig(*func,ig_type,0,0,fGSLMaxEval); + xsec = ig.Integral(kine_min, kine_max); + } + else { + double kine_min[2] = { -1., 0. }; + double kine_max[2] = { 1., 1. }; + func = new utils::gsl::d2Xsec_dn1dn2_E(model, interaction); + ROOT::Math::IntegratorMultiDim ig(*func,ig_type,0,0,fGSLMaxEval); + xsec = ig.Integral(kine_min, kine_max); + } - Interaction * interaction = new Interaction(*in); - interaction->SetBit(kISkipProcessChk); - //interaction->SetBit(kISkipKinematicChk); + delete func; - ROOT::Math::IBaseFunctionOneDim * func = - new utils::gsl::dXSec_dy_E(model, interaction); - ROOT::Math::IntegrationOneDim::Type ig_type = - utils::gsl::Integration1DimTypeFromString(fGSLIntgType); - ROOT::Math::Integrator ig(*func,ig_type,1,fGSLRelTol,fGSLMaxEval); - double xsec = ig.Integral(yl.min, yl.max) * (1E-38 * units::cm2); + xsec *= (1E-38 * units::cm2); - LOG("GLRESXSec", pDEBUG) << "*** XSec[ve-] (E=" << interaction->InitState().ProbeE(kRfLab) << ") = " << xsec; + LOG("HELeptonXSec", pDEBUG) << "*** XSec["<ExclTag().FinalLeptonPdg()<< "] (E=" << interaction->InitState().ProbeE(kRfLab) << ") = " << xsec / (1E-38 * units::cm2); delete interaction; - delete func; return xsec; } //____________________________________________________________________________ -void GLRESXSec::Configure(const Registry & config) +void HELeptonXSec::Configure(const Registry & config) { Algorithm::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ -void GLRESXSec::Configure(string config) +void HELeptonXSec::Configure(string config) { Algorithm::Configure(config); this->LoadConfig(); } //____________________________________________________________________________ -void GLRESXSec::LoadConfig(void) +void HELeptonXSec::LoadConfig(void) { // Get GSL integration type & relative tolerance - GetParamDef("gsl-integration-type", fGSLIntgType, string("adaptive") ) ; - GetParamDef("gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ; + GetParamDef("gsl-integration-type", fGSLIntgType, string("vegas") ) ; + int max_eval ; GetParamDef( "gsl-max-eval", max_eval, 500000 ) ; fGSLMaxEval = (unsigned int) max_eval ; diff --git a/src/Physics/HELepton/XSection/HELeptonXSec.h b/src/Physics/HELepton/XSection/HELeptonXSec.h new file mode 100644 index 000000000..8d01a6cc2 --- /dev/null +++ b/src/Physics/HELepton/XSection/HELeptonXSec.h @@ -0,0 +1,51 @@ +//____________________________________________________________________________ +/*! + +\class genie::HELeptonXSec + +\brief Total cross section integrator for neutrino-electron + +\author Alfonso Garcia + IFIC & Harvard University + +\ref Phys. Rev. D 100, 091301 (2019) + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _HE_LEPTON_XSEC_H_ +#define _HE_LEPTON_XSEC_H_ + +#include "Physics/XSectionIntegration/XSecIntegratorI.h" + +namespace genie { + +class XSecAlgorithmI; +class Interaction; + +class HELeptonXSec : public XSecIntegratorI { + +public: + HELeptonXSec(); + HELeptonXSec(string config); + virtual ~HELeptonXSec(); + + //! XSecIntegratorI interface implementation + double Integrate(const XSecAlgorithmI * model, const Interaction * i) const; + + //! Overload the Algorithm::Configure() methods to load private data + //! members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + void LoadConfig (void); +}; + +} // genie namespace +#endif // _HE_LEPTON_XSEC_H_ diff --git a/src/Physics/HELepton/XSection/HENuElPXSec.cxx b/src/Physics/HELepton/XSection/HENuElPXSec.cxx new file mode 100644 index 000000000..ddbc404c0 --- /dev/null +++ b/src/Physics/HELepton/XSection/HENuElPXSec.cxx @@ -0,0 +1,162 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include + +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/RefFrame.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/Utils/KineUtils.h" +#include "Physics/HELepton/XSection/HENuElPXSec.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +HENuElPXSec::HENuElPXSec() : +XSecAlgorithmI("genie::HENuElPXSec") +{ + born = new Born(); + +} +//____________________________________________________________________________ +HENuElPXSec::HENuElPXSec(string config) : +XSecAlgorithmI("genie::HENuElPXSec", config) +{ + +} +//____________________________________________________________________________ +HENuElPXSec::~HENuElPXSec() +{ + +} +//____________________________________________________________________________ +double HENuElPXSec::XSec( + const Interaction * interaction, KinePhaseSpace_t kps) const +{ + + if(! this -> ValidProcess (interaction) ) return 0.; + + const ProcessInfo & proc_info = interaction->ProcInfo(); + const InitialState & init_state = interaction -> InitState(); + const Kinematics & kinematics = interaction -> Kine(); + + bool isCC = proc_info.IsWeakCC(); + + int probepdg = init_state.ProbePdg(); + + double mlout = interaction->FSPrimLepton()->Mass(); //mass of outgoing charged lepton + double mlin = kElectronMass; //mass of incoming charged lepton + + double Enuin = init_state.ProbeE(kRfLab); + double s = born->GetS(mlin,Enuin); + + double n1 = kinematics.GetKV(kKVn1); + double n2 = kinematics.GetKV(kKVn2); + double t = born->GetT( mlin, mlout, s, n1 ); + if (t>0) return 0.; + + //nlo correction + double zeta = born->GetReAlpha()/kPi*(2.*TMath::Log(TMath::Sqrt(-t)/kElectronMass)-1.); + double omx = TMath::Power(n2, 1./zeta ); + double pdf_soft = TMath::Exp(zeta*(3./4.-TMath::EulerGamma()))/TMath::Gamma(1.+zeta) + omx*(omx-2.)/2./n2; + if ( omx<0. || omx>1. ) return 0.; + double s_r = s*(1. - omx); + double t_r = t*(1. - omx); + + //http://users.jyu.fi/~tulappi/fysh300sl11/l2.pdf [slide 22] + //remember we always define nuout as p4 + double Enuout = (mlin*mlin-t_r)/2./mlin; + if ( !born->IsInPhaseSpace(mlin,mlout,Enuin,Enuout) ) return 0.; + + double xsec = kPi/4./(s-mlin*mlin) * pdf_soft ; + + double ME = 0; + if ( pdg::IsNuE(init_state.ProbePdg()) ) ME = born->PXSecCCVNC(s_r,t_r,mlin,mlout); + else { + if (isCC) ME = born->PXSecCCV(s_r,t_r,mlin,mlout); + else { + if (probepdg>0) ME = born->PXSecNCVnu (s_r,t_r,mlin,mlout); + else ME = born->PXSecNCVnubar(s_r,t_r,mlin,mlout); + } + } + xsec *= TMath::Max(0.,ME); + + //----- If requested return the free electron xsec even for nuclear target + if( interaction->TestBit(kIAssumeFreeElectron) ) return xsec; + + //----- Scale for the number of scattering centers at the target + int Ne = init_state.Tgt().Z(); // num of scattering centers + xsec *= Ne; + + if(kps!=kPSn1n2fE) { + LOG("HENuElPXSec", pWARN) + << "Doesn't support transformation from " + << KinePhaseSpace::AsString(kPSn1n2fE) << " to " + << KinePhaseSpace::AsString(kps); + xsec = 0; + } + + LOG("HENuElPXSec", pINFO) << "dxsec/dn1dn2 (E= " << Enuin << ", n1= " << n1 << ", n2=" << n2 << ") = " << xsec; + + return xsec; + +} +//____________________________________________________________________________ +double HENuElPXSec::Integral(const Interaction * interaction) const +{ + double xsec = fXSecIntegrator->Integrate(this,interaction); + + return xsec; +} +//____________________________________________________________________________ +bool HENuElPXSec::ValidProcess(const Interaction* interaction) const +{ + if(interaction->TestBit(kISkipProcessChk)) return true; + + const ProcessInfo & proc_info = interaction->ProcInfo(); + if(!proc_info.IsGlashowResonance()) return false; + + const InitialState & init_state = interaction -> InitState(); + if(pdg::IsAntiNuE(init_state.ProbePdg())) return false; + if(pdg::IsNuE(init_state.ProbePdg()) && !proc_info.IsWeakCC()) return false; + if(pdg::IsAntiNuMu(init_state.ProbePdg()) && proc_info.IsWeakCC()) return false; + if(pdg::IsAntiNuTau(init_state.ProbePdg()) && proc_info.IsWeakCC()) return false; + + if(init_state.Tgt().HitNucIsSet()) return false; + + return true; +} +//____________________________________________________________________________ +void HENuElPXSec::Configure(const Registry & config) +{ + + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void HENuElPXSec::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void HENuElPXSec::LoadConfig(void) +{ + + //-- load the differential cross section integrator + fXSecIntegrator = dynamic_cast (this->SubAlg("XSec-Integrator")); + assert(fXSecIntegrator); + +} \ No newline at end of file diff --git a/src/Physics/HELepton/XSection/HENuElPXSec.h b/src/Physics/HELepton/XSection/HENuElPXSec.h new file mode 100644 index 000000000..a6d8b621a --- /dev/null +++ b/src/Physics/HELepton/XSection/HENuElPXSec.h @@ -0,0 +1,59 @@ +//____________________________________________________________________________ +/*! + +\class genie::HENuElPXSec + +\brief Differential cross section for neutrino-electron + +\author Alfonso Garcia + IFIC & Harvard University + +\ref Phys. Rev. D 100, 091301 (2019) + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _HE_NUEL_PXSEC_H_ +#define _HE_NUEL_PXSEC_H_ + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/HELepton/XSection/Born.h" + +namespace genie { + +class XSecIntegratorI; + +class HENuElPXSec : public XSecAlgorithmI { + +public: + HENuElPXSec (); + HENuElPXSec (string config); + virtual ~HENuElPXSec (); + + // XSecAlgorithmI interface implementation + double XSec (const Interaction * i, KinePhaseSpace_t k) const; + double Integral (const Interaction * i) const; + bool ValidProcess (const Interaction * i) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + void LoadConfig (void); + + const XSecIntegratorI * fXSecIntegrator; ///< diff. xsec integrator + + Born * born; + +}; + +} // genie namespace + +#endif // _HE_NUELECTRON_PXSEC_H_ diff --git a/src/Physics/GlashowResonance/XSection/LinkDef.h b/src/Physics/HELepton/XSection/LinkDef.h similarity index 52% rename from src/Physics/GlashowResonance/XSection/LinkDef.h rename to src/Physics/HELepton/XSection/LinkDef.h index b9a25d208..e1272f52f 100644 --- a/src/Physics/GlashowResonance/XSection/LinkDef.h +++ b/src/Physics/HELepton/XSection/LinkDef.h @@ -6,7 +6,11 @@ #pragma link C++ namespace genie; -#pragma link C++ class genie::GLRESXSec; +#pragma link C++ class genie::HELeptonXSec; + #pragma link C++ class genie::GLRESPXSec; +#pragma link C++ class genie::HENuElPXSec; +#pragma link C++ class genie::PhotonRESPXSec; +#pragma link C++ class genie::PhotonCOHPXSec; #endif diff --git a/src/Physics/GlashowResonance/EventGen/Makefile b/src/Physics/HELepton/XSection/Makefile similarity index 86% rename from src/Physics/GlashowResonance/EventGen/Makefile rename to src/Physics/HELepton/XSection/Makefile index 14c97b77a..339d6f916 100644 --- a/src/Physics/GlashowResonance/EventGen/Makefile +++ b/src/Physics/HELepton/XSection/Makefile @@ -12,8 +12,8 @@ MAKEFILE = Makefile # include $(GENIE)/src/make/Make.include -PACKAGE = Physics/GlashowResonance/EventGen -PACKAGE_ABBREV = PhGlwResEG +PACKAGE = Physics/HELepton/XSection +PACKAGE_ABBREV = PhHELptnXS DICTIONARY = _ROOT_DICT_$(PACKAGE_ABBREV) LIBNAME = libG$(PACKAGE_ABBREV) EXTRA_EXT_LIBS = diff --git a/src/Physics/HELepton/XSection/PhotonCOHPXSec.cxx b/src/Physics/HELepton/XSection/PhotonCOHPXSec.cxx new file mode 100644 index 000000000..82da79628 --- /dev/null +++ b/src/Physics/HELepton/XSection/PhotonCOHPXSec.cxx @@ -0,0 +1,165 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include "Physics/HELepton/XSection/PhotonCOHPXSec.h" +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/RefFrame.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/Utils/KineUtils.h" + +using namespace genie; +using namespace genie::constants; + +const double a = 0.523 * (units::fermi); +const double r0 = 1.126 * (units::fermi); + +//____________________________________________________________________________ +PhotonCOHPXSec::PhotonCOHPXSec() : +XSecAlgorithmI("genie::PhotonCOHPXSec") +{ + born = new Born(); +} +//____________________________________________________________________________ +PhotonCOHPXSec::PhotonCOHPXSec(string config) : +XSecAlgorithmI("genie::PhotonCOHPXSec", config) +{ + +} +//____________________________________________________________________________ +PhotonCOHPXSec::~PhotonCOHPXSec() +{ + +} +//____________________________________________________________________________ +double PhotonCOHPXSec::XSec( + const Interaction * interaction, KinePhaseSpace_t kps) const +{ + + if(! this -> ValidProcess (interaction) ) return 0.; + + const InitialState & init_state = interaction -> InitState(); + const Kinematics & kinematics = interaction -> Kine(); + + int probepdg = init_state.ProbePdg(); + double E = init_state.ProbeE(kRfLab); + + double mlout = 0; + if (pdg::IsNuE (TMath::Abs(probepdg))) mlout = kElectronMass; + else if (pdg::IsNuMu (TMath::Abs(probepdg))) mlout = kMuonMass; + else if (pdg::IsNuTau(TMath::Abs(probepdg))) mlout = kTauMass; + double mlout2 = mlout*mlout; + + int A = init_state.Tgt().A(); + int Z = init_state.Tgt().Z(); + double Mtgt = Z*kProtonMass + (A-Z)*kNeutronMass; + + double n1 = kinematics.GetKV(kKVn1); + double n2 = kinematics.GetKV(kKVn2); + double n3 = kinematics.GetKV(kKVn3); + + double mL = mlout+kMw; + + double Delta = TMath::Power(2.*E*Mtgt-mL*mL,2)-4.*TMath::Power(Mtgt*mL,2); + if (Delta<0) return 0.; + + double s12_min = E/(2.*E+Mtgt)*(mL*mL+2.*E*Mtgt-TMath::Sqrt(Delta)); + double s12_max = E/(2.*E+Mtgt)*(mL*mL+2.*E*Mtgt+TMath::Sqrt(Delta)); + double s12 = TMath::Exp( TMath::Log(s12_min)+(TMath::Log(s12_max)-TMath::Log(s12_min))*n2); + + double Q2_min = TMath::Power(s12,2)*Mtgt/2./E/(2.*E*Mtgt-s12); + double Q2_max = s12 - mL*mL; + double Q2 = TMath::Exp( TMath::Log(Q2_min) + (TMath::Log(Q2_max)-TMath::Log(Q2_min))*n3 ); + + double s = s12 - Q2; + double s13 = s12/2.*((1.+kMw2/s-mlout2/s)-TMath::Sqrt(born->Lambda(1.,kMw2/s,mlout2/s))*n1); + + double ME_T = born->PXSecPhoton_T(s12,s13,Q2,mlout2) * (1.-s12/2./E/Mtgt-TMath::Power(s12/2./E,2)/Q2); + double ME_L = born->PXSecPhoton_L(s12,s13,Q2,mlout2) * TMath::Power(1.-s12/4./E/Mtgt,2); + + double ME = ME_T+ME_L; + double dps2 = 1/32./kPi2/s12 * TMath::Sqrt( born->Lambda(1.,kMw2/s,mlout2/s) ) * (TMath::Log(Q2_max)-TMath::Log(Q2_min)) * (TMath::Log(s12_max)-TMath::Log(s12_min)); + double dP = born->GetReAlpha()*TMath::Power(Z,2)*F2_Q( TMath::Sqrt(Q2), r0*TMath::Power(A,1./3.) ); + + double xsec = ME * dps2 * dP; + + if(kps!=kPSn1n2n3fE) { + LOG("PhotonCOHPXSec", pWARN) + << "Doesn't support transformation from " + << KinePhaseSpace::AsString(kPSn1n2n3fE) << " to " + << KinePhaseSpace::AsString(kps); + xsec = 0; + } + + // If requested return the free nucleon xsec even for input nuclear tgt + if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec; + + LOG("PhotonCOHPXSec", pINFO) << "dxsec/dn1dn2dn3 (E= " << E << ", n1= " << n1 << ", n2=" << n2 << ", n3=" << n3 << ") = " << xsec; + + return xsec; +} +//____________________________________________________________________________ +double PhotonCOHPXSec::F2_Q(double Q, double r0) const +{ + // Analytic Woods-Saxon, A.3 of https://arxiv.org/pdf/1807.10973.pdf + double FF = 0.0; + double coth = 1./TMath::TanH(kPi*Q*a); + FF = 3.*kPi*a/(TMath::Power(r0,2)+TMath::Power(kPi*a,2)) / (Q*r0* TMath::SinH(kPi*Q*a)); + FF *= (kPi*a*coth*TMath::Sin(Q*r0) - r0 * TMath::Cos(Q*r0)); + return TMath::Power(FF,2); +} +//____________________________________________________________________________ +double PhotonCOHPXSec::Integral(const Interaction * interaction) const +{ + double xsec = fXSecIntegrator->Integrate(this,interaction); + return xsec; +} +//____________________________________________________________________________ +bool PhotonCOHPXSec::ValidProcess(const Interaction* interaction) const +{ + if(interaction->TestBit(kISkipProcessChk)) return true; + + const ProcessInfo & proc_info = interaction->ProcInfo(); + if(!proc_info.IsPhotonCoherent()) return false; + + const InitialState & init_state = interaction -> InitState(); + if(!pdg::IsLepton(init_state.ProbePdg())) return false; + + if(init_state.Tgt().HitNucIsSet()) return false; + + return true; +} +//____________________________________________________________________________ + + +//____________________________________________________________________________ +void PhotonCOHPXSec::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonCOHPXSec::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonCOHPXSec::LoadConfig(void) +{ + + //-- load the differential cross section integrator + fXSecIntegrator = dynamic_cast (this->SubAlg("XSec-Integrator")); + assert(fXSecIntegrator); + +} diff --git a/src/Physics/HELepton/XSection/PhotonCOHPXSec.h b/src/Physics/HELepton/XSection/PhotonCOHPXSec.h new file mode 100644 index 000000000..501b00c91 --- /dev/null +++ b/src/Physics/HELepton/XSection/PhotonCOHPXSec.h @@ -0,0 +1,61 @@ +//____________________________________________________________________________ +/*! + +\class genie::PhotonCOHPXSec + +\brief Differential cross section for W boson production + +\author Alfonso Garcia + IFIC & Harvard University + +\ref Phys. Rev. D 100, 091301 (2019) + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _PHOTON_COH_PXSEC_H_ +#define _PHOTON_COH_PXSEC_H_ + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/HELepton/XSection/Born.h" + +namespace genie { + +class XSecIntegratorI; + +class PhotonCOHPXSec : public XSecAlgorithmI { + +public: + PhotonCOHPXSec (); + PhotonCOHPXSec (string config); + virtual ~PhotonCOHPXSec (); + + // XSecAlgorithmI interface implementation + double XSec (const Interaction * i, KinePhaseSpace_t k) const; + double Integral (const Interaction * i) const; + bool ValidProcess (const Interaction * i) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + void LoadConfig (void); + double F2_Q (double Q, double r0) const; //EM Nuclear Form-factor + + + const XSecIntegratorI * fXSecIntegrator; ///< diff. xsec integrator + + Born * born; + +}; + +} // genie namespace + +#endif // _PHOTON_COH_PXSEC_H_ \ No newline at end of file diff --git a/src/Physics/HELepton/XSection/PhotonRESPXSec.cxx b/src/Physics/HELepton/XSection/PhotonRESPXSec.cxx new file mode 100644 index 000000000..1f0124bbd --- /dev/null +++ b/src/Physics/HELepton/XSection/PhotonRESPXSec.cxx @@ -0,0 +1,162 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include "Physics/HELepton/XSection/PhotonStrucFunc.h" +#include "Physics/HELepton/XSection/PhotonRESPXSec.h" +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/RefFrame.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/Utils/KineUtils.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +PhotonRESPXSec::PhotonRESPXSec() : +XSecAlgorithmI("genie::PhotonRESPXSec") +{ + born = new Born(); +} +//____________________________________________________________________________ +PhotonRESPXSec::PhotonRESPXSec(string config) : +XSecAlgorithmI("genie::PhotonRESPXSec", config) +{ + +} +//____________________________________________________________________________ +PhotonRESPXSec::~PhotonRESPXSec() +{ + +} +//____________________________________________________________________________ +double PhotonRESPXSec::XSec( + const Interaction * interaction, KinePhaseSpace_t kps) const +{ + + if(! this -> ValidProcess (interaction) ) return 0.; + + // Load SF tables + PhotonStrucFunc * sf_tbl = PhotonStrucFunc::Instance(); + + const InitialState & init_state = interaction -> InitState(); + const Kinematics & kinematics = interaction -> Kine(); + const XclsTag & xclstag = interaction -> ExclTag(); + + int probepdg = init_state.ProbePdg(); + int loutpdg = xclstag.FinalLeptonPdg(); + int tgtpdg = init_state.Tgt().HitNucPdg(); + + double mlout = interaction->FSPrimLepton()->Mass(); //mass of charged lepton + + double Mnuc = init_state.Tgt().HitNucMass(); + + double Enuin = init_state.ProbeE(kRfLab); + double s = born->GetS(Mnuc,Enuin); + + double n1 = kinematics.GetKV(kKVn1); + double n2 = kinematics.GetKV(kKVn2); + + double xmin = fQ2PDFmin/2./Enuin/Mnuc; + double x = TMath::Exp( TMath::Log(xmin) + (TMath::Log(1.0)-TMath::Log(xmin))*n2 ); + + if (xGetT(0.,mlout,s_r,n1); + + double xsec = kPi/4./(s_r-Mnuc*Mnuc) * sf_tbl->EvalSF(tgtpdg,probepdg,x) * (TMath::Log(1.0)-TMath::Log(xmin)) ; + + if ( pdg::IsPion(loutpdg) ) { + if ( TMath::Sqrt(s_r)PXSecCCRNC(s_r,t_r,0.,mlout); + else ME = born->PXSecCCR (s_r,t_r,0.,mlout); + xsec *= TMath::Max(0.,ME); + + if(kps!=kPSn1n2fE) { + LOG("PhotonRESPXSec", pWARN) + << "Doesn't support transformation from " + << KinePhaseSpace::AsString(kPSn1n2fE) << " to " + << KinePhaseSpace::AsString(kps); + xsec = 0; + } + + // If requested return the free nucleon xsec even for input nuclear tgt + if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec; + + // Compute nuclear cross section (simple scaling here, corrections must have been included in the structure functions) + int NNucl = (pdg::IsProton(tgtpdg)) ? init_state.Tgt().Z() : init_state.Tgt().N(); + xsec *= NNucl; + + LOG("PhotonRESPXSec", pINFO) << "dxsec/dn1dn2 (E= " << Enuin << ", n1= " << n1 << ", n2=" << n2 << ") = " << xsec; + + return xsec; + +} +//____________________________________________________________________________ +double PhotonRESPXSec::Integral(const Interaction * interaction) const +{ + double xsec = fXSecIntegrator->Integrate(this,interaction); + return xsec; +} +//____________________________________________________________________________ +bool PhotonRESPXSec::ValidProcess(const Interaction* interaction) const +{ + if(interaction->TestBit(kISkipProcessChk)) return true; + + const ProcessInfo & proc_info = interaction->ProcInfo(); + if(!proc_info.IsPhotonResonance()) return false; + + const InitialState & init_state = interaction -> InitState(); + if(!pdg::IsLepton(init_state.ProbePdg())) return false; + + if(! init_state.Tgt().HitNucIsSet()) return false; + + int hitnuc_pdg = init_state.Tgt().HitNucPdg(); + if(!pdg::IsNeutronOrProton(hitnuc_pdg)) return false; + + return true; +} +//____________________________________________________________________________ + + +//____________________________________________________________________________ +void PhotonRESPXSec::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonRESPXSec::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void PhotonRESPXSec::LoadConfig(void) +{ + + //-- load the differential cross section integrator + fXSecIntegrator = dynamic_cast (this->SubAlg("XSec-Integrator")); + assert(fXSecIntegrator); + + GetParam( "Xsec-Wmin", fWmin ) ; + + GetParam("Q2Grid-Min", fQ2PDFmin ); + GetParam("XGrid-Min", fxPDFmin ); + +} \ No newline at end of file diff --git a/src/Physics/HELepton/XSection/PhotonRESPXSec.h b/src/Physics/HELepton/XSection/PhotonRESPXSec.h new file mode 100644 index 000000000..8929a9336 --- /dev/null +++ b/src/Physics/HELepton/XSection/PhotonRESPXSec.h @@ -0,0 +1,65 @@ +//____________________________________________________________________________ +/*! + +\class genie::PhotonRESPXSec + +\brief Differential cross section for trident production + +\author Alfonso Garcia + IFIC & Harvard University + +\ref Phys. Rev. D 100, 091301 (2019) + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _PHOTON_RES_PXSEC_H_ +#define _PHOTON_RES_PXSEC_H_ + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Physics/HELepton/XSection/Born.h" + +namespace genie { + +class XSecIntegratorI; + +class PhotonRESPXSec : public XSecAlgorithmI { + +public: + PhotonRESPXSec (); + PhotonRESPXSec (string config); + virtual ~PhotonRESPXSec (); + + // XSecAlgorithmI interface implementation + double XSec (const Interaction * i, KinePhaseSpace_t k) const; + double Integral (const Interaction * i) const; + bool ValidProcess (const Interaction * i) const; + + // overload the Algorithm::Configure() methods to load private data + // members from configuration options + void Configure(const Registry & config); + void Configure(string config); + +private: + void LoadConfig (void); + + const XSecIntegratorI * fXSecIntegrator; ///< diff. xsec integrator + + double fWmin; ///< Minimum value of W + + double fQ2PDFmin; + double fxPDFmin; + + Born * born; + + +}; + +} // genie namespace + +#endif // _PHOTON_RES_PXSEC_H_ diff --git a/src/Physics/HELepton/XSection/PhotonStrucFunc.cxx b/src/Physics/HELepton/XSection/PhotonStrucFunc.cxx new file mode 100644 index 000000000..b2596e1e6 --- /dev/null +++ b/src/Physics/HELepton/XSection/PhotonStrucFunc.cxx @@ -0,0 +1,58 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Alfonso Garcia + IFIC & Harvard University +*/ +//____________________________________________________________________________ + +#include "Physics/HELepton/XSection/PhotonStrucFunc.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGCodes.h" + +#include "TSystem.h" + +using namespace genie; + +//_________________________________________________________________________ +PhotonStrucFunc * PhotonStrucFunc::fgInstance = 0; +//_________________________________________________________________________ +PhotonStrucFunc::PhotonStrucFunc() +{ + + int nucs[2] = { kPdgProton, kPdgNeutron }; + int pdgs[6] = { kPdgNuE, kPdgAntiNuE, kPdgNuMu, kPdgAntiNuMu, kPdgNuTau, kPdgAntiNuTau }; + + for (int k=0; k<2; k++) { + for(int j=0; j<6; j++) { + string SFname = string(gSystem->Getenv("GENIE")) + "/data/evgen/photon-sf/PhotonSF_hitnuc"+std::to_string(nucs[k])+"_hitlep"+std::to_string(pdgs[j])+".dat"; + if ( gSystem->AccessPathName( SFname.c_str()) ) { + LOG("PhotonStrucFunc", pWARN) << "File doesnt exist. SF table must be compute with gmkglressf."; + assert(0); + } + fSFTables[nucs[k]].Table[pdgs[j]] = new genie::Spline(); + fSFTables[nucs[k]].Table[pdgs[j]]->LoadFromAsciiFile(SFname); + } + } + + fgInstance = 0; + +} +//_________________________________________________________________________ +PhotonStrucFunc::~PhotonStrucFunc() +{ + +} +//_________________________________________________________________________ +PhotonStrucFunc * PhotonStrucFunc::Instance() +{ + if(fgInstance == 0) { + LOG("PhotonStrucFunc", pINFO) << "Late initialization"; + static PhotonStrucFunc::Cleaner cleaner; + cleaner.DummyMethodAndSilentCompiler(); + fgInstance = new PhotonStrucFunc(); + } + return fgInstance; +} \ No newline at end of file diff --git a/src/Physics/HELepton/XSection/PhotonStrucFunc.h b/src/Physics/HELepton/XSection/PhotonStrucFunc.h new file mode 100644 index 000000000..2d67df849 --- /dev/null +++ b/src/Physics/HELepton/XSection/PhotonStrucFunc.h @@ -0,0 +1,87 @@ +//____________________________________________________________________________ +/*! + +\class genie::PhotonStrucFunc + +\brief Structure function using photon PDFs of nucleons + +\author Alfonso Garcia + IFIC & Harvard University + +\created Dec 8, 2021 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + or see $GENIE/LICENSE +*/ +//____________________________________________________________________________ + +#ifndef _PHOTON_STRUC_FUNC_H_ +#define _PHOTON_STRUC_FUNC_H_ + +#include "Framework/Numerical/Spline.h" + +#include + +using std::map; + +namespace genie { + + struct SF_x { + double Fp1,Fm1; + double Fp2,Fm2; + double Fp3,Fm3; + }; + + class PhotonStrucFunc + { + public: + + // ................................................................ + // GLRES form factor type + // + + class PhotonStrucFuncTable + { + public: + PhotonStrucFuncTable() { } + ~PhotonStrucFuncTable() { /* note: should delete the grids! */ } + map< int, genie::Spline * > Table; + }; + + // ................................................................ + + static PhotonStrucFunc * Instance(void); + + double EvalSF ( int hitnuc, int hitlep, double x ) { return fSFTables[hitnuc].Table[hitlep]->Evaluate(x); } + + private: + + // Ctors & dtor + PhotonStrucFunc(); + PhotonStrucFunc(const PhotonStrucFunc &); + ~PhotonStrucFunc(); + + // Self + static PhotonStrucFunc * fgInstance; + + // These map holds all SF tables (interaction channel is the key) + map fSFTables; + + + // singleton cleaner + struct Cleaner { + void DummyMethodAndSilentCompiler(){} + ~Cleaner(){ + if (PhotonStrucFunc::fgInstance !=0){ + delete PhotonStrucFunc::fgInstance; + PhotonStrucFunc::fgInstance = 0; + } + } + }; + friend struct Cleaner; + }; + +} // genie namespace + +#endif // _PHOTON_STRUC_FUNC_H_ \ No newline at end of file diff --git a/src/Physics/Hadronization/LeptoHadronization.cxx b/src/Physics/Hadronization/LeptoHadronization.cxx index d6a5ea1af..35a4ac73e 100644 --- a/src/Physics/Hadronization/LeptoHadronization.cxx +++ b/src/Physics/Hadronization/LeptoHadronization.cxx @@ -86,7 +86,7 @@ void LeptoHadronization::ProcessEventRecord(GHepRecord * #ifdef __GENIE_PYTHIA6_ENABLED__ if(!this->Hadronize(event)) { - LOG("PythiaHad", pWARN) << "Hadronization failed!"; + LOG("LeptoHad", pWARN) << "Hadronization failed!"; event->EventFlags()->SetBitNumber(kHadroSysGenErr, true); genie::exceptions::EVGThreadException exception; exception.SetReason("Could not simulate the hadronic system"); @@ -434,7 +434,7 @@ bool LeptoHadronization::Hadronize(GHepRecord * } LongLorentzVector p4long( p->GetPx(), p->GetPy(), p->GetPz(), p->GetEnergy() ); - p4long.Boost(beta); + p4long.BoostZ(beta); p4long.Rotate(p4Hadlong); // Translate from long double to double @@ -502,6 +502,7 @@ void LeptoHadronization::LoadConfig(void) // It is, with quark masses added, used to define the minimum allowable energy of a colour-singlet parton system. GetParam( "Energy-Singlet", fMinESinglet ) ; +#ifdef __GENIE_PYTHIA6_ENABLED__ // PYTHIA parameters only valid for HEDIS GetParam( "Xsec-Wmin", fWmin ) ; fPythia->SetPARP(2, fWmin); //(D = 10. GeV) lowest c.m. energy for the event as a whole that the program will accept to simulate. (bellow 2GeV pythia crashes) @@ -510,10 +511,13 @@ void LeptoHadronization::LoadConfig(void) int errors; GetParam( "Errors", errors ) ; int qrk_mass; GetParam( "QuarkMass", qrk_mass ) ; -#ifdef __GENIE_PYTHIA6_ENABLED__ fPythia->SetMSTU(26, warnings); // (Default=10) maximum number of warnings that are printed fPythia->SetMSTU(22, errors); // (Default=10) maximum number of errors that are printed fPythia->SetMSTJ(93, qrk_mass); // light (d, u, s, c, b) quark masses are taken from PARF(101) - PARF(105) rather than PMAS(1,1) - PMAS(5,1). Diquark masses are given as sum of quark masses, without spin splitting term. + + fPythia->SetPMAS(24,1,kMw); //mass of the W boson (pythia=80.450 // genie=80.385) + fPythia->SetPMAS(24,2,0.); //set to 0 the width of the W boson to avoid problems with energy conservation + fPythia->SetPMAS(6,2,0.); //set to 0 the width of the top to avoid problems with energy conservation #endif } diff --git a/src/Physics/XSectionIntegration/GSLXSecFunc.cxx b/src/Physics/XSectionIntegration/GSLXSecFunc.cxx index b68dc2923..0ed084198 100644 --- a/src/Physics/XSectionIntegration/GSLXSecFunc.cxx +++ b/src/Physics/XSectionIntegration/GSLXSecFunc.cxx @@ -1084,4 +1084,100 @@ ROOT::Math::IBaseFunctionMultiDim * genie::utils::gsl::dXSec_Log_Wrapper::Clone return new dXSec_Log_Wrapper(fFn,fIfLog,fMins,fMaxes); } +//_____________________________________________________________________________ +genie::utils::gsl::d2Xsec_dn1dn2_E::d2Xsec_dn1dn2_E( + const XSecAlgorithmI * m, const Interaction * i, double scale) : +ROOT::Math::IBaseFunctionMultiDim(), +fModel(m), +fInteraction(i), +fScale(scale) +{ + +} +//____________________________________________________________________________ +genie::utils::gsl::d2Xsec_dn1dn2_E::~d2Xsec_dn1dn2_E() +{ + +} +//____________________________________________________________________________ +unsigned int genie::utils::gsl::d2Xsec_dn1dn2_E::NDim(void) const +{ + return 2; +} //____________________________________________________________________________ +double genie::utils::gsl::d2Xsec_dn1dn2_E::DoEval(const double * xin) const +{ +// inputs: +// n1 +// n2 +// outputs: +// differential cross section (hbar=c=1 units) +// + + double n1 = xin[0]; + double n2 = xin[1]; + + Kinematics * kinematics = fInteraction->KinePtr(); + kinematics->SetKV(kKVn1, n1); + kinematics->SetKV(kKVn2, n2); + + double xsec = fModel->XSec(fInteraction, kPSn1n2fE); + return fScale*xsec/(1E-38 * units::cm2); +} +//____________________________________________________________________________ +ROOT::Math::IBaseFunctionMultiDim * + genie::utils::gsl::d2Xsec_dn1dn2_E::Clone() const +{ + return + new genie::utils::gsl::d2Xsec_dn1dn2_E(fModel,fInteraction); +} +//_____________________________________________________________________________ +genie::utils::gsl::d2Xsec_dn1dn2dn3_E::d2Xsec_dn1dn2dn3_E( + const XSecAlgorithmI * m, const Interaction * i, double scale) : +ROOT::Math::IBaseFunctionMultiDim(), +fModel(m), +fInteraction(i), +fScale(scale) +{ + +} +//____________________________________________________________________________ +genie::utils::gsl::d2Xsec_dn1dn2dn3_E::~d2Xsec_dn1dn2dn3_E() +{ + +} +//____________________________________________________________________________ +unsigned int genie::utils::gsl::d2Xsec_dn1dn2dn3_E::NDim(void) const +{ + return 3; +} +//____________________________________________________________________________ +double genie::utils::gsl::d2Xsec_dn1dn2dn3_E::DoEval(const double * xin) const +{ +// inputs: +// n1 +// n2 +// n3 +// outputs: +// differential cross section (hbar=c=1 units) +// + + double n1 = xin[0]; + double n2 = xin[1]; + double n3 = xin[2]; + + Kinematics * kinematics = fInteraction->KinePtr(); + kinematics->SetKV(kKVn1, n1); + kinematics->SetKV(kKVn2, n2); + kinematics->SetKV(kKVn3, n3); + + double xsec = fModel->XSec(fInteraction, kPSn1n2n3fE); + return fScale*xsec/(1E-38 * units::cm2); +} +//____________________________________________________________________________ +ROOT::Math::IBaseFunctionMultiDim * + genie::utils::gsl::d2Xsec_dn1dn2dn3_E::Clone() const +{ + return + new genie::utils::gsl::d2Xsec_dn1dn2dn3_E(fModel,fInteraction); +} diff --git a/src/Physics/XSectionIntegration/GSLXSecFunc.h b/src/Physics/XSectionIntegration/GSLXSecFunc.h index 5dc313a76..4910735ce 100644 --- a/src/Physics/XSectionIntegration/GSLXSecFunc.h +++ b/src/Physics/XSectionIntegration/GSLXSecFunc.h @@ -472,6 +472,47 @@ class dXSec_Log_Wrapper: public ROOT::Math::IBaseFunctionMultiDim double * fMaxes; }; +//..................................................................................... +// +// genie::utils::gsl::d2Xsec_dn1dn2_E +// A 2-D cross section function: d2xsec/dn1dn2 = f(n1,n2)|(fixed E) +// +class d2Xsec_dn1dn2_E: public ROOT::Math::IBaseFunctionMultiDim +{ + public: + d2Xsec_dn1dn2_E(const XSecAlgorithmI * m, const Interaction * i, double scale=1. ); + ~d2Xsec_dn1dn2_E(); + // ROOT::Math::IBaseFunctionMultiDim interface + unsigned int NDim (void) const; + double DoEval (const double * xin) const; + ROOT::Math::IBaseFunctionMultiDim * Clone (void) const; + private: + const XSecAlgorithmI * fModel; + const Interaction * fInteraction; + double fScale; // can set to -1. for use with GSL minimizer +}; + +//..................................................................................... +// +// genie::utils::gsl::d2Xsec_dn1dn2dn3_E +// A 3-D cross section function: d2xsec/dn1dn2dn3 = f(n1,n2,n3)|(fixed E) +// +class d2Xsec_dn1dn2dn3_E: public ROOT::Math::IBaseFunctionMultiDim +{ + public: + d2Xsec_dn1dn2dn3_E(const XSecAlgorithmI * m, const Interaction * i, double scale=1. ); + ~d2Xsec_dn1dn2dn3_E(); + // ROOT::Math::IBaseFunctionMultiDim interface + unsigned int NDim (void) const; + double DoEval (const double * xin) const; + ROOT::Math::IBaseFunctionMultiDim * Clone (void) const; + private: + const XSecAlgorithmI * fModel; + const Interaction * fInteraction; + double fScale; // can set to -1. for use with GSL minimizer +}; + + } // gsl namespace } // utils namespace } // genie namespace diff --git a/src/Physics/XSectionIntegration/LinkDef.h b/src/Physics/XSectionIntegration/LinkDef.h index a807651a5..89440fead 100644 --- a/src/Physics/XSectionIntegration/LinkDef.h +++ b/src/Physics/XSectionIntegration/LinkDef.h @@ -25,5 +25,7 @@ #pragma link C++ class genie::utils::gsl::d3Xsec_dOmegaldThetapi; #pragma link C++ class genie::utils::gsl::dXSec_dElep_AR; #pragma link C++ class genie::utils::gsl::dXSec_Log_Wrapper; +#pragma link C++ class genie::utils::gsl::d2Xsec_dn1dn2_E; +#pragma link C++ class genie::utils::gsl::d2Xsec_dn1dn2dn3_E; #endif diff --git a/src/Tools/Geometry/ROOTGeomAnalyzer.cxx b/src/Tools/Geometry/ROOTGeomAnalyzer.cxx index 082536a4e..3b7245029 100644 --- a/src/Tools/Geometry/ROOTGeomAnalyzer.cxx +++ b/src/Tools/Geometry/ROOTGeomAnalyzer.cxx @@ -203,6 +203,42 @@ const PathLengthList & ROOTGeomAnalyzer::ComputePathLengths( return *fCurrPathLengthList; } +//___________________________________________________________________________ +std::vector< std::pair > ROOTGeomAnalyzer::ComputeMatLengths( + const TLorentzVector & x, const TLorentzVector & p) +{ + + // if trimming configure with neutrino ray's info + if ( fGeomVolSelector ) { + fGeomVolSelector->SetCurrentRay(x,p); + fGeomVolSelector->SetSI2Local(1/this->LengthUnits()); + } + + TVector3 udir = p.Vect().Unit(); // unit vector along direction + TVector3 pos = x.Vect(); // initial position + this->SI2Local(pos); // SI -> curr geom units + + if (!fMasterToTopIsIdentity) { + this->Master2Top(pos); // transform position (master -> top) + this->Master2TopDir(udir); // transform direction (master -> top) + } + + this->SwimOnce(pos,udir); + + std::vector> MatLengthList; + + const PathSegmentList::PathSegmentV_t& segments = fCurrPathSegmentList->GetPathSegmentV(); + + PathSegmentList::PathSegVCItr_t sitr; + for ( sitr = segments.begin(); sitr != segments.end(); ++sitr) { + const PathSegment& seg = *sitr; + double pl = seg.GetSummedStepRange(); + MatLengthList.push_back(std::make_pair(pl,seg.fMaterial)); + } + + return MatLengthList; +} + //___________________________________________________________________________ const TVector3 & ROOTGeomAnalyzer::GenerateVertex( const TLorentzVector & x, const TLorentzVector & p, int tgtpdg) diff --git a/src/Tools/Geometry/ROOTGeomAnalyzer.h b/src/Tools/Geometry/ROOTGeomAnalyzer.h index 60592145a..373af8e0a 100644 --- a/src/Tools/Geometry/ROOTGeomAnalyzer.h +++ b/src/Tools/Geometry/ROOTGeomAnalyzer.h @@ -66,9 +66,16 @@ public : virtual const PathLengthList & ComputePathLengths(const TLorentzVector & x, const TLorentzVector & p); + + virtual std::vector< std::pair > ComputeMatLengths(const TLorentzVector & x, + const TLorentzVector & p); + virtual const TVector3 & GenerateVertex(const TLorentzVector & x, const TLorentzVector & p, int tgtpdg); + virtual int GetTargetPdgCode (const TGeoMaterial * const m) const; + virtual int GetTargetPdgCode (const TGeoMixture * const m, int ielement) const; + /// set geometry driver's configuration options virtual void SetScannerNPoints (int np) { fNPoints = np; } /* box scanner */ @@ -123,8 +130,6 @@ public : virtual void Load (TGeoManager * gm); virtual void BuildListOfTargetNuclei (void); - virtual int GetTargetPdgCode (const TGeoMaterial * const m) const; - virtual int GetTargetPdgCode (const TGeoMixture * const m, int ielement) const; virtual double GetWeight (const TGeoMaterial * mat, int pdgc); virtual double GetWeight (const TGeoMixture * mixt, int pdgc); virtual double GetWeight (const TGeoMixture * mixt, int ielement, int pdgc); diff --git a/src/make/Make.include b/src/make/Make.include index 1552afc49..5c66cfbac 100644 --- a/src/make/Make.include +++ b/src/make/Make.include @@ -527,6 +527,7 @@ CPP_INCLUDES := \ $(LOG_INCLUDES) \ $(ROOT_INCLUDES) \ $(LHAPDF_INCLUDES) \ + $(APFEL_INCLUDES) \ $(GSL_INCLUDES) \ $(GENIE_INCLUDES) @@ -542,6 +543,7 @@ LIBRARIES := $(SYSLIBS) \ $(ROOT_LIBRARIES) \ $(PYTHIA6_LIBRARIES) \ $(LHAPDF_LIBRARIES) \ + $(APFEL_LIBRARIES) \ $(XML_LIBRARIES) \ $(LOG_LIBRARIES) \ $(GSL_LIBRARIES) \ diff --git a/src/scripts/setup/genie-config b/src/scripts/setup/genie-config index ff633a8d7..55b2e83ea 100755 --- a/src/scripts/setup/genie-config +++ b/src/scripts/setup/genie-config @@ -35,7 +35,7 @@ fmwk_libs=" -lGFwMsg -lGFwReg -lGFwAlg -lGFwInt -lGFwGHEP -lGFwNum -lGFwUtl -lGF # (standard) phys_libs=" -lGPhXSIg -lGPhPDF -lGPhNuclSt -lGPhCmn -lGPhDcy -lGPhHadTransp -lGPhHadnz -lGPhDeEx \ -lGPhAMNGXS -lGPhAMNGEG -lGPhChmXS -lGPhCohXS -lGPhCohEG -lGPhDISXS -lGPhDISEG \ - -lGPhDfrcXS -lGPhDfrcEG -lGPhGlwResXS -lGPhGlwResEG -lGPhIBDXS -lGPhIBDEG -lGPhHadTens \ + -lGPhDfrcXS -lGPhDfrcEG -lGPhHELptnXS -lGPhHELptnEG -lGPhIBDXS -lGPhIBDEG -lGPhHadTens \ -lGPhMNucXS -lGPhMNucEG -lGPhMEL -lGPhNuElXS -lGPhNuElEG \ -lGPhQELXS -lGPhQELEG -lGPhResXS -lGPhResEG -lGPhStrXS -lGPhStrEG -lGPhHEDISXS -lGPhHEDISEG" # (optional)