diff --git a/README.md b/README.md index e27d9417c3..e0dedcbb94 100644 --- a/README.md +++ b/README.md @@ -1,6 +1,5 @@ # The GENIE Event Generator -
@@ -8,14 +7,14 @@

- The GENIE Generator is a leading simulation tool used by nearly all modern neutrino experiments. - It features a modular framework with state-of-the-art physics for neutrino and charged-lepton interactions, - and several BSM channels. It incorporates results from GENIE’s global data analysis and includes multiple tuned models. - GENIE supports all neutrino types, targets, and energy scales from MeV to PeV, + The GENIE Generator is a leading simulation tool used by nearly all modern neutrino experiments. + It features a modular framework with state-of-the-art physics for neutrino and charged-lepton interactions, + and several BSM channels. It incorporates results from GENIE’s global data analysis and includes multiple tuned models. + GENIE supports all neutrino types, targets, and energy scales from MeV to PeV, and provides tools for flux handling, geometry, event generation, and reweighting.

- For more information, visit + For more information, visit http://genie-mc.org | https://genie-mc.github.io.

@@ -27,7 +26,7 @@ Luis Alvarez-Ruso (*IFIC*), Costas Andreopoulos (*Liverpool*), Adi Ashkenazi (*Tel Aviv*), Joshua Barrow (*Minnesota*), Steve Dytman (*Pittsburgh*), Hugh Gallagher (*Tufts*), Alfonso Andres Garcia Soto (*IFIC*), Steven Gardiner (*Fermilab*), Matan Goldenberg (*Tel Aviv*), Robert Hatcher (*Fermilab*), Or Hen (*MIT*), Igor Kakorin (*JINR*), Konstantin Kuzmin (*ITEP and JINR*), Liang Liu (*Fermilab*), Xianguo Lu (*Warwick*), Anselmo Meregaglia (*Bordeaux, CNRS/IN2P3*), Vadim Naumov (*JINR*), Afroditi Papadopoulou (*Argonne*), Gabriel Perdue (*Fermilab*), Komninos-John Plows (*Oxford*), Marco Roda (*Liverpool*), Alon Sportes (*Tel Aviv*), Júlia Tena Vidal (*Tel Aviv*), Jeremy Wolcott (*Tufts*), Qiyu Yan (*UCAS and Warwick*) -**Past authors:** Christopher Barry (*Liverpool*), Steve Dennis (*Liverpool*), Walter Giele (*Fermilab*), Timothy Hobbs (*Fermilab*), Libo Jiang (*Pittsburgh*), Rhiannon Jones (*Liverpool*), Weijun Li (*Oxford*), Donna Naples (*Pittsburgh*), Beth Slater (*Liverpool*), Noah Steinberg (*Fermilab*), Vladyslav Syrotenko (*Tufts*), Julia Yarba (*Fermilab*) +**Past authors:** Christopher Barry (*Liverpool*), Steve Dennis (*Liverpool*), Walter Giele (*Fermilab*), Timothy Hobbs (*Fermilab*), Libo Jiang (*Pittsburgh*), Rhiannon Jones (*Liverpool*), Weijun Li (*Oxford*), Donna Naples (*Pittsburgh*), Beth Slater (*Liverpool*), Noah Steinberg (*Fermilab*), Vladyslav Syrotenko (*Tufts*), Julia Yarba (*Fermilab*) For more details on the GENIE collaboration please visit [this page](https://genie-mc.github.io/collaboration.html). @@ -54,7 +53,7 @@ DOIs for recent releases of the GENIE Event Generator are listed below: ## Physics & User manual -For installation and usage information, as well as information on the GENIE framework, event generator modules and tuning, +For installation and usage information, as well as information on the GENIE framework, event generator modules and tuning, see the [latest version of the GENIE Physics & User Manual](https://www.overleaf.com/read/rqmbqzzvsvmb#5ab475), originally published as arXiv:1510.05494. ## Public releases and physics tunes @@ -63,7 +62,7 @@ For a list of public releases and a summary information, see [this page](https:/ A list of model configurations and tunes supported in each release is maintained [here](https://genie-mc.github.io/tunes.html). Details on the naming conventions for releases, model configurations and tunes can be found [here](https://genie-mc.github.io/naming_conventions.html). -[Recent publications and talks](https://genie-mc.github.io/pub.html) +[Recent publications and talks](https://genie-mc.github.io/pub.html) by GENIE authors highlight key modeling advances and results from our global analysis of scattering data. ## Contribution guidelines diff --git a/config/CommonParam.xml b/config/CommonParam.xml index 2d24e8e061..2134a0b17c 100644 --- a/config/CommonParam.xml +++ b/config/CommonParam.xml @@ -276,18 +276,21 @@ Or changing the name of this parameter set ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ Set the coupling of dark matter to nucleons --> - 0.1 + + 1.0 1.0 - 1.0 + -1.0 1.0 + 1.0 - 1.0 + -1.0 1.0 - 1.0 + -1.0 1.0 - 1.0 + -1.0 1.0 - 1.0 + -1.0 + 1.0 1.0 1.441 diff --git a/config/DMRESHadronicSystemGenerator.xml b/config/DMRESHadronicSystemGenerator.xml new file mode 100644 index 0000000000..05c1dfffea --- /dev/null +++ b/config/DMRESHadronicSystemGenerator.xml @@ -0,0 +1,22 @@ + + + + + + + + genie::UnstableParticleDecayer/BeforeHadronTransport + + + + diff --git a/config/DMRESInteractionListGenerator.xml b/config/DMRESInteractionListGenerator.xml new file mode 100644 index 0000000000..77f20ea85f --- /dev/null +++ b/config/DMRESInteractionListGenerator.xml @@ -0,0 +1,24 @@ + + + + + + + + Resonances + + + + true + + + diff --git a/config/DMRESKinematicsGenerator.xml b/config/DMRESKinematicsGenerator.xml new file mode 100644 index 0000000000..23676f1f1e --- /dev/null +++ b/config/DMRESKinematicsGenerator.xml @@ -0,0 +1,35 @@ + + + + + + + + NonResBackground + + + + + + 1.400 + 0.500 + + + + 1.400 + + + diff --git a/config/DMRESOutgoingDarkGenerator.xml b/config/DMRESOutgoingDarkGenerator.xml new file mode 100644 index 0000000000..549e3c7a23 --- /dev/null +++ b/config/DMRESOutgoingDarkGenerator.xml @@ -0,0 +1,19 @@ + + + + + + + + Lepton + + + + diff --git a/config/DMRESPXSec.xml b/config/DMRESPXSec.xml new file mode 100644 index 0000000000..767076139a --- /dev/null +++ b/config/DMRESPXSec.xml @@ -0,0 +1,64 @@ + + + + + + + + WeakInt,NonResBackground,CKM,FermiGas,BoostedDarkMatter + + genie::DMRESXSecFast/Default + + 0 + 1.000 + 1.05 + 0.76338 + 1.120 + 0.840 + + + + WeakInt,NonResBackground,CKM,FermiGas,BoostedDarkMatter + + genie::DMRESXSecFast/Default + + 2 + 1.000 + 1.05 + 0.76338 + 1.120 + 0.840 + + + + diff --git a/config/DMRESXSec.xml b/config/DMRESXSec.xml new file mode 100644 index 0000000000..0bf5b2961e --- /dev/null +++ b/config/DMRESXSec.xml @@ -0,0 +1,19 @@ + + + + + + + + + Resonances + 500 + adaptive + 1000000000 + 1e-9 + + false + + diff --git a/config/DMRESXSecFast.xml b/config/DMRESXSecFast.xml new file mode 100644 index 0000000000..10ea7583f7 --- /dev/null +++ b/config/DMRESXSecFast.xml @@ -0,0 +1,19 @@ + + + + + + + + + Resonances + 500 + adaptive + 1000000000 + 1e-9 + + false + + diff --git a/config/EventGenerator.xml b/config/EventGenerator.xml index 9187896319..a35feaf690 100644 --- a/config/EventGenerator.xml +++ b/config/EventGenerator.xml @@ -114,7 +114,6 @@ XSecModel alg Yes Cross section model used at the thread genie::DMDISInteractionListGenerator/DM-Default - @@ -689,6 +688,22 @@ XSecModel alg Yes Cross section model used at the thread genie::DMEInteractionListGenerator/DME + + + 10 + genie::InitialStateAppender/Default + genie::VertexGenerator/Default + genie::FermiMover/Default + genie::DMRESKinematicsGenerator/RES + genie::DMRESOutgoingDarkGenerator/Default + genie::DMRESHadronicSystemGenerator/Default + genie::NucDeExcitationSim/Default + genie::HadronTransporter/Default + genie::NucBindEnergyAggregator/Default + genie::UnstableParticleDecayer/AfterHadronTransport + genie::DMRESInteractionListGenerator/DM-Default + + @@ -847,7 +862,7 @@ XSecModel alg Yes Cross section model used at the thread - + 4 @@ -858,7 +873,7 @@ XSecModel alg Yes Cross section model used at the thread genie::DummyHNLInteractionListGenerator/Default - + 3 @@ -868,5 +883,5 @@ XSecModel alg Yes Cross section model used at the thread genie::NormInteractionListGenerator/Default genie::NormXSec/Default - + diff --git a/config/EventGeneratorListAssembler.xml b/config/EventGeneratorListAssembler.xml index 8f98d98d40..b2e6e78495 100644 --- a/config/EventGeneratorListAssembler.xml +++ b/config/EventGeneratorListAssembler.xml @@ -7,31 +7,31 @@ .......................................................................................................................... NOTE: -The GENIE authors would like to caution users against using anything but a comprehensive +The GENIE authors would like to caution users against using anything but a comprehensive GENIE mode (`Default' setting below) for physics studies. No detector measures generator-level reaction modes like CCQE or NCRES. -Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+}, {0mu-, 1pi0}, +Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+}, {0mu-, 1pi0}, {1 track, 1 shower}, {1 mu-like ring} etc depending on granularity, thresholds and PID capabilities. No final state / topology is a proxy for any particular reaction mode. Intranuclear re-scattering in particular causes significant migration between states (see Table 8.1 in the Physics and User manual). -The Default configuration depends on the tune and it is defined in TuneGeneratorList.xml, +The Default configuration depends on the tune and it is defined in TuneGeneratorList.xml, which is defined in each Tune subdirectory. Examples: - {1mu-,0pi} is mostly numuCCQE, but can also be caused by numu resonance production followed by pion absorption. - numuCCQE gives mostly {1mu-,0pi} but occasionaly can give {1mu-,1pi} if the recoil nucleon re-interacts. -- NC1pi0 final states can be caused by all a) NC elastic followed by nucleon rescattering, +- NC1pi0 final states can be caused by all a) NC elastic followed by nucleon rescattering, b) NC resonance neutrino-production, c) NC non-resonance background, d) low-W NC DIS e) NC coherent scattering. Each source contributes differently in the pion momentum distribution. -We also recommend that you treat the generator-level reaction modes largely as as an internal degree of freedom. +We also recommend that you treat the generator-level reaction modes largely as as an internal degree of freedom. Consequently, try to define your selection efficiencies and purities in terms of observable final states _only_. For example, if you define as `numuCCQE-like' := {1mu-,0pi} then define your selection efficiency as: efficiency = (`numuCCQE-like' events passing cuts) / (all true {1mu-,0pi} events) -and not as: +and not as: efficiency = (`numuCCQE-like' events passing cuts) / (all true numuCCQE events) .......................................................................................................................... @@ -47,44 +47,44 @@ Generator-%d alg No --> - + 1 genie::EventGenerator/AM-NUGAMMA - + 1 genie::EventGenerator/CEvNS - + 2 genie::EventGenerator/DFR-CC genie::EventGenerator/DFR-NC - + 1 genie::EventGenerator/DFR-CC - + 1 genie::EventGenerator/DFR-NC - + 2 genie::EventGenerator/QEL-CC-CHARM genie::EventGenerator/DIS-CC-CHARM - + 1 genie::EventGenerator/DIS-CC-CHARM - + 1 genie::EventGenerator/QEL-CC-CHARM @@ -103,48 +103,48 @@ Generator-%d alg No 1 genie::EventGenerator/IBD - - + + 2 genie::EventGenerator/IMD genie::EventGenerator/IMD-ANH - + 1 genie::EventGenerator/NUE-EL - + 3 genie::EventGenerator/NUE-EL genie::EventGenerator/IMD genie::EventGenerator/IMD-ANH - + 2 genie::EventGenerator/QEL-CC genie::EventGenerator/QEL-NC - + 1 genie::EventGenerator/QEL-CC - + 1 genie::EventGenerator/QEL-NC - + 1 genie::EventGenerator/MEC-CC - + 1 genie::EventGenerator/MEC-NC @@ -156,13 +156,13 @@ Generator-%d alg No - + 2 genie::EventGenerator/MEC-CC genie::EventGenerator/QEL-CC - + 2 genie::EventGenerator/MEC-NC genie::EventGenerator/QEL-NC @@ -174,55 +174,55 @@ Generator-%d alg No genie::EventGenerator/QEL-EM - + 2 genie::EventGenerator/RES-CC genie::EventGenerator/RES-NC - + 1 genie::EventGenerator/RES-CC - + 1 genie::EventGenerator/RES-NC - + 2 genie::EventGenerator/COH-CC-PION genie::EventGenerator/COH-NC-PION - + 1 genie::EventGenerator/COH-CC-PION - + 1 genie::EventGenerator/COH-NC-PION - + 2 genie::EventGenerator/DIS-CC genie::EventGenerator/DIS-NC - + 1 genie::EventGenerator/DIS-CC - + 1 genie::EventGenerator/DIS-NC - + 6 genie::EventGenerator/QEL-NC genie::EventGenerator/RES-NC @@ -232,7 +232,7 @@ Generator-%d alg No genie::EventGenerator/DFR-NC - + 9 genie::EventGenerator/QEL-CC genie::EventGenerator/RES-CC @@ -245,7 +245,7 @@ Generator-%d alg No genie::EventGenerator/QEL-CC-LAMBDA - + 5 genie::EventGenerator/QEL-CC genie::EventGenerator/RES-CC @@ -266,26 +266,26 @@ Generator-%d alg No genie::EventGenerator/NUE-EL - - + 1 genie::EventGenerator/QEL-EM - + 1 genie::EventGenerator/MEC-EM - + 1 genie::EventGenerator/RES-EM - + 1 genie::EventGenerator/DIS-EM @@ -293,7 +293,7 @@ Generator-%d alg No - + 4 genie::EventGenerator/QEL-EM genie::EventGenerator/RES-EM @@ -301,11 +301,11 @@ Generator-%d alg No genie::EventGenerator/MEC-EM - - + 2 genie::EventGenerator/RES-CC genie::EventGenerator/DIS-CC @@ -313,7 +313,7 @@ Generator-%d alg No - + 9 genie::EventGenerator/QEL-CC @@ -323,8 +323,8 @@ Generator-%d alg No genie::EventGenerator/IMD-ANH genie::EventGenerator/QEL-CC-LAMBDA genie::EventGenerator/DIS-CC-SINGLEK - genie::EventGenerator/QEL-CC-CHARM - genie::EventGenerator/DIS-CC-CHARM + genie::EventGenerator/QEL-CC-CHARM + genie::EventGenerator/DIS-CC-CHARM @@ -379,7 +379,7 @@ Generator-%d alg No genie::EventGenerator/QEL-CC-LAMBDA genie::EventGenerator/DIS-CC-SINGLEK genie::EventGenerator/MEC-CC - genie::EventGenerator/MEC-NC + genie::EventGenerator/MEC-NC genie::EventGenerator/QEL-CC-CHARM genie::EventGenerator/DIS-CC-CHARM @@ -397,6 +397,11 @@ Generator-%d alg No genie::EventGenerator/DMDIS + + 1 + genie::EventGenerator/DMRES + + 1 genie::EventGenerator/DME @@ -409,10 +414,11 @@ Generator-%d alg No - 3 + 4 genie::EventGenerator/DMEL genie::EventGenerator/DMDIS genie::EventGenerator/DME + genie::EventGenerator/DMRES @@ -420,10 +426,16 @@ Generator-%d alg No genie::EventGenerator/COHDNu + + 3 + genie::EventGenerator/DMEL + genie::EventGenerator/DMDIS + genie::EventGenerator/DME + - + 13 genie::EventGenerator/HEDIS-CC genie::EventGenerator/HEDIS-NC @@ -440,7 +452,7 @@ Generator-%d alg No genie::EventGenerator/PhotonCOH - + 11 genie::EventGenerator/GLRES-Mu genie::EventGenerator/GLRES-Tau @@ -455,7 +467,7 @@ Generator-%d alg No genie::EventGenerator/PhotonCOH - + 6 genie::EventGenerator/HEDIS-CC genie::EventGenerator/HEDIS-NC @@ -465,7 +477,7 @@ Generator-%d alg No genie::EventGenerator/GLRES-Had - + 5 genie::EventGenerator/HEDIS-CC genie::EventGenerator/GLRES-Mu @@ -474,23 +486,23 @@ Generator-%d alg No genie::EventGenerator/GLRES-Had - + 2 genie::EventGenerator/HEDIS-CC genie::EventGenerator/HEDIS-NC - + 1 genie::EventGenerator/HEDIS-CC - + 1 genie::EventGenerator/HEDIS-NC - + 4 genie::EventGenerator/GLRES-Mu genie::EventGenerator/GLRES-Tau @@ -498,43 +510,43 @@ Generator-%d alg No 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 @@ -542,41 +554,40 @@ Generator-%d alg No 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 - + 1 genie::EventGenerator/NORM - - + + 2 genie::EventGenerator/QEL-CC genie::EventGenerator/NORM - diff --git a/config/GDM18_00a/ModelConfiguration.xml b/config/GDM18_00a/ModelConfiguration.xml index e562596913..3b37566680 100644 --- a/config/GDM18_00a/ModelConfiguration.xml +++ b/config/GDM18_00a/ModelConfiguration.xml @@ -4,7 +4,7 @@ - + - - - - + - + genie::ReinSehgalRESPXSec/NoPauliBlock genie::ReinSehgalRESPXSec/NoPauliBlock - genie::ReinSehgalRESPXSec/Default - + genie::ReinSehgalRESPXSec/EM-NoPauliBlock + genie::KNOTunedQPMDISPXSec/Default genie::KNOTunedQPMDISPXSec/Default genie::KNOTunedQPMDISPXSec/Default - genie::ReinSehgalCOHPiPXSec/Default genie::ReinSehgalCOHPiPXSec/Default @@ -138,6 +138,8 @@ University of Liverpool genie::AhrensDMELPXSec/Velocity0 genie::QPMDMDISPXSec/Velocity0 genie::DMElectronPXSec/Velocity0 + genie::DMRESPXSec/Velocity0 + diff --git a/config/GDM18_00a/TuneGeneratorList.xml b/config/GDM18_00a/TuneGeneratorList.xml index 60161680c5..b7b9946f45 100644 --- a/config/GDM18_00a/TuneGeneratorList.xml +++ b/config/GDM18_00a/TuneGeneratorList.xml @@ -9,10 +9,10 @@ This xml files contains the default configuration for the tune .......................................................................................................................... NOTE: -The GENIE authors would like to caution users against using anything but a comprehensive +The GENIE authors would like to caution users against using anything but a comprehensive GENIE mode (`Default' setting below) for physics studies. No detector measures generator-level reaction modes like CCQE or NCRES. -Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+}, {0mu-, 1pi0}, +Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+}, {0mu-, 1pi0}, {1 track, 1 shower}, {1 mu-like ring} etc depending on granularity, thresholds and PID capabilities. No final state / topology is a proxy for any particular reaction mode. Intranuclear re-scattering in particular causes significant migration between states @@ -21,16 +21,16 @@ Intranuclear re-scattering in particular causes significant migration between st Examples: - {1mu-,0pi} is mostly numuCCQE, but can also be caused by numu resonance production followed by pion absorption. - numuCCQE gives mostly {1mu-,0pi} but occasionaly can give {1mu-,1pi} if the recoil nucleon re-interacts. -- NC1pi0 final states can be caused by all a) NC elastic followed by nucleon rescattering, +- NC1pi0 final states can be caused by all a) NC elastic followed by nucleon rescattering, b) NC resonance neutrino-production, c) NC non-resonance background, d) low-W NC DIS e) NC coherent scattering. Each source contributes differently in the pion momentum distribution. -We also recommend that you treat the generator-level reaction modes largely as as an internal degree of freedom. +We also recommend that you treat the generator-level reaction modes largely as as an internal degree of freedom. Consequently, try to define your selection efficiencies and purities in terms of observable final states _only_. For example, if you define as `numuCCQE-like' := {1mu-,0pi} then define your selection efficiency as: efficiency = (`numuCCQE-like' events passing cuts) / (all true {1mu-,0pi} events) -and not as: +and not as: efficiency = (`numuCCQE-like' events passing cuts) / (all true numuCCQE events) .......................................................................................................................... @@ -46,12 +46,42 @@ Generator-%d alg No --> - + + 4 + genie::EventGenerator/DMEL + genie::EventGenerator/DMDIS + genie::EventGenerator/DME + genie::EventGenerator/DMRES + + + 3 genie::EventGenerator/DMEL - genie::EventGenerator/DMDIS - genie::EventGenerator/DME + genie::EventGenerator/DMDIS + genie::EventGenerator/DME + + + + 18 + genie::EventGenerator/QEL-CC + genie::EventGenerator/QEL-NC + genie::EventGenerator/RES-CC + genie::EventGenerator/RES-NC + genie::EventGenerator/DIS-CC + genie::EventGenerator/DIS-NC + genie::EventGenerator/COH-CC-PION + genie::EventGenerator/COH-NC-PION + genie::EventGenerator/DIS-CC-CHARM + genie::EventGenerator/QEL-CC-CHARM + genie::EventGenerator/NUE-EL + genie::EventGenerator/IMD + genie::EventGenerator/IMD-ANH + genie::EventGenerator/DFR-CC + genie::EventGenerator/DFR-NC + genie::EventGenerator/QEL-CC-LAMBDA + genie::EventGenerator/MEC-CC + genie::EventGenerator/MEC-NC - + diff --git a/config/GDM18_00b/ModelConfiguration.xml b/config/GDM18_00b/ModelConfiguration.xml index 3d7c2ff72e..0dee032efe 100644 --- a/config/GDM18_00b/ModelConfiguration.xml +++ b/config/GDM18_00b/ModelConfiguration.xml @@ -4,7 +4,7 @@ - + - - - - + - + genie::ReinSehgalRESPXSec/NoPauliBlock genie::ReinSehgalRESPXSec/NoPauliBlock genie::ReinSehgalRESPXSec/Default - + genie::KNOTunedQPMDISPXSec/Default genie::KNOTunedQPMDISPXSec/Default genie::KNOTunedQPMDISPXSec/Default - genie::ReinSehgalCOHPiPXSec/Default genie::ReinSehgalCOHPiPXSec/Default @@ -138,6 +138,8 @@ University of Liverpool genie::AhrensDMELPXSec/Velocity2 genie::QPMDMDISPXSec/Velocity2 genie::DMElectronPXSec/Velocity2 + genie::DMRESPXSec/Velocity2 + diff --git a/config/GDM18_00b/TuneGeneratorList.xml b/config/GDM18_00b/TuneGeneratorList.xml index 60161680c5..3805e309ab 100644 --- a/config/GDM18_00b/TuneGeneratorList.xml +++ b/config/GDM18_00b/TuneGeneratorList.xml @@ -9,10 +9,10 @@ This xml files contains the default configuration for the tune .......................................................................................................................... NOTE: -The GENIE authors would like to caution users against using anything but a comprehensive +The GENIE authors would like to caution users against using anything but a comprehensive GENIE mode (`Default' setting below) for physics studies. No detector measures generator-level reaction modes like CCQE or NCRES. -Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+}, {0mu-, 1pi0}, +Detectors measure final states / topologies like {1mu-,0pi}, {1mu-,1pi+}, {0mu-, 1pi0}, {1 track, 1 shower}, {1 mu-like ring} etc depending on granularity, thresholds and PID capabilities. No final state / topology is a proxy for any particular reaction mode. Intranuclear re-scattering in particular causes significant migration between states @@ -21,16 +21,16 @@ Intranuclear re-scattering in particular causes significant migration between st Examples: - {1mu-,0pi} is mostly numuCCQE, but can also be caused by numu resonance production followed by pion absorption. - numuCCQE gives mostly {1mu-,0pi} but occasionaly can give {1mu-,1pi} if the recoil nucleon re-interacts. -- NC1pi0 final states can be caused by all a) NC elastic followed by nucleon rescattering, +- NC1pi0 final states can be caused by all a) NC elastic followed by nucleon rescattering, b) NC resonance neutrino-production, c) NC non-resonance background, d) low-W NC DIS e) NC coherent scattering. Each source contributes differently in the pion momentum distribution. -We also recommend that you treat the generator-level reaction modes largely as as an internal degree of freedom. +We also recommend that you treat the generator-level reaction modes largely as as an internal degree of freedom. Consequently, try to define your selection efficiencies and purities in terms of observable final states _only_. For example, if you define as `numuCCQE-like' := {1mu-,0pi} then define your selection efficiency as: efficiency = (`numuCCQE-like' events passing cuts) / (all true {1mu-,0pi} events) -and not as: +and not as: efficiency = (`numuCCQE-like' events passing cuts) / (all true numuCCQE events) .......................................................................................................................... @@ -46,12 +46,41 @@ Generator-%d alg No --> - + + 4 + genie::EventGenerator/DMEL + genie::EventGenerator/DMDIS + genie::EventGenerator/DME + genie::EventGenerator/DMRES + + + 3 genie::EventGenerator/DMEL - genie::EventGenerator/DMDIS - genie::EventGenerator/DME + genie::EventGenerator/DMDIS + genie::EventGenerator/DME - + + 18 + genie::EventGenerator/QEL-CC + genie::EventGenerator/QEL-NC + genie::EventGenerator/RES-CC + genie::EventGenerator/RES-NC + genie::EventGenerator/DIS-CC + genie::EventGenerator/DIS-NC + genie::EventGenerator/COH-CC-PION + genie::EventGenerator/COH-NC-PION + genie::EventGenerator/DIS-CC-CHARM + genie::EventGenerator/QEL-CC-CHARM + genie::EventGenerator/NUE-EL + genie::EventGenerator/IMD + genie::EventGenerator/IMD-ANH + genie::EventGenerator/DFR-CC + genie::EventGenerator/DFR-NC + genie::EventGenerator/QEL-CC-LAMBDA + genie::EventGenerator/MEC-CC + genie::EventGenerator/MEC-NC + + diff --git a/config/Messenger.xml b/config/Messenger.xml index ddc0283797..68723e2a0d 100644 --- a/config/Messenger.xml +++ b/config/Messenger.xml @@ -86,6 +86,13 @@ NOTICE WARN NOTICE + INFO + INFO + WARN + WARN + NOTICE + NOTICE + WARN WARN NOTICE INFO @@ -231,7 +238,7 @@ WARN WARN ERROR - FATAL + FATAL NOTICE NOTICE diff --git a/config/Messenger_laconic.xml b/config/Messenger_laconic.xml index 11d274bc63..68ea0e1163 100644 --- a/config/Messenger_laconic.xml +++ b/config/Messenger_laconic.xml @@ -21,7 +21,7 @@ same stream is listed twice with conflicting priority then the one found last is used --> - WARN + WARN WARN WARN WARN @@ -79,6 +79,13 @@ WARN WARN WARN + INFO + INFO + WARN + WARN + NOTICE + NOTICE + WARN WARN WARN WARN diff --git a/config/Messenger_rambling.xml b/config/Messenger_rambling.xml index ccc7743ea0..96086459dc 100644 --- a/config/Messenger_rambling.xml +++ b/config/Messenger_rambling.xml @@ -80,6 +80,13 @@ NOTICE NOTICE NOTICE + INFO + INFO + WARN + WARN + NOTICE + NOTICE + WARN NOTICE INFO INFO @@ -93,7 +100,7 @@ WARN NOTICE NOTICE - NOTICE + NOTICE INFO INFO WARN diff --git a/config/Messenger_whisper.xml b/config/Messenger_whisper.xml index 6a815610cc..82732e5bd5 100644 --- a/config/Messenger_whisper.xml +++ b/config/Messenger_whisper.xml @@ -80,6 +80,13 @@ NOTICE WARN NOTICE + INFO + INFO + WARN + WARN + NOTICE + NOTICE + WARN FATAL FATAL FATAL diff --git a/config/RSHelicityAmplModelDMn.xml b/config/RSHelicityAmplModelDMn.xml new file mode 100644 index 0000000000..49ca01898d --- /dev/null +++ b/config/RSHelicityAmplModelDMn.xml @@ -0,0 +1,18 @@ + + + + + + + + BoostedDarkMatter + + + diff --git a/config/RSHelicityAmplModelDMp.xml b/config/RSHelicityAmplModelDMp.xml new file mode 100644 index 0000000000..4ba3606acd --- /dev/null +++ b/config/RSHelicityAmplModelDMp.xml @@ -0,0 +1,18 @@ + + + + + + + + BoostedDarkMatter + + + diff --git a/config/master_config.xml b/config/master_config.xml index c4986d8c87..0dabf07fe2 100644 --- a/config/master_config.xml +++ b/config/master_config.xml @@ -27,6 +27,7 @@ QELPrimaryLeptonGenerator.xml DISPrimaryLeptonGenerator.xml RESPrimaryLeptonGenerator.xml + DMRESOutgoingDarkGenerator.xml NuEPrimaryLeptonGenerator.xml DMEOutgoingDarkGenerator.xml COHPrimaryLeptonGenerator.xml @@ -45,6 +46,7 @@ SKKinematicsGenerator.xml DMELKinematicsGenerator.xml DMDISKinematicsGenerator.xml + DMRESKinematicsGenerator.xml HELeptonKinematicsGenerator.xml HEDISKinematicsGenerator.xml IBDHadronicSystemGenerator.xml @@ -53,6 +55,7 @@ DFRHadronicSystemGenerator.xml DISHadronicSystemGenerator.xml RESHadronicSystemGenerator.xml + DMRESHadronicSystemGenerator.xml RSPPHadronicSystemGenerator.xml SKHadronicSystemGenerator.xml NuETargetRemnantGenerator.xml @@ -93,6 +96,7 @@ SKInteractionListGenerator.xml DMELInteractionListGenerator.xml DMDISInteractionListGenerator.xml + DMRESInteractionListGenerator.xml IBDInteractionListGenerator.xml CEvNSInteractionListGenerator.xml COHDNuInteractionListGenerator.xml @@ -156,6 +160,8 @@ RSHelicityAmplModelNCn.xml RSHelicityAmplModelEMp.xml RSHelicityAmplModelEMn.xml + RSHelicityAmplModelDMn.xml + RSHelicityAmplModelDMp.xml EngelFormFactor.xml @@ -172,6 +178,7 @@ AlamSimoAtharVacasSKXSec.xml IMDXSec.xml RESXSec.xml + DMRESXSec.xml MECXSec.xml NuElectronXSec.xml DMElectronXSec.xml @@ -209,6 +216,7 @@ PattonCEvNSPXSec.xml NuElectronPXSec.xml DMElectronPXSec.xml + DMRESPXSec.xml GLRESPXSec.xml HENuElPXSec.xml PhotonCOHPXSec.xml @@ -232,6 +240,7 @@ Default.xml Default.xml ReinSehgalRESXSecFast.xml + DMRESXSecFast.xml BertuzzoDNuCOHPXSec.xml HybridXSecAlgorithm.xml HEDISPXSec.xml diff --git a/src/Apps/gAtmoEvGen.cxx b/src/Apps/gAtmoEvGen.cxx index 1872f3170e..57d74dbdfc 100644 --- a/src/Apps/gAtmoEvGen.cxx +++ b/src/Apps/gAtmoEvGen.cxx @@ -376,6 +376,7 @@ int main(int argc, char** argv) * option. */ mcj_driver->ForceSingleProbScale(); + // initialize an ntuple writer NtpWriter ntpw(kDefOptNtpFormat, gOptRunNu, gOptRanSeed); ntpw.CustomizeFilenamePrefix(gOptEvFilePrefix); diff --git a/src/Apps/gEvGenDM.cxx b/src/Apps/gEvGenDM.cxx index 6f48039bd1..9d2015cac5 100644 --- a/src/Apps/gEvGenDM.cxx +++ b/src/Apps/gEvGenDM.cxx @@ -300,14 +300,17 @@ void GenerateEventsAtFixedInitState(void) // Create init state InitialState init_state(target, dark_matter); - + +// +/* bool unitary = CheckUnitarityLimit(); if (!unitary) { LOG("gevgen_dm", pFATAL) << "Cross-section risks exceeding unitarity limit - Exiting"; exit(1); } - +*/ +// // Create/config event generation driver GEVGDriver evg_driver; diff --git a/src/Framework/Conventions/KinePhaseSpace.h b/src/Framework/Conventions/KinePhaseSpace.h index ed83495be2..a9613e829b 100644 --- a/src/Framework/Conventions/KinePhaseSpace.h +++ b/src/Framework/Conventions/KinePhaseSpace.h @@ -67,7 +67,7 @@ typedef enum EKinePhaseSpace { // TODO: rename this value when the correct variables are identified kPSTAfE, kPSEgTlOgfE, - kPSDMELEvGen, // Equivalent to kPSQELEvGen for Dark Matter scattering + kPSDMELEvGen, // Equivalent to kPSQELEvGen for Dark Matter scattering kPSxQ2fE, kPSlog10xlog10Q2fE, kPSEDNufE, // Used for Dark Neutrinos, two body final state diff --git a/src/Framework/EventGen/GMCJDriver.cxx b/src/Framework/EventGen/GMCJDriver.cxx index 4e9e2041d4..87d633705a 100644 --- a/src/Framework/EventGen/GMCJDriver.cxx +++ b/src/Framework/EventGen/GMCJDriver.cxx @@ -130,7 +130,7 @@ void GMCJDriver::SetXSecSplineNbins(int nbins) fXSecSplineNbins = nbins; LOG("GMCJDriver", pNOTICE) - << "Number of bins in energy used for the xsec splines: " + << "Number of bins in energy used for the xsec splines: " << fXSecSplineNbins; } //___________________________________________________________________________ @@ -147,7 +147,7 @@ void GMCJDriver::SetPmaxNbins(int nbins) fPmaxNbins = nbins; LOG("GMCJDriver", pNOTICE) - << "Number of bins in energy used for maximum int. probability: " + << "Number of bins in energy used for maximum int. probability: " << fPmaxNbins; } //___________________________________________________________________________ @@ -654,11 +654,11 @@ void GMCJDriver::BootstrapXSecSplineSummation(void) // knots with zero y values (although the GENIE Spline object handles it) double min = rE.min; double max = TMath::Min(rE.max, fEmax); - + // Because of edge issue (see GENIE docdb 297) these lines are commented out // double dE = fEmax/10.; // double max = (fEmax+dE < rE.max) ? fEmax+dE : rE.max; - + // in the logaritmic binning is important to have a narrow binning to // describe better the glashow resonance peak. evgdriver->CreateXSecSumSpline(fXSecSplineNbins,min,max,true); @@ -747,7 +747,7 @@ void GMCJDriver::ComputeProbScales(void) // get xsec sum over all modelled processes for given neutrino+target) double sxsecLow = evgdriver->XSecSumSpline()->Evaluate(EvLow); - double sxsecHigh = evgdriver->XSecSumSpline()->Evaluate(EvHigh); + double sxsecHigh = evgdriver->XSecSumSpline()->Evaluate(EvHigh); // get max{path-length x density} double plmax = fMaxPathLengths.PathLength(target_pdgc); @@ -860,11 +860,11 @@ EventRecord * GMCJDriver::GenerateEvent1Try(void) } if (fForceInteraction) { - bool pl_ok = this->ComputePathLengths(); - if(!pl_ok) { - LOG("GMCJDriver", pFATAL) << "** Cannot calculate path lenths!"; - exit(1); - } + bool pl_ok = this->ComputePathLengths(); + if(!pl_ok) { + LOG("GMCJDriver", pFATAL) << "** Cannot calculate path lenths!"; + exit(1); + } if(fCurPathLengths.AreAllZero()) { LOG("GMCJDriver", pNOTICE) << "** Rejecting current flux neutrino (misses generation volume)"; @@ -877,53 +877,53 @@ EventRecord * GMCJDriver::GenerateEvent1Try(void) } else { // Compute the interaction probabilities assuming max. path lengths - // and decide whether the neutrino would interact -- - // Many flux neutrinos should be rejected here, drastically reducing - // the number of neutrinos that I need to propagate through the - // actual detector geometry (this is skipped when using + // and decide whether the neutrino would interact -- + // Many flux neutrinos should be rejected here, drastically reducing + // the number of neutrinos that I need to propagate through the + // actual detector geometry (this is skipped when using // pre-calculated flux interaction probabilities) if(fPreSelect) { - LOG("GMCJDriver", pNOTICE) + LOG("GMCJDriver", pNOTICE) << "Computing interaction probabilities for max. path lengths"; Psum = this->ComputeInteractionProbabilities(true /* <- max PL*/); Pno = 1-Psum; LOG("GMCJDriver", pNOTICE) - << "The no-interaction probability (max. path lengths) is: " + << "The no-interaction probability (max. path lengths) is: " << 100*Pno << " %"; if(Pno<0.) { - LOG("GMCJDriver", pFATAL) + LOG("GMCJDriver", pFATAL) << "Negative no-interaction probability! (P = " << 100*Pno << " %)" << " Particle E=" << fFluxDriver->Momentum().E() << " type=" << fFluxDriver->PdgCode() << "Psum=" << Psum; gAbortingInErr=true; exit(1); } if(R>=1-Pno) { - LOG("GMCJDriver", pNOTICE) + LOG("GMCJDriver", pNOTICE) << "** Rejecting current flux neutrino"; return 0; } - } // preselect + } // preselect bool pl_ok = false; - // If possible use pre-generated flux neutrino interaction probabilities + // If possible use pre-generated flux neutrino interaction probabilities if(fFluxIntTree){ - Psum = this->PreGenFluxInteractionProbability(); - } + Psum = this->PreGenFluxInteractionProbability(); + } // Else compute them in the usual manner else { // Compute (pathLength x density x weight fraction) for all materials // in the input geometry, for the neutrino generated by the flux driver pl_ok = this->ComputePathLengths(); if(!pl_ok) { - LOG("GMCJDriver", pERROR) + LOG("GMCJDriver", pERROR) << "** Rejecting current flux neutrino (err computing path-lengths)"; return 0; } if(fCurPathLengths.AreAllZero()) { - LOG("GMCJDriver", pNOTICE) + LOG("GMCJDriver", pNOTICE) << "** Rejecting current flux neutrino (misses generation volume)"; return 0; } @@ -935,14 +935,14 @@ EventRecord * GMCJDriver::GenerateEvent1Try(void) LOG("GMCJDriver", pNOTICE) << "** Rejecting current flux neutrino (has null interaction probability)"; return 0; - } + } // Now decide whether the current neutrino interacts Pno = 1-Psum; LOG("GMCJDriver", pNOTICE) << "The actual 'no interaction' probability is: " << 100*Pno << " %"; if(Pno<0.) { - LOG("GMCJDriver", pFATAL) + LOG("GMCJDriver", pFATAL) << "Negative no interactin probability! (P = " << 100*Pno << " %)"; // print info about what caused the problem @@ -966,27 +966,27 @@ EventRecord * GMCJDriver::GenerateEvent1Try(void) exit(1); } if(R>=1-Pno) { - LOG("GMCJDriver", pNOTICE) + LOG("GMCJDriver", pNOTICE) << "** Rejecting current flux neutrino"; return 0; } // - // The flux neutrino interacts! + // The flux neutrino interacts! // - // Calculate path lengths for first time and check potential mismatch if + // Calculate path lengths for first time and check potential mismatch if // used pre-generated flux interaction probabilities if(fFluxIntTree){ - pl_ok = this->ComputePathLengths(); - if(!pl_ok) { - LOG("GMCJDriver", pFATAL) << "** Cannot calculate path lenths!"; - exit(1); - } + pl_ok = this->ComputePathLengths(); + if(!pl_ok) { + LOG("GMCJDriver", pFATAL) << "** Cannot calculate path lenths!"; + exit(1); + } double Psum_curr = this->ComputeInteractionProbabilities(false /* <- actual PL */); - bool mismatch = TMath::Abs(Psum-Psum_curr) > controls::kASmallNum; + bool mismatch = TMath::Abs(Psum-Psum_curr) > controls::kASmallNum; if(mismatch){ - LOG("GMCJDriver", pFATAL) << + LOG("GMCJDriver", pFATAL) << "** Mismatch between pre-calculated and current interaction "<< "probabilities!"; exit(1); diff --git a/src/Framework/Interaction/Interaction.cxx b/src/Framework/Interaction/Interaction.cxx index f3456f6676..55cb4bddf8 100644 --- a/src/Framework/Interaction/Interaction.cxx +++ b/src/Framework/Interaction/Interaction.cxx @@ -915,7 +915,7 @@ Interaction * Interaction::MECNC( Interaction * Interaction::MECEM(int tgt, int probe, double E) { - Interaction * interaction = + Interaction * interaction = Interaction::Create(tgt, probe, kScMEC, kIntEM); InitialState * init_state = interaction->InitStatePtr(); @@ -1093,6 +1093,31 @@ Interaction * Interaction::DMDI( return interaction; } //___________________________________________________________________________ +Interaction * Interaction::RESDM(int target, int hitnuc, int probe, double E) +{ + Interaction * interaction = + Interaction::Create(target,probe,kScDarkMatterResonant, kIntDarkMatter); + + InitialState * init_state = interaction->InitStatePtr(); + init_state->SetProbeE(E); + init_state->TgtPtr()->SetHitNucPdg(hitnuc); + + return interaction; +} +//___________________________________________________________________________ +Interaction * Interaction::RESDM( + int target, int hitnuc, int probe, const TLorentzVector & p4probe) +{ + Interaction * interaction = + Interaction::Create(target,probe,kScDarkMatterResonant, kIntDarkMatter); + + InitialState * init_state = interaction->InitStatePtr(); + init_state->SetProbeP4(p4probe); + init_state->TgtPtr()->SetHitNucPdg(hitnuc); + + return interaction; +} +//___________________________________________________________________________ Interaction * Interaction::HNL(int probe, double E, int decayed_mode) { Interaction * interaction = diff --git a/src/Framework/Interaction/Interaction.h b/src/Framework/Interaction/Interaction.h index 0336a21076..e7d46bc5a5 100644 --- a/src/Framework/Interaction/Interaction.h +++ b/src/Framework/Interaction/Interaction.h @@ -11,7 +11,7 @@ \author Costas Andreopoulos University of Liverpool - Changes required to implement the GENIE Boosted Dark Matter module + Changes required to implement the GENIE Boosted Dark Matter module were installed by Josh Berger (Univ. of Wisconsin) \created April 25, 2004 @@ -48,16 +48,16 @@ const UInt_t kISkipProcessChk = 1<<17; ///< if set, skip process validity c const UInt_t kISkipKinematicChk = 1<<16; ///< if set, skip kinematic validity checks const UInt_t kIAssumeFreeNucleon = 1<<15; ///< const UInt_t kIAssumeFreeElectron = 1<<15; ///< -const UInt_t kINoNuclearCorrection = 1<<14; ///< if set, inhibit nuclear corrections +const UInt_t kINoNuclearCorrection = 1<<14; ///< if set, inhibit nuclear corrections class Interaction; -ostream & operator << (ostream & stream, const Interaction & i); +ostream & operator << (ostream & stream, const Interaction & i); class Interaction : public TObject { public: using TObject::Print; // suppress clang 'hides overloaded virtual function [-Woverloaded-virtual]' warnings - using TObject::Copy; // + using TObject::Copy; // Interaction(); Interaction(const InitialState & init, const ProcessInfo & proc); @@ -88,7 +88,7 @@ class Interaction : public TObject { int FSPrimLeptonPdg (void) const; ///< final state primary lepton pdg int RecoilNucleonPdg (void) const; ///< recoil nucleon pdg TParticlePDG * FSPrimLepton (void) const; ///< final state primary lepton - TParticlePDG * RecoilNucleon (void) const; ///< recoil nucleon + TParticlePDG * RecoilNucleon (void) const; ///< recoil nucleon // Copy, reset, print itself and build string code void Reset (void); @@ -132,10 +132,10 @@ class Interaction : public TObject { static Interaction * DFRCC (int tgt, int nuc, int probe, double E=0); static Interaction * DFRCC (int tgt, int nuc, int probe, const TLorentzVector & p4probe); static Interaction * COHCC (int tgt, int probe, unsigned int prod_pdg, double E=0); - static Interaction * COHCC (int tgt, int probe, unsigned int prod_pdg, + static Interaction * COHCC (int tgt, int probe, unsigned int prod_pdg, const TLorentzVector & p4probe); static Interaction * COHNC (int tgt, int probe, unsigned int prod_pdg, double E=0); - static Interaction * COHNC (int tgt, int probe, unsigned int prod_pdg, + static Interaction * COHNC (int tgt, int probe, unsigned int prod_pdg, const TLorentzVector & p4probe); static Interaction * CEvNS (int tgt, int probe, double E=0); static Interaction * CEvNS (int tgt, int probe, const TLorentzVector & p4probe); @@ -164,6 +164,8 @@ class Interaction : public TObject { static Interaction * DMDI (int tgt, int nuc, int qrk, bool sea, int probe, double E=0); static Interaction * DMDI (int tgt, int nuc, int probe, const TLorentzVector & p4probe); static Interaction * DMDI (int tgt, int nuc, int qrk, bool sea, int probe, const TLorentzVector & p4probe); + static Interaction * RESDM (int tgt, int nuc, int probe, double E=0); + static Interaction * RESDM (int tgt, int nuc, int probe, const TLorentzVector & p4probe); static Interaction * HNL (int probe, double E=0, int decayed_mode=-1); private: @@ -181,7 +183,7 @@ class Interaction : public TObject { Kinematics * fKinematics; ///< kinematical variables XclsTag * fExclusiveTag; ///< Additional info for exclusive channels KPhaseSpace * fKinePhSp; ///< Kinematic phase space - + ClassDef(Interaction,2) }; diff --git a/src/Framework/Interaction/KPhaseSpace.cxx b/src/Framework/Interaction/KPhaseSpace.cxx index 06aa68e41d..5535bc732f 100644 --- a/src/Framework/Interaction/KPhaseSpace.cxx +++ b/src/Framework/Interaction/KPhaseSpace.cxx @@ -8,6 +8,7 @@ Changes required to implement the GENIE Boosted Dark Matter module were installed by Josh Berger (Univ. of Wisconsin) + and Zachary Orr (Colorado State Univ.) */ //____________________________________________________________________________ @@ -133,6 +134,7 @@ double KPhaseSpace::Threshold(void) const pi.IsDarkMatterElastic() || pi.IsInverseBetaDecay() || pi.IsResonant() || + pi.IsDarkMatterResonant() || pi.IsDeepInelastic() || pi.IsDarkMatterDeepInelastic() || pi.IsDiffractive()) @@ -168,7 +170,7 @@ double KPhaseSpace::Threshold(void) const double smin = TMath::Power(Wmin+ml,2.); double Ethr = 0.5*(smin-Mn2)/Mn; // threshold is different for dark matter case - if(pi.IsDarkMatterElastic() || pi.IsDarkMatterDeepInelastic()) { + if(pi.IsDarkMatterElastic() || pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant()) { // Correction to minimum kinematic variables Wmin = Mn; smin = TMath::Power(Wmin+ml,2); @@ -291,6 +293,7 @@ bool KPhaseSpace::IsAboveThreshold(void) const pi.IsDarkMatterElastic() || pi.IsInverseBetaDecay() || pi.IsResonant() || + pi.IsDarkMatterResonant() || pi.IsDeepInelastic() || pi.IsDarkMatterDeepInelastic() || pi.IsDiffractive() || @@ -330,7 +333,7 @@ bool KPhaseSpace::IsAllowed(void) const // RES // Check the running W vs the W limits // & the running Q2 vs Q2 limits for the given W - if(pi.IsResonant()) { + if(pi.IsResonant() || pi.IsDarkMatterResonant()) { Range1D_t Wl = this->WLim(); Range1D_t Q2l = this->Q2Lim_W(); double W = kine.W(); @@ -438,7 +441,7 @@ Range1D_t KPhaseSpace::WLim(void) const bool is_em = pi.IsEM(); bool is_qel = pi.IsQuasiElastic() || pi.IsInverseBetaDecay() || pi.IsDarkMatterElastic(); bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive(); - bool is_dmdis = pi.IsDarkMatterDeepInelastic(); + bool is_dminel = pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant(); if(is_qel) { double MR = fInteraction->RecoilNucleon()->Mass(); @@ -469,7 +472,7 @@ Range1D_t KPhaseSpace::WLim(void) const return Wl; } - if(is_dmdis) { + if(is_dminel) { const InitialState & init_state = fInteraction->InitState(); double Ev = init_state.ProbeE(kRfHitNucRest); double M = init_state.Tgt().HitNucP4Ptr()->M(); //can be off m/shell @@ -514,9 +517,9 @@ Range1D_t KPhaseSpace::Q2Lim_W(void) const bool is_inel = pi.IsDeepInelastic() || pi.IsResonant() || pi.IsDiffractive(); bool is_coh = pi.IsCoherentProduction(); bool is_dme = pi.IsDarkMatterElastic(); - bool is_dmdis = pi.IsDarkMatterDeepInelastic(); + bool is_dminel = pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant(); - if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dmdis) return Q2l; + if(!is_qel && !is_inel && !is_coh && !is_dme && !is_dminel) return Q2l; if(is_coh) { return Q2Lim(); @@ -533,7 +536,7 @@ Range1D_t KPhaseSpace::Q2Lim_W(void) const if (pi.IsInverseBetaDecay()) { Q2l = kinematics::InelQ2Lim_W(Ev,M,ml,W, controls::kMinQ2Limit_VLE); - } else if (is_dme || is_dmdis) { + } else if (is_dme || is_dminel) { Q2l = kinematics::DarkQ2Lim_W(Ev,M,ml,W); } else { Q2l = is_em ? kinematics::electromagnetic::InelQ2Lim_W(Ev,ml,M,W) : kinematics::InelQ2Lim_W(Ev,M,ml,W); @@ -571,9 +574,9 @@ Range1D_t KPhaseSpace::Q2Lim(void) const bool is_coh = pi.IsCoherentProduction(); bool is_cevns = pi.IsCoherentElastic(); bool is_dme = pi.IsDarkMatterElastic(); - bool is_dmdis = pi.IsDarkMatterDeepInelastic(); + bool is_dminel = pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant(); - if(!is_qel && !is_inel && !is_coh && !is_cevns && !is_dme && !is_dmdis) return Q2l; + if(!is_qel && !is_inel && !is_coh && !is_cevns && !is_dme && !is_dminel) return Q2l; const InitialState & init_state = fInteraction->InitState(); double Ev = init_state.ProbeE(kRfHitNucRest); @@ -652,7 +655,7 @@ Range1D_t KPhaseSpace::Q2Lim(void) const return Q2l; } - if (is_dmdis) { + if (is_dminel) { Q2l = kinematics::DarkQ2Lim(Ev,M,ml); return Q2l; } @@ -695,8 +698,8 @@ Range1D_t KPhaseSpace::XLim(void) const return xl; } //DMDIS - bool is_dmdis = pi.IsDarkMatterDeepInelastic(); - if(is_dmdis) { + bool is_dminel = pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant(); + if(is_dminel) { const InitialState & init_state = fInteraction->InitState(); double Ev = init_state.ProbeE(kRfHitNucRest); double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell @@ -747,8 +750,8 @@ Range1D_t KPhaseSpace::YLim(void) const return yl; } //DMDIS - bool is_dmdis = pi.IsDarkMatterDeepInelastic(); - if(is_dmdis) { + bool is_dminel = pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant(); + if(is_dminel) { const InitialState & init_state = fInteraction->InitState(); double Ev = init_state.ProbeE(kRfHitNucRest); double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell @@ -820,8 +823,8 @@ Range1D_t KPhaseSpace::YLim_X(void) const return yl; } //DMDIS - bool is_dmdis = pi.IsDarkMatterDeepInelastic(); - if(is_dmdis) { + bool is_dminel = pi.IsDarkMatterDeepInelastic() || pi.IsDarkMatterResonant(); + if(is_dminel) { const InitialState & init_state = fInteraction->InitState(); double Ev = init_state.ProbeE(kRfHitNucRest); double M = init_state.Tgt().HitNucP4Ptr()->M(); // can be off m/shell @@ -1006,7 +1009,7 @@ double KPhaseSpace::Threshold_SPP_iso(bool isMassless) const { const InitialState & init_state = fInteraction->InitState(); PDGLibrary * pdglib = PDGLibrary::Instance(); - + // imply isospin symmetry double mpi = (pdglib->Find(kPdgPiP)->Mass() + pdglib->Find(kPdgPi0)->Mass() + pdglib->Find(kPdgPiM)->Mass())/3; double M = (pdglib->Find(kPdgProton)->Mass() + pdglib->Find(kPdgNeutron)->Mass())/2; @@ -1055,7 +1058,7 @@ Range1D_t KPhaseSpace::WLim_SPP(bool isMassless) const Wl.min *= 1. + std::numeric_limits::epsilon(); Wl.max *= 1. - std::numeric_limits::epsilon(); } - + return Wl; } //____________________________________________________________________________ @@ -1091,7 +1094,7 @@ Range1D_t KPhaseSpace::WLim_SPP_iso(bool isMassless) const Wl.min *= 1. + std::numeric_limits::epsilon(); Wl.max *= 1. - std::numeric_limits::epsilon(); } - + return Wl; } //____________________________________________________________________________ @@ -1182,8 +1185,7 @@ Range1D_t KPhaseSpace::Q2Lim_W_SPP_iso(bool isMassless) const Q2l.min *= 1. + std::numeric_limits::epsilon(); Q2l.max *= 1. - std::numeric_limits::epsilon(); } - + return Q2l; } //____________________________________________________________________________ - diff --git a/src/Framework/Interaction/ProcessInfo.cxx b/src/Framework/Interaction/ProcessInfo.cxx index 1456b9c383..01aba4c372 100644 --- a/src/Framework/Interaction/ProcessInfo.cxx +++ b/src/Framework/Interaction/ProcessInfo.cxx @@ -8,6 +8,7 @@ Changes required to implement the GENIE Boosted Dark Matter module were installed by Josh Berger (Univ. of Wisconsin) + and Zachary Orr (Colorado State Univ.) Changes required to implement the GENIE Dark Neutrino module were installed by Iker de Icaza (Univ. of Sussex) @@ -101,6 +102,11 @@ bool ProcessInfo::IsResonant(void) const return (fScatteringType == kScResonant); } //____________________________________________________________________________ +bool ProcessInfo::IsDarkMatterResonant(void) const +{ + return (fScatteringType == kScDarkMatterResonant); +} +//____________________________________________________________________________ bool ProcessInfo::IsCoherentProduction(void) const { return (fScatteringType == kScCoherentProduction); diff --git a/src/Framework/Interaction/ProcessInfo.h b/src/Framework/Interaction/ProcessInfo.h index cf2548089c..69d2cbfe49 100644 --- a/src/Framework/Interaction/ProcessInfo.h +++ b/src/Framework/Interaction/ProcessInfo.h @@ -13,6 +13,7 @@ Changes required to implement the GENIE Boosted Dark Matter module were installed by Josh Berger (Univ. of Wisconsin) + and Zachary Orr (Colorado State Univ.) Changes required to implement the GENIE Dark Neutrino module were installed by Iker de Icaza (Univ. of Sussex) @@ -64,6 +65,7 @@ class ProcessInfo : public TObject { bool IsDeepInelastic (void) const; bool IsDarkMatterDeepInelastic (void) const; bool IsResonant (void) const; + bool IsDarkMatterResonant (void) const; bool IsCoherentProduction (void) const; bool IsCoherentElastic (void) const; bool IsSinglePion (void) const; diff --git a/src/Framework/Interaction/ScatteringType.h b/src/Framework/Interaction/ScatteringType.h index bd4f80b231..9acfab8e73 100644 --- a/src/Framework/Interaction/ScatteringType.h +++ b/src/Framework/Interaction/ScatteringType.h @@ -34,7 +34,7 @@ namespace genie { // append to the end of the total list and set a new enum counter value. typedef enum EScatteringType { - kScUnknown = -100, + kScUnknown = -100, kScNull = 0, kScQuasiElastic, kScSingleKaon, @@ -56,7 +56,8 @@ typedef enum EScatteringType { kScDarkMatterElastic = 101, kScDarkMatterDeepInelastic, kScDarkMatterElectron, - kScNorm + kScDarkMatterResonant, + kScNorm = 110 } ScatteringType_t; class ScatteringType @@ -89,6 +90,7 @@ class ScatteringType case(kScDarkMatterElastic) : return "DMEL"; break; case(kScDarkMatterDeepInelastic) : return "DMDIS"; break; case(kScDarkMatterElectron) : return "DME"; break; + case(kScDarkMatterResonant) : return "DMRES"; break; default : return "Unknown"; break; } return "Unknown"; diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.cxx b/src/Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.cxx new file mode 100644 index 0000000000..8441a688a9 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.cxx @@ -0,0 +1,236 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Zachary W. Orr + Colorado State University + + based on code by + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory +*/ +//____________________________________________________________________________ + +// #include +// #if ROOT_VERSION_CODE >= ROOT_VERSION(5,15,6) +// #include +// #else +// #include +// #endif + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/GHEP/GHepStatus.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Interaction/SppChannel.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/BaryonResonance.h" +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/Utils/PrintUtils.h" +// #include "Physics/Decay/DecayModelI.h" +#include "Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.h" + +using namespace genie; +using namespace genie::constants; + +//___________________________________________________________________________ +DMRESHadronicSystemGenerator::DMRESHadronicSystemGenerator() : +HadronicSystemGenerator("genie::DMRESHadronicSystemGenerator") +{ + +} +//___________________________________________________________________________ +DMRESHadronicSystemGenerator::DMRESHadronicSystemGenerator(string config): +HadronicSystemGenerator("genie::DMRESHadronicSystemGenerator", config) +{ + +} +//___________________________________________________________________________ +DMRESHadronicSystemGenerator::~DMRESHadronicSystemGenerator() +{ + +} +//___________________________________________________________________________ +void DMRESHadronicSystemGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ +// This method generates the final state hadronic system + + // Get the right resonance PDG code so that the selected resonance + // conserves charge + int pdgc = GetResonancePdgCode(evrec); + + // Add the selected resonance + this->AddResonance(evrec,pdgc); + + // Decay the resonance (and its decay products, if they include resonances) + fResonanceDecayer->ProcessEventRecord(evrec); + + // Add the baryon resonance decay products at the event record + //this->AddResonanceDecayProducts(evrec,pdgc); + + // Handle resonance decay channels to other resonances or short-living + // partices + //LOG("RESHadronicVtx", pNOTICE) + // << "Decay any resonance in the initial resonance decay products"; + //this->PreHadronTransportDecays(evrec); +} +//___________________________________________________________________________ +int DMRESHadronicSystemGenerator::GetResonancePdgCode(GHepRecord * evrec) const +{ +// In the RES thread the resonance is specifed when selecting interaction +// This method adds it to the GHEP record. + + Interaction * interaction = evrec->Summary(); + + // Get resonance id + const XclsTag & xcls = interaction->ExclTag(); + assert(xcls.KnownResonance()); + Resonance_t res = xcls.Resonance(); + + // Get resonance charge + int q_res = this->ResonanceCharge(evrec); + + // Find resonance PDG code from resonance charge and id + int pdgc = utils::res::PdgCode(res, q_res); + + LOG("RESHadronicVtx", pNOTICE) + << "Selected event has RES with PDGC = " << pdgc << ", Q = " << q_res; + + return pdgc; +} +//___________________________________________________________________________ +void DMRESHadronicSystemGenerator::AddResonance( + GHepRecord * evrec, int pdgc) const +{ + //-- Compute RES p4 = p4(neutrino) + p4(hit nucleon) - p4(primary lepton) + TLorentzVector p4 = this->Hadronic4pLAB(evrec); + + //-- Add the resonance at the EventRecord + GHepStatus_t ist = kIStPreDecayResonantState; + int mom = evrec->HitNucleonPosition(); + + //-- Get vtx position + GHepParticle * neutrino = evrec->Probe(); + const TLorentzVector & vtx = *(neutrino->X4()); + + evrec->AddParticle(pdgc, ist, mom,-1,-1,-1, p4, vtx); +} +//___________________________________________________________________________ +// void RESHadronicSystemGenerator::AddResonanceDecayProducts( +// GHepRecord * evrec, int pdgc) const +// { +// // Decay the baryon resonance, take the decay products, boost them in the LAB +// // and add them in the GHEP record. +// // Unlike the SPP thread where the resonance decay products are determined +// // from the selected SPP channel, in the RES thread we can any of the the +// // resonance's kinematically available(the RES is not on the mass shell)decay +// // channels +// +// // find the resonance position +// int irpos = evrec->ParticlePosition(pdgc, kIStPreDecayResonantState, 0); +// assert(irpos>0); +// +// // access the GHEP entry +// GHepParticle * resonance = evrec->Particle(irpos); +// assert(resonance); +// +// // resonance location +// const TLorentzVector & x4 = *(resonance->X4()); +// +// // prepare the decayer inputs +// DecayerInputs_t dinp; +// dinp.PdgCode = pdgc; +// dinp.P4 = resonance->P4(); +// +// // do the decay +// TClonesArray * decay_products = fResonanceDecayer->Decay(dinp); +// if(!decay_products) { +// LOG("RESHadronicVtx", pWARN) << "Got an empty decay product list!"; +// LOG("RESHadronicVtx", pWARN) +// << "Quitting the current event generation thread"; +// +// evrec->EventFlags()->SetBitNumber(kHadroSysGenErr, true); +// +// genie::exceptions::EVGThreadException exception; +// exception.SetReason("Not enough phase space for hadronizer"); +// exception.SwitchOnFastForward(); +// throw exception; +// +// return; +// } +// +// // get the decay weight (if any) +// double wght = fResonanceDecayer->Weight(); +// +// // update the event weight +// evrec->SetWeight(wght * evrec->Weight()); +// +// // decide the istatus of decay products +// GHepParticle * nuc = evrec->TargetNucleus(); +// GHepStatus_t dpist = (nuc) ? kIStHadronInTheNucleus : kIStStableFinalState; +// +// // if the list is not empty, boost and copy the decay products in GHEP +// if(decay_products) { +// +// // first, mark the resonance as decayed +// resonance->SetStatus(kIStDecayedState); +// +// // loop over the daughter and add them to the event record +// TMCParticle * dpmc = 0; +// TObjArrayIter decay_iter(decay_products); +// +// while( (dpmc = (TMCParticle *) decay_iter.Next()) ) { +// +// int dppdg = dpmc->GetKF(); +// double px = dpmc->GetPx(); +// double py = dpmc->GetPy(); +// double pz = dpmc->GetPz(); +// double E = dpmc->GetEnergy(); +// TLorentzVector p4(px,py,pz,E); +// +// //-- Only add the decay products - the mother particle already exists +// if(dpmc->GetKS()==1) { +// evrec->AddParticle(dppdg,dpist,irpos,-1,-1,-1, p4, x4); +// } +// } +// +// // done, release the original list +// decay_products->Delete(); +// delete decay_products; +// }// !=0 +// } +//___________________________________________________________________________ +void DMRESHadronicSystemGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//___________________________________________________________________________ +void DMRESHadronicSystemGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//___________________________________________________________________________ +void DMRESHadronicSystemGenerator::LoadConfig(void) +{ + fResonanceDecayer = 0; + //fPreINukeDecayer = 0; + + // Get the specified decayers + fResonanceDecayer = + dynamic_cast (this->SubAlg("Decayer")); + assert(fResonanceDecayer); + // fPreINukeDecayer = + // dynamic_cast (this->SubAlg("PreTransportDecayer")); + // assert(fPreINukeDecayer); +} +//___________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.h b/src/Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.h new file mode 100644 index 0000000000..960a8e7188 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESHadronicSystemGenerator.h @@ -0,0 +1,67 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESHadronicSystemGenerator + +\brief Generates the 'final state' hadronic system in DM RES interactions. (same as v RES) + It adds the remnant nucleus (if any), the pre-selected resonance + and the resonance decay products at the GHEP record. + Unlike the SPP thread, in the RES thread the resonance is specified + at the time an interaction is selected but its decay products not + (semi-inclusive resonance reactions). The off the mass-shell baryon + resonance is decayed using a phase space generator. All kinematically + available decay channels are being used (not just 1 pi channels). + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Zachary W. Orr + Colorado State University + + based on code by + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + + +\created November 23, 2004 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DMRES_HADRONIC_SYSTEM_GENERATOR_H_ +#define _DMRES_HADRONIC_SYSTEM_GENERATOR_H_ + +#include "Physics/Common/HadronicSystemGenerator.h" + +namespace genie { + +class Decayer; + +class DMRESHadronicSystemGenerator : public HadronicSystemGenerator { + +public : + DMRESHadronicSystemGenerator(); + DMRESHadronicSystemGenerator(string config); + ~DMRESHadronicSystemGenerator(); + + // implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) 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); + int GetResonancePdgCode (GHepRecord * evrec) const; + void AddResonance (GHepRecord * evrec, int pdgc) const; + // void AddResonanceDecayProducts (GHepRecord * evrec, int pdgc) const; + + const EventRecordVisitorI * fResonanceDecayer; +}; + +} // genie namespace + +#endif // _DMRES_HADRONIC_SYSTEM_GENERATOR_H_ diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.cxx b/src/Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.cxx new file mode 100644 index 0000000000..f71d47c700 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.cxx @@ -0,0 +1,169 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Zachary W. Orr + Colorado State University + + based on code by + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory +*/ +//____________________________________________________________________________ + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/EventGen/InteractionList.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.h" + +using namespace genie; + +//___________________________________________________________________________ +DMRESInteractionListGenerator::DMRESInteractionListGenerator() : +InteractionListGeneratorI("genie::DMRESInteractionListGenerator") +{ + +} +//___________________________________________________________________________ +DMRESInteractionListGenerator::DMRESInteractionListGenerator(string config) : +InteractionListGeneratorI("genie::DMRESInteractionListGenerator", config) +{ + +} +//___________________________________________________________________________ +DMRESInteractionListGenerator::~DMRESInteractionListGenerator() +{ + +} +//___________________________________________________________________________ +InteractionList * DMRESInteractionListGenerator::CreateInteractionList( + const InitialState & init_state) const +{ + LOG("IntLst", pINFO) << "InitialState = " << init_state.AsString(); + + // In the thread generating interactions from the list produced here (RES), + // we simulate (for free and nuclear targets) semi-inclusive resonance + // interactions: v + N -> v(l) + R -> v(l) + X + // Specifically, the RES thread generates: + // + // CC: + // nu + p (A) -> l- R (A), for all resonances with Q=+2 + // nu + n (A) -> l- R (A), for all resonances with Q=+1 + // \bar{nu} + p (A) -> l+ R (A), for all resonances with Q= 0 + // \bar{nu} + n (A) -> l+ R (A), for all resonances with Q=-1 + // NC: + // nu + p (A) -> nu R (A), for all resonances with Q=+1 + // nu + n (A) -> nu R (A), for all resonances with Q= 0 + // \bar{nu} + p (A) -> \bar{nu} R (A), for all resonances with Q=+1 + // \bar{nu} + n (A) -> \bar{nu} R (A), for all resonances with Q= 0 + // + // and then the resonance R should be allowed to decay to get the full + // hadronic final state X. All decay channels kinematically accessible + // to the (off the mass-shell produced) resonance can be allowed. + + // specify the requested interaction type + InteractionType_t inttype; + // Only accept DM scattering here + if (fIsDM) inttype = kIntDarkMatter; + else { + LOG("IntLst", pWARN) + << "Unknown InteractionType! Returning NULL InteractionList " + << "for init-state: " << init_state.AsString(); + return 0; + } + + // create a process information object + ProcessInfo proc_info(kScDarkMatterResonant, inttype); + + // learn whether the input nuclear or free target has avail. p and n + const Target & inp_target = init_state.Tgt(); + bool hasP = (inp_target.Z() > 0); + bool hasN = (inp_target.N() > 0); + + // possible hit nucleons + const int hit_nucleon[2] = {kPdgProton, kPdgNeutron}; + + // create an interaction list + InteractionList * intlist = new InteractionList; + + // loop over all baryon resonances considered in current MC job + unsigned int nres = fResList.NResonances(); + for(unsigned int ires = 0; ires < nres; ires++) { + + //get current resonance + Resonance_t res = fResList.ResonanceId(ires); + + // loop over hit nucleons + for(int i=0; i<2; i++) { + + // proceed only if the hit nucleon exists in the current init state + if(hit_nucleon[i]==kPdgProton && !hasP) continue; + if(hit_nucleon[i]==kPdgNeutron && !hasN) continue; + + // proceed only if the current resonance conserves charge + // (the only problematic case is when the RES charge has to be +2 + // because then only Delta resonances are possible) +// bool skip_res = proc_info.IsWeakCC() && +// pdg::IsNeutrino(init_state.ProbePdg()) && +// (hit_nucleon[i]==kPdgProton) && +// (!utils::res::IsDelta(res)); +// if(skip_res) continue; + + // create an interaction + Interaction * interaction = new Interaction(init_state, proc_info); + + // add the struck nucleon + Target * target = interaction->InitStatePtr()->TgtPtr(); + target->SetHitNucPdg(hit_nucleon[i]); + + // add the baryon resonance in the exclusive tag + XclsTag * xcls = interaction->ExclTagPtr(); + xcls->SetResonance(res); + + // add the interaction at the interaction list + intlist->push_back(interaction); + + }//hit nucleons + } //resonances + + if(intlist->size() == 0) { + LOG("IntLst", pERROR) + << "Returning NULL InteractionList for init-state: " + << init_state.AsString(); + delete intlist; + return 0; + } + + return intlist; +} +//___________________________________________________________________________ +void DMRESInteractionListGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfigData(); +} +//____________________________________________________________________________ +void DMRESInteractionListGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfigData(); +} +//____________________________________________________________________________ +void DMRESInteractionListGenerator::LoadConfigData(void) +{ + string resonances = ""; + this->GetParam("ResonanceNameList", resonances); + SLOG("IntLst", pDEBUG) << "Resonance list: " << resonances; + + fResList.Clear(); + fResList.DecodeFromNameList(resonances); + LOG("IntLst", pDEBUG) << fResList; + + this->GetParamDef("is-DM", fIsDM, false); +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.h b/src/Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.h new file mode 100644 index 0000000000..558ab01cd2 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESInteractionListGenerator.h @@ -0,0 +1,54 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESInteractionListGenerator + +\brief Creates a list of all the interactions that can be generated by the + RES thread (generates semi-inclusive resonance reactions). + Concrete implementations of the InteractionListGeneratorI interface. + follows same procedure as RESInteractionListGenerator.h + +\author Zach Orr + Colorado State University + + +\created Jan 27, 2024 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DMRES_INTERACTION_LIST_GENERATOR_H_ +#define _DMRES_INTERACTION_LIST_GENERATOR_H_ + +#include "Framework/ParticleData/BaryonResList.h" +#include "Framework/EventGen/InteractionListGeneratorI.h" + +namespace genie { + +class DMRESInteractionListGenerator : public InteractionListGeneratorI { + +public : + DMRESInteractionListGenerator(); + DMRESInteractionListGenerator(string config); + ~DMRESInteractionListGenerator(); + + // 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 fIsDM; + BaryonResList fResList; +}; + +} // genie namespace +#endif // _DMRES_INTERACTION_LIST_GENERATOR_H_ diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESKinematicsGenerator.cxx b/src/Physics/BoostedDarkMatter/EventGen/DMRESKinematicsGenerator.cxx new file mode 100644 index 0000000000..bac01c3af7 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESKinematicsGenerator.cxx @@ -0,0 +1,372 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + /author Zach Orr + Colorado State University + follows same procedure as RESKinematicsGenerator.cxx + + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory +*/ +//____________________________________________________________________________ + +#include +#include +#include + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/KineVar.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/Numerical/MathUtils.h" +#include "Framework/ParticleData/BaryonResonance.h" +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Utils/KineUtils.h" +#include "Physics/BoostedDarkMatter/EventGen/DMRESKinematicsGenerator.h" + +using namespace genie; +using namespace genie::controls; +using namespace genie::utils; + +//___________________________________________________________________________ +DMRESKinematicsGenerator::DMRESKinematicsGenerator() : +KineGeneratorWithCache("genie::DMRESKinematicsGenerator") +{ + fEnvelope = 0; +} +//___________________________________________________________________________ +DMRESKinematicsGenerator::DMRESKinematicsGenerator(string config) : +KineGeneratorWithCache("genie::DMRESKinematicsGenerator", config) +{ + fEnvelope = 0; +} +//___________________________________________________________________________ +DMRESKinematicsGenerator::~DMRESKinematicsGenerator() +{ + if(fEnvelope) delete fEnvelope; +} +//___________________________________________________________________________ +void DMRESKinematicsGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ + if(fGenerateUniformly) { + LOG("DMRESKinematics", 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(); + + //-- Get the interaction from the GHEP record + Interaction * interaction = evrec->Summary(); + interaction->SetBit(kISkipProcessChk); + + //-- DM process + bool is_DM = interaction->ProcInfo().IsDarkMatterResonant(); + + //-- Compute the W limits + // (the physically allowed W's, unless an external cut is imposed) + const KPhaseSpace & kps = interaction->PhaseSpace(); + Range1D_t W = kps.Limits(kKVW); + + if(W.max <=0 || W.min>=W.max) { + LOG("DMRESKinematics", pWARN) << "No available phase space"; + evrec->EventFlags()->SetBitNumber(kKineGenErr, true); + genie::exceptions::EVGThreadException exception; + exception.SetReason("No available phase space"); + exception.SwitchOnFastForward(); + throw exception; + } + + const InitialState & init_state = interaction -> InitState(); + double E = init_state.ProbeE(kRfHitNucRest); + + //-- 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 = (fGenerateUniformly) ? -1 : this->MaxXSec(evrec); + + //-- Try to select a valid W, Q2 pair using the rejection method + double dW = W.max - W.min; + double xsec = -1; + + unsigned int iter = 0; + bool accept = false; + while(1) { + iter++; + if(iter > kRjMaxIterations) { + LOG("DMRESKinematics", pWARN) + << "*** Could not select a valid (W,Q^2) pair after " + << iter << " iterations"; + evrec->EventFlags()->SetBitNumber(kKineGenErr, true); + genie::exceptions::EVGThreadException exception; + exception.SetReason("Couldn't select kinematics"); + exception.SwitchOnFastForward(); + throw exception; + } + + double gW = 0; // current hadronic invariant mass + double gQ2 = 0; // current momentum transfer + double gQD2 = 0; // tranformed Q2 to take out dipole form + + if(fGenerateUniformly) + { + //-- Generate a W uniformly in the kinematically allowed range. + // For the generated W, compute the Q2 range and generate a value + // uniformly over that range + gW = W.min + dW * rnd->RndKine().Rndm(); + interaction->KinePtr()->SetW(gW); + Range1D_t Q2 = kps.Q2Lim_W(); + if(Q2.max<=0. || Q2.min>=Q2.max) continue; + gQ2 = Q2.min + (Q2.max-Q2.min) * rnd->RndKine().Rndm(); + interaction->SetBit(kISkipKinematicChk); + } + else + { + + // dm scattering + // Selecting unweighted event kinematics using an importance sampling + // method. Q2 with be transformed to QD2 to take out the dipole form. + interaction->KinePtr()->SetW(W.min); + Range1D_t Q2 = kps.Q2Lim_W(); + double Q2min = -99.; + if (is_DM) + Q2min = Q2.min + kASmallNum; + else + Q2min = 0 + kASmallNum; + double Q2max = Q2.max - kASmallNum; + + // In unweighted mode - use transform that takes out the dipole form + double QD2min = utils::kinematics::Q2toQD2(Q2max); + double QD2max = utils::kinematics::Q2toQD2(Q2min); + + gW = W.min + dW * rnd->RndKine().Rndm(); + gQD2 = QD2min + (QD2max - QD2min) * rnd->RndKine().Rndm(); + + // QD2 -> Q2 + gQ2 = utils::kinematics::QD2toQ2(gQD2); + } // uniformly over phase space? + + LOG("DMRESKinematics", pINFO) << "Trying: W = " << gW << ", Q2 = " << gQ2; + + //-- Set kinematics for current trial + interaction->KinePtr()->SetW(gW); + interaction->KinePtr()->SetQ2(gQ2); + + //-- Computing cross section for the current kinematics + xsec = fXSecModel->XSec(interaction, kPSWQD2fE); + //-- Decide whether to accept the current kinematics + if(!fGenerateUniformly) + { + // unified neutrino / electron scattering + double t = xsec_max * rnd->RndKine().Rndm(); + this->AssertXSecLimits(interaction, xsec, xsec_max); + accept = (t < xsec); + } // charged lepton or neutrino scattering? + else + { + accept = (xsec>0); + } // uniformly over phase space + + //-- If the generated kinematics are accepted, finish-up module's job + if(accept) { + LOG("DMRESKinematics", pINFO) + << "Selected: W = " << gW << ", Q2 = " << gQ2; + // reset 'trust' bits + interaction->ResetBit(kISkipProcessChk); + interaction->ResetBit(kISkipKinematicChk); + + // compute x,y for selected W,Q2 + // note: hit nucleon can be off the mass-shell + double gx=-1, gy=-1; + double M = init_state.Tgt().HitNucP4().M(); + kinematics::WQ2toXY(E,M,gW,gQ2,gx,gy); + + // set the cross section for the selected kinematics + evrec->SetDiffXSec(xsec,kPSWQ2fE); + + // for uniform kinematics, compute an event weight as + // wght = (phase space volume)*(differential xsec)/(event total xsec) + if(fGenerateUniformly) { + double vol = kinematics::PhaseSpaceVolume(interaction,kPSWQ2fE); + double totxsec = evrec->XSec(); + double wght = (vol/totxsec)*xsec; + LOG("DMRESKinematics", pNOTICE) << "Kinematics wght = "<< wght; + + // apply computed weight to the current event weight + wght *= evrec->Weight(); + LOG("DMRESKinematics", pNOTICE) << "Current event wght = " << wght; + evrec->SetWeight(wght); + } + + // lock selected kinematics & clear running values + interaction->KinePtr()->SetQ2(gQ2, true); + interaction->KinePtr()->SetW (gW, true); + interaction->KinePtr()->Setx (gx, true); + interaction->KinePtr()->Sety (gy, true); + interaction->KinePtr()->ClearRunningValues(); + + return; + } // accept + } // iterations +} +//___________________________________________________________________________ +void DMRESKinematicsGenerator::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void DMRESKinematicsGenerator::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void DMRESKinematicsGenerator::LoadConfig(void) +{ + // Safety factor for the maximum differential cross section + this->GetParamDef("MaxXSec-SafetyFactor", fSafetyFactor, 1.25); + + // Minimum energy for which max xsec would be cached, forcing explicit + // calculation for lower eneries + this->GetParamDef("Cache-MinEnergy", fEMin, 0.5); + + // Load Wcut used in DIS/RES join scheme + this->GetParam("Wcut", fWcut); + + // Maximum allowed fractional cross section deviation from maxim cross + // section used in rejection method + this->GetParamDef("MaxXSec-DiffTolerance", fMaxXSecDiffTolerance, 999999.); + assert(fMaxXSecDiffTolerance>=0); + + // Generate kinematics uniformly over allowed phase space and compute + // an event weight? + this->GetParamDef("UniformOverPhaseSpace", fGenerateUniformly, false); + + // Envelope employed when importance sampling is used + // (initialize with dummy range) + if(fEnvelope) delete fEnvelope; + fEnvelope = new TF2("res-envelope", + kinematics::RESImportanceSamplingEnvelope,0.01,1,0.01,1,4); + // stop ROOT from deleting this object of its own volition + gROOT->GetListOfFunctions()->Remove(fEnvelope); +} + +double DMRESKinematicsGenerator::ComputeMaxXSec( + const Interaction * interaction) const +{ +// Computes the maximum differential cross section in the requested phase +// space. This method overloads KineGeneratorWithCache::ComputeMaxXSec +// method and the value is cached at a circular cache branch for retrieval +// during subsequent event generation. +// The computed max differential cross section does not need to be the exact +// maximum. The number used in the rejection method will be scaled up by a +// safety factor. But this needs to be fast - do not use a very fine grid. + + double max_xsec = 0.; + + const InitialState & init_state = interaction -> InitState(); + double E = init_state.ProbeE(kRfHitNucRest); + bool is_em = interaction->ProcInfo().IsEM(); + double Q2Thres = is_em ? utils::kinematics::electromagnetic::kMinQ2Limit : controls::kMinQ2Limit; + + double md; + if(!interaction->ExclTag().KnownResonance()) md=1.23; + else { + Resonance_t res = interaction->ExclTag().Resonance(); + md=res::Mass(res); + } + + // ** 2-D Scan + const KPhaseSpace & kps = interaction->PhaseSpace(); + Range1D_t rW = kps.WLim(); + + int NW = 20; + double Wmin = rW.min + kASmallNum; + double Wmax = rW.max - kASmallNum; + + Wmax = TMath::Min(Wmax,fWcut); + + Wmin = TMath::Max(Wmin, md-.3); + Wmax = TMath::Min(Wmax, md+.3); + + if(Wmax-Wmin<0.05) { NW=1; Wmin=Wmax; } + + double dW = (NW>1) ? (Wmax-Wmin)/(NW-1) : 0.; + + for(int iw=0; iwKinePtr()->SetW(W); + + int NQ2 = 25; + int NQ2b = 4; + + Range1D_t rQ2 = kps.Q2Lim_W(); + if( rQ2.max < Q2Thres || rQ2.min <=0 ) continue; + if( rQ2.max-rQ2.min<0.02 ) {NQ2=5; NQ2b=3;} + + double logQ2min = TMath::Log(rQ2.min+kASmallNum); + double logQ2max = TMath::Log(rQ2.max-kASmallNum); + double dlogQ2 = (logQ2max - logQ2min) /(NQ2-1); + double xseclast = -1; + bool increasing = true; + + for(int iq2=0; iq2KinePtr()->SetQ2(Q2); + double xsec = fXSecModel->XSec(interaction, kPSWQD2fE); + LOG("DMRESKinematics", pDEBUG) + << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec; + max_xsec = TMath::Max(xsec, max_xsec); + increasing = xsec-xseclast>=0; + xseclast=xsec; + + // once the cross section stops increasing, I reduce the step size and + // step backwards a little bit to handle cases that the max cross section + // is grossly underestimated (very peaky distribution & large step) + if(!increasing) + { + dlogQ2/=NQ2b; + for(int iq2b=0; iq2bKinePtr()->SetQ2(Q2); + xsec = fXSecModel->XSec(interaction, kPSWQD2fE); + LOG("DMRESKinematics", pDEBUG) + << "xsec(W= " << W << ", Q2= " << Q2 << ") = " << xsec; + max_xsec = TMath::Max(xsec, max_xsec); + } + break; + } + } // Q2 + }//W + + // Apply safety factor, since value retrieved from the cache might + // correspond to a slightly different energy + // Apply larger safety factor for smaller energies. + max_xsec *= ( (El+Resonance) kinematics. + Is a concrete implementation of the EventRecordVisitorI interface. + Follows same procedure as RESKinematicsGenerator.h + +\author Zach Orr + Colorado State University + +\created November 18, 2004 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DMRES_KINEMATICS_GENERATOR_H_ +#define _DMRES_KINEMATICS_GENERATOR_H_ + +#include "Framework/Utils/Range1.h" +#include "Physics/Common/KineGeneratorWithCache.h" + +class TF2; + +namespace genie { + +class DMRESKinematicsGenerator : public KineGeneratorWithCache { + +public : + DMRESKinematicsGenerator(); + DMRESKinematicsGenerator(string config); + ~DMRESKinematicsGenerator(); + + // implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) 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 ComputeMaxXSec (const Interaction * interaction) const; + + mutable TF2 * fEnvelope; ///< 2-D envelope used for importance sampling + double fWcut; ///< Wcut parameter in DIS/RES join scheme +}; + +} // genie namespace +#endif // _DMRES_KINEMATICS_GENERATOR_H_ diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.cxx b/src/Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.cxx new file mode 100644 index 0000000000..4b0d2bbcb9 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.cxx @@ -0,0 +1,51 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + + Author: Zach Orr + Colorado State University + + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory +*/ +//____________________________________________________________________________ + +#include + +#include "Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.h" +#include "Framework/EventGen/EVGThreadException.h" +#include "Framework/GHEP/GHepRecord.h" +#include "Framework/GHEP/GHepParticle.h" +#include "Framework/GHEP/GHepFlags.h" +#include "Framework/Messenger/Messenger.h" + +using namespace genie; + +//___________________________________________________________________________ +DMRESOutgoingDarkGenerator::DMRESOutgoingDarkGenerator() : +OutgoingDarkGenerator("genie::DMRESOutgoingDarkGenerator") +{ + +} +//___________________________________________________________________________ +DMRESOutgoingDarkGenerator::DMRESOutgoingDarkGenerator(string config) : +OutgoingDarkGenerator("genie::DMRESOutgoingDarkGenerator", config) +{ + +} +//___________________________________________________________________________ +DMRESOutgoingDarkGenerator::~DMRESOutgoingDarkGenerator() +{ + +} +//___________________________________________________________________________ +void DMRESOutgoingDarkGenerator::ProcessEventRecord(GHepRecord * evrec) const +{ +// This method generates the final state primary lepton in RES events + + // no modification is required to the std implementation + OutgoingDarkGenerator::ProcessEventRecord(evrec); +} +//___________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.h b/src/Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.h new file mode 100644 index 0000000000..df482b736f --- /dev/null +++ b/src/Physics/BoostedDarkMatter/EventGen/DMRESOutgoingDarkGenerator.h @@ -0,0 +1,37 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESPrimaryLeptonGenerator + +\brief Generates the final state primary lepton in DM RES interactions. + Is a concrete implementation of the EventRecordVisitorI interface. + +\author Zach Orr Colorado State + +\created Jan 27 2024 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DMRES_OUTGOING_DARK_GENERATOR_H_ +#define _DMRES_OUTGOING_DARK_GENERATOR_H_ + +#include "Physics/Common/OutgoingDarkGenerator.h" + +namespace genie { + + class DMRESOutgoingDarkGenerator : public OutgoingDarkGenerator { + +public : + DMRESOutgoingDarkGenerator(); + DMRESOutgoingDarkGenerator(string config); + ~DMRESOutgoingDarkGenerator(); + + // implement the EventRecordVisitorI interface + void ProcessEventRecord(GHepRecord * event_rec) const; +}; + +} // genie namespace +#endif // _DMRES_OUTGOING_DARK_GENERATOR_H_ diff --git a/src/Physics/BoostedDarkMatter/EventGen/LinkDef.h b/src/Physics/BoostedDarkMatter/EventGen/LinkDef.h index 1a89a5ca5f..06dcfac82d 100644 --- a/src/Physics/BoostedDarkMatter/EventGen/LinkDef.h +++ b/src/Physics/BoostedDarkMatter/EventGen/LinkDef.h @@ -10,9 +10,16 @@ #pragma link C++ class genie::DMELInteractionListGenerator; #pragma link C++ class genie::DMELKinematicsGenerator; #pragma link C++ class genie::DMELOutgoingDarkGenerator; + #pragma link C++ class genie::DMDISOutgoingDarkGenerator; #pragma link C++ class genie::DMDISInteractionListGenerator; #pragma link C++ class genie::DMDISKinematicsGenerator; + +#pragma link C++ class genie::DMRESInteractionListGenerator; +#pragma link C++ class genie::DMRESKinematicsGenerator; +#pragma link C++ class genie::DMRESOutgoingDarkGenerator; +#pragma link C++ class genie::DMRESHadronicSystemGenerator; + #pragma link C++ class genie::DMEOutgoingDarkGenerator; #pragma link C++ class genie::DMEInteractionListGenerator; #pragma link C++ class genie::DMEKinematicsGenerator; diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESPXSec.cxx b/src/Physics/BoostedDarkMatter/XSection/DMRESPXSec.cxx new file mode 100644 index 0000000000..224c401e94 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESPXSec.cxx @@ -0,0 +1,514 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Zachary W. Orr + Colorado State University + +based on code by + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory +*/ +//____________________________________________________________________________ + +#include +#include + +#include "Framework/Algorithm/AlgFactory.h" +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/RefFrame.h" +#include "Framework/Conventions/KineVar.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/Numerical/Spline.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGLibrary.h" +#include "Framework/Utils/KineUtils.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/Utils/Range1.h" +#include "Framework/Utils/BWFunc.h" +#include "Physics/BoostedDarkMatter/XSection/DMRESPXSec.h" +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h" +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h" +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Physics/NuclearState/FermiMomentumTablePool.h" +#include "Physics/NuclearState/FermiMomentumTable.h" +#include "Physics/NuclearState/NuclearUtils.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +DMRESPXSec::DMRESPXSec() : +XSecAlgorithmI("genie::DMRESPXSec") +{ + +} +//____________________________________________________________________________ +DMRESPXSec::DMRESPXSec(string config) : +XSecAlgorithmI("genie::DMRESPXSec", config) +{ + +} +//____________________________________________________________________________ +DMRESPXSec::~DMRESPXSec() +{ + +} +//____________________________________________________________________________ +double DMRESPXSec::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(); + const ProcessInfo & proc_info = interaction -> ProcInfo(); + const Target & target = init_state.Tgt(); + + const Kinematics & kinematics = interaction -> Kine(); + LOG("DMRESPXSec", pDEBUG) << "Using v^" << fVelMode << " dependence"; + +//Get Kinematic Parameters: + double W = kinematics.W(); + double W2 = TMath::Power(W, 2); //W^2 + double q2 = kinematics.q2(); + +//Rejection Sampling from Kinematics: +// Under the DIS/RES joining scheme, xsec(RES)=0 for W>=Wcut + if(fUsingDisResJoin) { + if(W>=fWcut) { +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DMRES", pDEBUG) + << "RES/DIS Join Scheme: XSec[RES, W=" << W + << " >= Wcut=" << fWcut << "] = 0"; +#endif + return 0; + } + } +//Pass: Resonance XSec + +//Probe Res Interaction: +// Get the input baryon resonance: + Resonance_t resonance = interaction->ExclTag().Resonance(); + string resname = utils::res::AsString(resonance); + bool is_delta = utils::res::IsDelta (resonance); + +// Get the dark matter, target nucleon & DM current: + int nucpdgc = target.HitNucPdg(); + int probepdgc = init_state.ProbePdg(); + bool is_dm = pdg::IsDarkMatter (probepdgc); + bool is_dmbar = pdg::IsAntiDarkMatter (probepdgc); + + + bool is_p = pdg::IsProton (nucpdgc); + bool is_n = pdg::IsNeutron (nucpdgc); + bool is_dmres = proc_info.IsDarkMatterResonant(); + + //Need break for is not DM + if(!is_dmres) { + return 0; + } + + // Get baryon resonance parameters + int IR = utils::res::ResonanceIndex (resonance); + int LR = utils::res::OrbitalAngularMom (resonance); + double MR = utils::res::Mass (resonance); + double WR = utils::res::Width (resonance); + double NR = fNormBW?utils::res::BWNorm (resonance,fN0ResMaxNWidths,fN2ResMaxNWidths,fGnResMaxNWidths):1; + + // Following NeuGEN, avoid problems with underlying unphysical + // model assumptions by restricting the allowed W phase space + // around the resonance peak + if (fNormBW) { + if (W > MR + fN0ResMaxNWidths * WR && IR==0) return 0.; + else if (W > MR + fN2ResMaxNWidths * WR && IR==2) return 0.; + else if (W > MR + fGnResMaxNWidths * WR) return 0.; + } + + //Auxillary Kinematic factors + double E = init_state.ProbeE(kRfHitNucRest); //E1 + double E2 = TMath::Power(E,2); //E1^2 + double Mnuc = target.HitNucMass(); //Nucleon mass + double Mnuc2 = TMath::Power(Mnuc, 2); //m^2 + + //Scalar Charge + double fQchiS2 = TMath::Power(fQchiS,2); //QS^2 + // L and R charges + double fQchiL = fQchiV-fQchiA; + double fQchiR = fQchiV+fQchiA; + + // Compute kinematic factors + double k = 0.5 * (W2 - Mnuc2)/Mnuc; + double v = k - 0.5 * q2/Mnuc; + double v2 = TMath::Power(v, 2); + double Q2 = v2 - q2; + double Q = TMath::Sqrt(Q2); + double Eprime = E - v; + + double U = 0.5 * (E + Eprime + Q) / E; + double V = 0.5 * (E + Eprime - Q) / E; + double U2 = TMath::Power(U, 2); + double V2 = TMath::Power(V, 2); + double UV = U*V; + + +// Masses + double mchi = init_state.GetProbeP4(kRfHitNucRest)->M(); //DM mass: mx + double mchi2 = TMath::Power(mchi, 2); // mx^2 + double mZprime2 = TMath::Power(fMedMass, 2); // mZ'^2 + double mZprime4 = TMath::Power(mZprime2, 2); //mZ'^4 +//DM charge auxillary factors + double fQchiLmR = fQchiL-fQchiR; // QL - QR + double fQchiLR = fQchiL * fQchiR; // QL*QR + double fQchiLmR2 = TMath::Power(fQchiLmR, 2); //(QL - QR)^2 + double fQchiL2 = TMath::Power(fQchiL, 2); //QL^2 + double fQchiR2 = TMath::Power(fQchiR, 2); //QR^2 +//DM mass auxillary factors + double mZ_q2 = TMath::Power(mZprime2 - q2, 2)/mZprime4; //(mZ'^2 - q^2)^2 / mZ'^4 + double mchiTerm = (mchi2 * Q2)/(E2 * q2); //mchi^2 * Q^2 / E^2 * q^2 + + //Spin average over incident DM = 0.5 + double AL = U2*fQchiL2 + V2*fQchiR2 + 2*mchiTerm*fQchiLR; + double AR = V2*fQchiL2 + U2*fQchiR2 + 2*mchiTerm*fQchiLR; + double AS = 2*UV*(fQchiL2 + fQchiR2) + mchiTerm*fQchiLmR2; + //DM term + double AZ = 0.0; + if (mZprime2 != 0.0){ + AZ = -mchiTerm*mZ_q2*fQchiLmR2; + } + + //For Scalar DM: spin average = 1 + if(fVelMode == 2){ + AL = fQchiS2 * (2*UV + 2*mchiTerm); + AR = fQchiS2 * (2*UV + 2*mchiTerm); + AS = fQchiS2 * (TMath::Power((U+V), 2)); + AZ = 0.0; +} + // Calculate the Feynman-Kislinger-Ravndall parameters + double Go = TMath::Power(1 - 0.25 * q2/Mnuc2, 0.5-IR); + double GV = Go * TMath::Power( 1./(1-q2/fMv2), 2); + double GA = Go * TMath::Power( 1./(1-q2/fMa2), 2); + if(mZprime2 == 0.){ + GA = 0; + } + + double d = TMath::Power(W+Mnuc,2.) - q2; + double d2 = W2 - Mnuc2 + q2; + double sq2omg = TMath::Sqrt(2./fOmega); + double nomg = IR * fOmega; + double mq_w = Mnuc*Q/W; + double BC = (2*W*mq_w) / (W2 - Mnuc2 + q2); + fFKRDM.Lamda = sq2omg * mq_w; + fFKRDM.Tv = GV / (3.*W*sq2omg); + fFKRDM.Rv = kSqrt2 * mq_w*(W+Mnuc)*GV / d; + fFKRDM.S = (-q2/Q2) * (3*W*Mnuc + q2 - Mnuc2) * GV / (6*Mnuc2); + fFKRDM.Ta = (2./3.) * (fZeta/sq2omg) * mq_w * GA / d; + fFKRDM.Ra = (kSqrt2/6.) * fZeta * (GA/W) * (W+Mnuc + 2*nomg*W/d ); + fFKRDM.Bs = fZeta/(3.*W*sq2omg) * (1 + (W2-Mnuc2+q2)/ d) * GA; + fFKRDM.Cs = fZeta/(6.*Q) * (W2 - Mnuc2 + nomg*(W2-Mnuc2+q2)/d) * (GA/Mnuc); + fFKRDM.Bz = ((fZeta * GA) * (W - Mnuc)) / (2 * W * mq_w * sq2omg); + fFKRDM.Cz = ((fZeta * GA)/3) * (1 + ((3*q2 + nomg) / d)); + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("FKR", pDEBUG) + << "FKR params for DMRES = " << resname << " : " << fFKRDM; +#endif + + // Calculate the Rein-Sehgal Helicity Amplitudes + + const RSHelicityAmplModelDMI * hamplmod = 0; + if(is_dmres) { + if (is_p) { hamplmod = fHAmplModelDMp;} + else { hamplmod = fHAmplModelDMn;} + } + assert(hamplmod); + + const RSHelicityAmplDM & hampl = hamplmod->Compute(resonance, fFKRDM); + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DMHAmpl", pDEBUG) + << "Helicity Amplitudes for DMRES = " << resname << " : " << hampl; +#endif + + // Compute the cross section structure factors + double gZp4 = TMath::Power(fgZp, 4.0); + double XoEn = 1.0-(mchi2/E2); //1 - (mchi/E)^2 + double XoPropMass = TMath::Power(q2 - mZprime2, 2); //(q^2 - mZ'^2)^2 + double Xo = gZp4/(XoEn*XoPropMass); + double sig0 = 0.0625*(Xo/kPi)*(-q2/Q2)*(W/Mnuc); + double scLR = 0.5*W/Mnuc; + double scS = 0.5*(Mnuc/W)*(Q2/(-q2)); + double sigL = 0.0; + double sigR = 0.0; + double sigS = 0.0; + double sigZ = 0.0; + +//Including Hadron spin avg and res-rest frame transformations + sigL = scLR* (hampl.Amp2Plus3() + hampl.Amp2Plus1()); + sigR = scLR* (hampl.Amp2Minus3 () + hampl.Amp2Minus1 ()); + sigS = scS * (hampl.Amp20Plus () + hampl.Amp20Minus()); + sigZ = scS *(hampl.Ampz20Plus () + hampl.Ampz20Minus()); + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DMRES", pDEBUG) << "sig_{0} = " << sig0; + LOG("DMRES", pDEBUG) << "sig_{L} = " << sigL; + LOG("DMRES", pDEBUG) << "sig_{R} = " << sigR; + LOG("DMRES", pDEBUG) << "sig_{S} = " << sigS; + LOG("DMRES", pDEBUG) << "sig_{Z} = " << sigZ; +#endif + +//DM Cross Section Calculation: + +//XSec calc: + double xsec = 0.0; + if (is_dm) { + xsec = sig0*(AL*sigL + AR*sigR + AS*sigS + AZ*sigZ); + } + else + if (is_dmbar) { + xsec = sig0*(AL*sigR + AR*sigL + AS*sigS + AZ*sigZ); + } + xsec = TMath::Max(0.,xsec); + + +//________________________________________________________________________________________ +//X-Sec scaling by BreitWigner,Jacobian,free-nucleon,Pauli-Blocking +//________________________________________________________________________________________ + + + // Check whether the cross section is to be weighted with a + // Breit-Wigner distribution (default: true) + double bw = 1.0; + if(fWghtBW) { + bw = utils::bwfunc::BreitWignerL(W,LR,MR,WR,NR); + + } +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DMRES", pDEBUG) + << "BreitWigner(RES=" << resname << ", W=" << W << ") = " << bw; +#endif + xsec *= bw; + + + //Apply given scaling factor + double xsec_scale = 1.; + if (is_dmres) { xsec_scale = fXSecScaleDM; } + xsec *= xsec_scale; + +#ifdef __GENIE_LOW_LEVEL_MESG_ENABLED__ + LOG("DMRES", pINFO) + << "\n d2xsec/dQ2dW" << "[" << interaction->AsString() + << "](W=" << W << ", q2=" << q2 << ", E=" << E << ") = " << xsec; +#endif + + // The algorithm computes d^2xsec/dWdQ2 + // Check whether variable tranformation is needed + if(kps!=kPSWQ2fE) { + double J = utils::kinematics::Jacobian(interaction,kPSWQ2fE,kps); + xsec *= J; + } + + // If requested return the free nucleon xsec even for input nuclear tgt + if( interaction->TestBit(kIAssumeFreeNucleon) ) return xsec; + + int Z = target.Z(); + int A = target.A(); + int N = A-Z; + + // Take into account the number of scattering centers in the target + int NNucl = (is_p) ? Z : N; + + xsec*=NNucl; // nuclear xsec (no nuclear suppression factor) + + if (fUsePauliBlocking && A!=1) + { + // Calculation of Pauli blocking according references: + // + // [1] S.L. Adler, S. Nussinov, and E.A. Paschos, "Nuclear + // charge exchange corrections to leptonic pion production + // in the (3,3) resonance region," Phys. Rev. D 9 (1974) + // 2125-2143 [Erratum Phys. Rev. D 10 (1974) 1669]. + // [2] J.Y. Yu, "Neutrino interactions and nuclear effects in + // oscillation experiments and the nonperturbative disper- + // sive sector in strong (quasi-)abelian fields," Ph. D. + // Thesis, Dortmund U., Dortmund, 2002 (unpublished). + // [3] E.A. Paschos, J.Y. Yu, and M. Sakuda, "Neutrino pro- + // duction of resonances," Phys. Rev. D 69 (2004) 014013 + // [arXiv: hep-ph/0308130]. + + double P_Fermi = 0.0; + + // Maximum value of Fermi momentum of target nucleon (GeV) + if (A<6 || !fUseRFGParametrization) + { + // Look up the Fermi momentum for this target + FermiMomentumTablePool * kftp = FermiMomentumTablePool::Instance(); + const FermiMomentumTable * kft = kftp->GetTable(fKFTable); + P_Fermi = kft->FindClosestKF(pdg::IonPdgCode(A, Z), nucpdgc); + } + else { + // Define the Fermi momentum for this target + P_Fermi = utils::nuclear::FermiMomentumForIsoscalarNucleonParametrization(target); + // Correct the Fermi momentum for the struck nucleon + if(is_p) { P_Fermi *= TMath::Power( 2.*Z/A, 1./3); } + else { P_Fermi *= TMath::Power( 2.*N/A, 1./3); } + } + + double FactorPauli_RES = 1.0; + + double k0 = 0., q = 0., q0 = 0.; + + if (P_Fermi > 0.) + { + k0 = (W2-Mnuc2-Q2)/(2*W); + k = TMath::Sqrt(k0*k0+Q2); // previous value of k is overridden + q0 = (W2-Mnuc2+kPionMass2)/(2*W); + q = TMath::Sqrt(q0*q0-kPionMass2); + } + + if (2*P_Fermi < k-q) + FactorPauli_RES = 1.0; + if (2*P_Fermi >= k+q) + FactorPauli_RES = ((3*k*k+q*q)/(2*P_Fermi)-(5*TMath::Power(k,4)+TMath::Power(q,4)+10*k*k*q*q)/(40*TMath::Power(P_Fermi,3)))/(2*k); + if (2*P_Fermi >= k-q && 2*P_Fermi <= k+q) + FactorPauli_RES = ((q+k)*(q+k)-4*P_Fermi*P_Fermi/5-TMath::Power(k-q, 3)/(2*P_Fermi)+TMath::Power(k-q, 5)/(40*TMath::Power(P_Fermi, 3)))/(4*q*k); + + xsec *= FactorPauli_RES; + } + + return xsec; +} +//________________________________________________________________________________________ +//________________________________________________________________________________________ + + +//____________________________________________________________________________ +double DMRESPXSec::Integral(const Interaction * interaction) const +{ + double xsec = fXSecIntegrator->Integrate(this,interaction); + return xsec; +} +//____________________________________________________________________________ +bool DMRESPXSec::ValidProcess(const Interaction * interaction) const +{ + if(interaction->TestBit(kISkipProcessChk)) return true; + + const InitialState & init_state = interaction->InitState(); + const ProcessInfo & proc_info = interaction->ProcInfo(); + const XclsTag & xcls = interaction->ExclTag(); + + if(!proc_info.IsDarkMatterResonant()) return false; + if(!xcls.KnownResonance()) return false; + + int hitnuc = init_state.Tgt().HitNucPdg(); + bool is_pn = (pdg::IsProton(hitnuc) || pdg::IsNeutron(hitnuc)); + + if (!is_pn) return false; + + int probe = init_state.ProbePdg(); + bool is_dmres = proc_info.IsDarkMatterResonant(); + bool is_dm = proc_info.IsDarkMatter(); + + if (!is_dm && !is_dmres) return false; + + return true; +} +//____________________________________________________________________________ +void DMRESPXSec::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void DMRESPXSec::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void DMRESPXSec::LoadConfig(void) +{ + // Cross section scaling factor + this->GetParam( "RES-DM-XSecScale", fXSecScaleDM ) ; + + //FKR Params + this->GetParam( "RES-Zeta", fZeta ) ; + this->GetParam( "RES-Omega", fOmega ) ; + + // velocity dependence of interaction: fermion/scalar DM + this->GetParamDef("velocity-mode", fVelMode, 0 ); + //DM mass + double ma, mv ; + this->GetParam( "RES-Ma", ma ) ; + this->GetParam( "RES-Mv", mv ) ; + fMa2 = TMath::Power(ma,2); + fMv2 = TMath::Power(mv,2); + //DM Charges + double QchiL, QchiR; + this->GetParam( "DarkLeftCharge", QchiL ) ; + this->GetParam( "DarkRightCharge", QchiR ) ; + this->GetParam( "DarkScalarCharge", fQchiS ) ; + fQchiV = 0.5*(QchiL + QchiR); + fQchiA = 0.5*(- QchiL + QchiR); + // DM Mediator + this->GetParam("ZpCoupling", fgZp ) ; + fMedMass = PDGLibrary::Instance()->Find(kPdgMediator)->Mass(); + + this->GetParamDef( "BreitWignerWeight", fWghtBW, true ) ; + this->GetParamDef( "BreitWignerNorm", fNormBW, true); + + + this->GetParam("FermiMomentumTable", fKFTable); + this->GetParam("RFG-UseParametrization", fUseRFGParametrization); + this->GetParam("UsePauliBlockingForRES", fUsePauliBlocking); + + // Load all the sub-algorithms needed + + fHAmplModelDMp = 0; + fHAmplModelDMn = 0; + + AlgFactory * algf = AlgFactory::Instance(); + + fHAmplModelDMp = dynamic_cast ( + algf->GetAlgorithm("genie::RSHelicityAmplModelDMp","Default")); //DM + p + fHAmplModelDMn = dynamic_cast ( + algf->GetAlgorithm("genie::RSHelicityAmplModelDMn","Default")); //DM + n + + assert( fHAmplModelDMp ); + assert( fHAmplModelDMn ); + + // Use algorithm within a DIS/RES join scheme. If yes get Wcut + this->GetParam( "UseDRJoinScheme", fUsingDisResJoin ) ; + fWcut = 999999; + if(fUsingDisResJoin) { + this->GetParam( "Wcut", fWcut ) ; + } + + double thw ; + this->GetParam( "WeinbergAngle", thw ) ; + fSin48w = TMath::Power( TMath::Sin(thw), 4 ); + // NeuGEN limits in the allowed resonance phase space: + // W < min{ Wmin(physical), (res mass) + x * (res width) } + // It limits the integration area around the peak and avoids the + // problem with huge xsec increase at low Q2 and high W. + // In correspondence with Hugh, Rein said that the underlying problem + // are unphysical assumptions in the model. + this->GetParamDef( "MaxNWidthForN2Res", fN2ResMaxNWidths, 2.0 ) ; + this->GetParamDef( "MaxNWidthForN0Res", fN0ResMaxNWidths, 6.0 ) ; + this->GetParamDef( "MaxNWidthForGNRes", fGnResMaxNWidths, 4.0 ) ; + +//______________________________________________________________________________________ + + + // load the differential cross section integrator + fXSecIntegrator = + dynamic_cast (this->SubAlg("XSec-Integrator")); + assert(fXSecIntegrator); +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESPXSec.h b/src/Physics/BoostedDarkMatter/XSection/DMRESPXSec.h new file mode 100644 index 0000000000..3ee8456c21 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESPXSec.h @@ -0,0 +1,107 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESPXSec + +\brief Computes the double differential cross section for resonance + DM-production according to the Rein-Sehgal model. + + The computed cross section is the d^2 xsec/ dQ^2 dW \n + + where \n + \li \c Q^2 : momentum transfer ^ 2 + \li \c W : invariant mass of the final state hadronic system + + Is a concrete implementation of the XSecAlgorithmI interface. + +\ref D.Rein and L.M.Sehgal, Neutrino Excitation of Baryon Resonances + and Single Pion Production, Ann.Phys.133, 79 (1981) + +\author Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + + Changes made for Dark Matter - Zach Orr, Colorado State University + +\created May 05, 2004 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DM_RES_PXSEC_H_ +#define _DM_RES_PXSEC_H_ + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Framework/ParticleData/BaryonResonance.h" +#include "Physics/BoostedDarkMatter/XSection/FKRDM.h" + +namespace genie { + +class RSHelicityAmplModelDMI; +class Spline; +class XSecIntegratorI; + +class DMRESPXSec : public XSecAlgorithmI { + +public: + DMRESPXSec(); + DMRESPXSec(string config); + virtual ~DMRESPXSec(); + + // implement the XSecAlgorithmI interface + 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); + + mutable FKRDM fFKRDM; + + const RSHelicityAmplModelDMI * fHAmplModelDMp; + const RSHelicityAmplModelDMI * fHAmplModelDMn; + + // configuration data + double fQchiV; ///< DM vector charge + double fQchiA; ///< DM axial charge + double fQchiS; ///< DM scalar charge + int fVelMode; ///< DM (scalar/fermion) + double fMedMass; ///< DM mediator mass + double fgZp; ///< DM mediator coupling + bool fWghtBW; ///< weight with resonance breit-wigner? + bool fNormBW; ///< normalize resonance breit-wigner to 1? + double fZeta; ///< FKR parameter Zeta + double fOmega; ///< FKR parameter Omega + double fMa2; ///< (axial mass)^2 + double fMv2; ///< (vector mass)^2 + bool fUsingDisResJoin; ///< use a DIS/RES joining scheme? + double fWcut; ///< apply DIS/RES joining scheme < Wcut + double fN2ResMaxNWidths; ///< limits allowed phase space for n=2 res + double fN0ResMaxNWidths; ///< limits allowed phase space for n=0 res + double fGnResMaxNWidths; ///< limits allowed phase space for other res + string fKFTable; ///< table of Fermi momentum (kF) constants for various nuclei + bool fUseRFGParametrization; ///< use parametrization for fermi momentum insted of table? + bool fUsePauliBlocking; ///< account for Pauli blocking? + + double fSin48w; ///< sin^4(Weingberg angle) + +// bool fUsingNuTauScaling; ///< use NeuGEN nutau xsec reduction factors? +// Spline * fNuTauRdSpl; ///< xsec reduction spline for nu_tau +// Spline * fNuTauBarRdSpl; ///< xsec reduction spline for nu_tau_bar + + double fXSecScaleDM; ///< external DM xsec scaling factor + + + const XSecIntegratorI * fXSecIntegrator; +}; + +} // genie namespace + +#endif // _DM_RES_PXSEC_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESXSec.cxx b/src/Physics/BoostedDarkMatter/XSection/DMRESXSec.cxx new file mode 100644 index 0000000000..ec6712987d --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESXSec.cxx @@ -0,0 +1,121 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + based on code of + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + + + Changes made for RES DM by Zach Orr, Colorado State University + */ + //____________________________________________________________________________ + + #include + #include + #include + + #include "Framework/Conventions/GBuild.h" + #include "Framework/Conventions/Constants.h" + #include "Framework/Conventions/Units.h" + #include "Framework/Conventions/KineVar.h" + #include "Framework/Messenger/Messenger.h" + #include "Framework/ParticleData/PDGUtils.h" + #include "Framework/Numerical/MathUtils.h" + #include "Framework/Utils/KineUtils.h" + #include "Framework/Numerical/GSLUtils.h" + #include "Physics/BoostedDarkMatter/XSection/DMRESXSec.h" + #include "Physics/XSectionIntegration/GSLXSecFunc.h" + + using namespace genie; + using namespace genie::constants; + + //____________________________________________________________________________ + DMRESXSec::DMRESXSec() : + XSecIntegratorI("genie::DMRESXSec") + { + + } + //____________________________________________________________________________ + DMRESXSec::DMRESXSec(string config) : + XSecIntegratorI("genie::DMRESXSec", config) + { + + } + //____________________________________________________________________________ + DMRESXSec::~DMRESXSec() + { + + } + //____________________________________________________________________________ + double DMRESXSec::Integrate( + const XSecAlgorithmI * model, const Interaction * in) const + { + if(! model->ValidProcess(in) ) return 0.; + + const KPhaseSpace & kps = in->PhaseSpace(); + if(!kps.IsAboveThreshold()) { + LOG("COHXSec", pDEBUG) << "*** Below energy threshold"; + return 0; + } + Range1D_t Q2l = kps.Limits(kKVQ2); + Range1D_t Wl = kps.Limits(kKVW); + + Interaction * interaction = new Interaction(*in); + interaction->SetBit(kISkipProcessChk); + //interaction->SetBit(kISkipKinematicChk); + + ROOT::Math::IBaseFunctionMultiDim * func = + new utils::gsl::d2XSec_dWdQ2_E(model, interaction); + + ROOT::Math::IntegrationMultiDim::Type ig_type = + utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); + double abstol = 1E-16; //We mostly care about relative tolerance. + + ROOT::Math::IntegratorMultiDim ig(*func, ig_type, abstol, fGSLRelTol, fGSLMaxEval); + + double kine_min[2] = { Wl.min, Q2l.min }; + double kine_max[2] = { Wl.max, Q2l.max }; + double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); + + LOG("DMRESXSec", pERROR) << "Integrator opt / Integrator = " << ig.Options().Integrator(); + + if(xsec < 0) { + LOG("DMRESXSec", pERROR) << "Algorithm " << *model << " returns a negative cross-section (xsec = " << xsec << " 1E-38 * cm2)"; + LOG("DMRESXSec", pERROR) << "for process" << *interaction; + LOG("DMRESXSec", pERROR) << "Integrator status code = " << ig.Status(); + LOG("DMRESXSec", pERROR) << "Integrator error code = " << ig.Error(); + } + + //LOG("RESXSec", pINFO) << "XSec[RES] (Ev = " << Ev << " GeV) = " << xsec; + + delete interaction; + delete func; + return xsec; + } + //____________________________________________________________________________ + void DMRESXSec::Configure(const Registry & config) + { + Algorithm::Configure(config); + this->LoadConfig(); + } + //____________________________________________________________________________ + void DMRESXSec::Configure(string config) + { + Algorithm::Configure(config); + this->LoadConfig(); + } + //____________________________________________________________________________ + void DMRESXSec::LoadConfig(void) + { + // Get GSL integration type & relative tolerance + GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ; + GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 1E-2 ) ; + int max, min ; + GetParamDef( "gsl-max-eval", max, 500000 ) ; + GetParamDef( "gsl-min-eval", min, 5000 ) ; + fGSLMaxEval = (unsigned int) max ; + fGSLMinEval = (unsigned int) min ; + } + //____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESXSec.h b/src/Physics/BoostedDarkMatter/XSection/DMRESXSec.h new file mode 100644 index 0000000000..37360aeb33 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESXSec.h @@ -0,0 +1,48 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESXSec + +\brief Computes the DM RES Cross Section.\n + Is a concrete implementation of the XSecIntegratorI interface.\n + +\author Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + + updated for DM RES by Zach Orr, Colorado State University + +\created May 04, 2004 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DMRES_XSEC_H_ +#define _DMRES_XSEC_H_ + +#include "Physics/XSectionIntegration/XSecIntegratorI.h" + +namespace genie { + +class DMRESXSec : public XSecIntegratorI { + +public: + DMRESXSec(); + DMRESXSec(string param_set); + virtual ~DMRESXSec(); + + //! 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 // _DMRES_XSEC_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESXSecFast.cxx b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecFast.cxx new file mode 100644 index 0000000000..4a4e0f2061 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecFast.cxx @@ -0,0 +1,227 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2026, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Igor Kakorin + Joint Institute for Nuclear Research + + based on code of + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + + adjusted for Dark Matter by + Zachary W. Orr + Colorado State University + */ +//____________________________________________________________________________ + +#include +#include +#include + +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Conventions/KineVar.h" +#include "Physics/XSectionIntegration/GSLXSecFunc.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Physics/BoostedDarkMatter/XSection/DMRESXSecFast.h" +#include "Framework/Utils/RunOpt.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/Utils/KineUtils.h" +#include "Framework/Utils/Cache.h" +#include "Framework/Utils/CacheBranchFx.h" +#include "Framework/Utils/XSecSplineList.h" +#include "Framework/Numerical/GSLUtils.h" + +using namespace genie; +using namespace genie::constants; +using namespace genie::units; + +//____________________________________________________________________________ +DMRESXSecFast::DMRESXSecFast() : +DMRESXSecWithCacheFast("genie::DMRESXSecFast") +{ + +} +//____________________________________________________________________________ +DMRESXSecFast::DMRESXSecFast(string config) : +DMRESXSecWithCacheFast("genie::DMRESXSecFast", config) +{ + +} +//____________________________________________________________________________ +DMRESXSecFast::~DMRESXSecFast() +{ + +} +//____________________________________________________________________________ +double DMRESXSecFast::Integrate( + const XSecAlgorithmI * model, const Interaction * interaction) const +{ + if(! model->ValidProcess(interaction) ) return 0.; + fSingleResXSecModel = model; + + const KPhaseSpace & kps = interaction->PhaseSpace(); + if(!kps.IsAboveThreshold()) { + LOG("DMRESXSecFast", pDEBUG) << "*** Below energy threshold"; + return 0; + } + + //-- Get init state and process information + const InitialState & init_state = interaction->InitState(); + const ProcessInfo & proc_info = interaction->ProcInfo(); + const Target & target = init_state.Tgt(); + + InteractionType_t it = proc_info.InteractionTypeId(); + int nucleon_pdgc = target.HitNucPdg(); + int nu_pdgc = init_state.ProbePdg(); + + //-- Get neutrino energy in the struck nucleon rest frame + double Ev = init_state.ProbeE(kRfHitNucRest); + + //-- Get the requested resonance + Resonance_t res = interaction->ExclTag().Resonance(); + + // If the input interaction is off a nuclear target, then chek whether + // the corresponding free nucleon cross section already exists at the + // cross section spline list. + // If yes, calculate the nuclear cross section based on that value. + // + XSecSplineList * xsl = XSecSplineList::Instance(); + if(init_state.Tgt().IsNucleus() && !xsl->IsEmpty() ) { + Interaction * in = new Interaction(*interaction); + if(pdg::IsProton(nucleon_pdgc)) { + in->InitStatePtr()->TgtPtr()->SetId(kPdgTgtFreeP); + } else { + in->InitStatePtr()->TgtPtr()->SetId(kPdgTgtFreeN); + } + if(xsl->SplineExists(model,in)) { + const Spline * spl = xsl->GetSpline(model, in); + double xsec = spl->Evaluate(Ev); + SLOG("DMResTF", pNOTICE) + << "XSec[RES/" << utils::res::AsString(res)<< "/free] (Ev = " + << Ev << " GeV) = " << xsec/(1E-38 *cm2)<< " x 1E-38 cm^2"; + if(! interaction->TestBit(kIAssumeFreeNucleon) ) { + int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N(); + xsec *= NNucl; + } + delete in; + return xsec; + } + delete in; + } + + // There was no corresponding free nucleon spline saved in XSecSplineList that + // could be used to speed up this calculation. + // Check whether local caching of free nucleon cross sections is allowed. + // If yes, store free nucleon cross sections at a cache branch and use those + // at any subsequent call. + // + bool bare_xsec_pre_calc = RunOpt::Instance()->BareXSecPreCalc(); + if(bare_xsec_pre_calc && !fUsePauliBlocking) { + Cache * cache = Cache::Instance(); + string key = this->CacheBranchName(res, it, nu_pdgc, nucleon_pdgc); + LOG("DMResTF", pINFO) + << "Finding cache branch with key: " << key; + CacheBranchFx * cache_branch = + dynamic_cast (cache->FindCacheBranch(key)); + if(!cache_branch) { + LOG("DMResTF", pWARN) + << "No cached RES v-production data for input neutrino" + << " (pdgc: " << nu_pdgc << ")"; + LOG("DMResTF", pWARN) + << "Wait while computing/caching RES production xsec first..."; + + this->CacheResExcitationXSec(interaction); + + LOG("DMResTF", pINFO) << "Done caching resonance xsec data"; + LOG("DMResTF", pINFO) + << "Finding newly created cache branch with key: " << key; + cache_branch = + dynamic_cast (cache->FindCacheBranch(key)); + assert(cache_branch); + } + const CacheBranchFx & cbranch = (*cache_branch); + + //-- Get cached resonance neutrinoproduction xsec + // (If E>Emax, assume xsec = xsec(Emax) - but do not evaluate the + // cross section spline at the end of its energy range-) + double rxsec = (EvTestBit(kIAssumeFreeNucleon) ) return rxsec; + + int NNucl = (pdg::IsProton(nucleon_pdgc)) ? target.Z() : target.N(); + rxsec*=NNucl; // nuclear xsec + return rxsec; + } // disable local caching + + // Just go ahead and integrate the input differential cross section for the + // specified interaction. + else { + + Range1D_t rW = Range1D_t(0.0,1.0); + Range1D_t rQ2 = Range1D_t(0.0,1.0); + + LOG("DMResTF", pINFO) + << "*** Integrating d^2 XSec/dWdQ^2 for R: " + << utils::res::AsString(res) << " at Ev = " << Ev; + LOG("DMResTF", pINFO) + << "{W} = " << rW.min << ", " << rW.max; + LOG("DMResTF", pINFO) + << "{Q^2} = " << rQ2.min << ", " << rQ2.max; + + ROOT::Math::IBaseFunctionMultiDim * func= new utils::gsl::d2XSecRESFast_dWQ2_E(model, interaction); + ROOT::Math::IntegrationMultiDim::Type ig_type = utils::gsl::IntegrationNDimTypeFromString(fGSLIntgType); + ROOT::Math::IntegratorMultiDim ig(ig_type,0,fGSLRelTol,fGSLMaxEval); + ig.SetFunction(*func); + double kine_min[2] = { rW.min, rQ2.min }; + double kine_max[2] = { rW.max, rQ2.max }; + double xsec = ig.Integral(kine_min, kine_max) * (1E-38 * units::cm2); + + delete func; + return xsec; + } + return 0; +} +//____________________________________________________________________________ +void DMRESXSecFast::Configure(const Registry & config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void DMRESXSecFast::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void DMRESXSecFast::LoadConfig(void) +{ + + // Get GSL integration type & relative tolerance + GetParamDef( "gsl-integration-type", fGSLIntgType, string("adaptive") ) ; + GetParamDef( "gsl-relative-tolerance", fGSLRelTol, 0.01 ) ; + GetParamDef( "gsl-max-eval", fGSLMaxEval, 100000 ) ; + GetParam("UsePauliBlockingForRES", fUsePauliBlocking); + // Get upper E limit on res xsec spline (=f(E)) before assuming xsec=const + GetParamDef( "ESplineMax", fEMax, 100. ) ; + fEMax = TMath::Max(fEMax, 20.); // don't accept user Emax if less than 20 GeV + + // Create the baryon resonance list specified in the config. + fResList.Clear(); + string resonances ; + GetParam( "ResonanceNameList", resonances ) ; + fResList.DecodeFromNameList(resonances); +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESXSecFast.h b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecFast.h new file mode 100644 index 0000000000..3581c03f33 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecFast.h @@ -0,0 +1,69 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESXSecFast + +\brief Computes the cross section for an exclusive 1pi reaction through + resonance neutrinoproduction according to the Rein-Sehgal model. + + This algorithm produces in principle what you could also get from + the genie::RESXSec algorithm (RES cross section integrator) by + specifying the genie::DMRESPXSec as the differential + cross section model. However, DMRESXSecFast offers a faster + alternative. Before computing any RES cross section this algorithm + computes and caches splines for resonance neutrino-production cross + sections. This improves the speed of the GENIE spline construction + phase if splines for multiple nuclear targets are to be computed. + Also this class integrates cross sections faster, than + DMRESXSec because of integration area transformation. + + Is a concrete implementation of the XSecAlgorithmI interface.\n + +\ref D.Rein and L.M.Sehgal, Neutrino Excitation of Baryon Resonances + and Single Pion Production, Ann.Phys.133, 79 (1981) + +\author Igor Kakorin + Joint Institute for Nuclear Research + + based on code of + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + +\created March 01, 2017 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + +*/ +//____________________________________________________________________________ + +#ifndef _DM_RES_XSEC_FAST_H_ +#define _DM_RES_XSEC_FAST_H_ + +#include "Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.h" + +namespace genie { + +class DMRESXSecFast : public DMRESXSecWithCacheFast { + +public: + DMRESXSecFast(); + DMRESXSecFast(string param_set); + virtual ~DMRESXSecFast(); + + // 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); + + bool fUsePauliBlocking; ///< account for Pauli blocking? +}; + +} // genie namespace +#endif // _DM_RES_XSEC_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.cxx b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.cxx new file mode 100644 index 0000000000..e5e3c0f540 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.cxx @@ -0,0 +1,314 @@ +//____________________________________________________________________________ +/* + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + + Igor Kakorin + Joint Institute for Nuclear Research + + based on code of + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory +*/ +//____________________________________________________________________________ + +#include +#include + +#include +#include +#include + +#include "Framework/EventGen/XSecAlgorithmI.h" +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Conventions/GBuild.h" +#include "Framework/Conventions/Controls.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Conventions/Units.h" +#include "Framework/Conventions/KinePhaseSpace.h" +#include "Framework/Conventions/KineVar.h" +#include "Physics/XSectionIntegration/GSLXSecFunc.h" +#include "Framework/Interaction/Interaction.h" +#include "Framework/Messenger/Messenger.h" +#include "Framework/ParticleData/PDGUtils.h" +#include "Framework/ParticleData/PDGCodes.h" +#include "Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.h" +#include "Framework/Numerical/MathUtils.h" +#include "Framework/Utils/KineUtils.h" +#include "Framework/Utils/Cache.h" +#include "Framework/Utils/CacheBranchFx.h" +#include "Framework/Numerical/GSLUtils.h" +#include "Framework/Utils/Range1.h" + + + + + + + +using std::ostringstream; + +using namespace genie; +using namespace genie::controls; +using namespace genie::constants; +//using namespace genie::units; + +//____________________________________________________________________________ +DMRESXSecWithCacheFast::DMRESXSecWithCacheFast() : +XSecIntegratorI() +{ + +} +//____________________________________________________________________________ +DMRESXSecWithCacheFast::DMRESXSecWithCacheFast(string nm) : +XSecIntegratorI(nm) +{ + +} +//____________________________________________________________________________ +DMRESXSecWithCacheFast::DMRESXSecWithCacheFast(string nm,string conf): +XSecIntegratorI(nm,conf) +{ + +} +//____________________________________________________________________________ +DMRESXSecWithCacheFast::~DMRESXSecWithCacheFast() +{ + +} +//____________________________________________________________________________ +void DMRESXSecWithCacheFast::CacheResExcitationXSec( + const Interaction * in) const +{ +// Cache resonance neutrino production data from free nucleons + + Cache * cache = Cache::Instance(); + + assert(fSingleResXSecModel); +// assert(fIntegrator); + + // Compute the number of spline knots - use at least 10 knots per decade + // && at least 40 knots in the full energy range + const double Emin = 0.01; + const int nknots_min = (int) (10*(TMath::Log(fEMax)-TMath::Log(Emin))); + const int nknots = TMath::Max(100, nknots_min); + double * E = new double[nknots]; // knot 'x' + + TLorentzVector p4(0,0,0,0); + //*************************************** + //Get Dark Matter mass for calculating p4 + double mchi = in->InitState().GetProbeP4(kRfHitNucRest)->M(); //DM mass: mx + //*************************************** + int nu_code = in->InitState().ProbePdg(); + int nuc_code = in->InitState().Tgt().HitNucPdg(); + int tgt_code = (nuc_code==kPdgProton) ? kPdgTgtFreeP : kPdgTgtFreeN; + + Interaction * interaction = new Interaction(*in); + interaction->InitStatePtr()->SetPdgs(tgt_code, nu_code); + interaction->InitStatePtr()->TgtPtr()->SetHitNucPdg(nuc_code); + + InteractionType_t wkcur = interaction->ProcInfo().InteractionTypeId(); + unsigned int nres = fResList.NResonances(); + for(unsigned int ires = 0; ires < nres; ires++) { + + // Get next resonance from the resonance list + Resonance_t res = fResList.ResonanceId(ires); + + interaction->ExclTagPtr()->SetResonance(res); + + // Get a unique cache branch name + string key = this->CacheBranchName(res, wkcur, nu_code, nuc_code); + + // Make sure the cache branch does not already exists + CacheBranchFx * cache_branch = + dynamic_cast (cache->FindCacheBranch(key)); + assert(!cache_branch); + + // Create the new cache branch + LOG("DMResCF", pNOTICE) + << "\n ** Creating cache branch - key = " << key; + cache_branch = new CacheBranchFx("RES Excitation XSec"); + cache->AddCacheBranch(key, cache_branch); + assert(cache_branch); + + const KPhaseSpace & kps = interaction->PhaseSpace(); + double Ethr = kps.Threshold(); + LOG("DMResCF", pNOTICE) << "E threshold = " << Ethr; + + // Distribute the knots in the energy range as is being done in the + // XSecSplineList so that the energy threshold is treated correctly + // in the spline - see comments there in. + int nkb = (Ethr>Emin) ? 5 : 0; // number of knots < threshold + int nka = nknots-nkb; // number of knots >= threshold + // knots < energy threshold + double dEb = (Ethr>Emin) ? (Ethr - Emin) / nkb : 0; + for(int i=0; i= energy threshold + double E0 = TMath::Max(Ethr,Emin); + double dEa = (TMath::Log10(fEMax) - TMath::Log10(E0)) /(nka-1); + for(int i=0; iInitStatePtr()->SetProbeP4(p4); + + if(Ev>Ethr+kASmallNum) { + // Get integration ranges + Range1D_t rW = Range1D_t(0.0,1.0); + Range1D_t rQ2 = Range1D_t(0.0,1.0); + + LOG("DMResCF", pINFO) + << "*** Integrating d^2 XSec/dWdQ^2 for R: " + << utils::res::AsString(res) << " at Ev = " << Ev; + LOG("DMResCF", pINFO) + << "{W} = " << rW.min << ", " << rW.max; + LOG("DMResCF", pINFO) + << "{Q^2} = " << rQ2.min << ", " << rQ2.max; + + if(rW.maxCreateSpline(); + }//ires + + delete [] E; + delete interaction; +} +//____________________________________________________________________________ +string DMRESXSecWithCacheFast::CacheBranchName( + Resonance_t res, InteractionType_t it, int nupdgc, int nucleonpdgc) const +{ +// Build a unique name for the cache branch + + Cache * cache = Cache::Instance(); + string res_name = utils::res::AsString(res); + string it_name = InteractionType::AsString(it); + string nc_nuc = ((nucleonpdgc==kPdgProton) ? "p" : "n"); + + ostringstream intk; + intk << "ResExcitationXSec/R:" << res_name << ";nu:" << nupdgc + << ";int:" << it_name << nc_nuc; + + string algkey = fSingleResXSecModel->Id().Key(); + string ikey = intk.str(); + string key = cache->CacheBranchKey(algkey, ikey); + + return key; +} +//____________________________________________________________________________ +// GSL wrappers +//____________________________________________________________________________ +genie::utils::gsl::d2XSecRESFast_dWQ2_E::d2XSecRESFast_dWQ2_E( + const XSecAlgorithmI * m, const Interaction * i) : +ROOT::Math::IBaseFunctionMultiDim(), +fModel(m), +fInteraction(i) +{ + kps = fInteraction->PhaseSpacePtr(); + Range1D_t Wl = kps->WLim(); + fWmin = Wl.min; + fWmax = Wl.max; + Registry fConfig = (const_cast(fModel))->GetConfig(); + bool fUsingDisResJoin = fConfig.GetBool("UseDRJoinScheme"); + double fWcut = 999999; + if(fUsingDisResJoin) + { + fWcut = fConfig.GetDouble("Wcut"); + } + fWmax=TMath::Min(fWcut, fWmax); + if (fWcutExclTag().Resonance(); + int IR = utils::res::ResonanceIndex (resonance); + double MR = utils::res::Mass (resonance); + double WR = utils::res::Width (resonance); + if (IR==0) + fWcut = MR + fN0ResMaxNWidths * WR; + else if (IR==2) + fWcut = MR + fN2ResMaxNWidths * WR; + else + fWcut = MR + fGNResMaxNWidths * WR; + fWmax=TMath::Min(fWcut, fWmax); + if (fWcutKinePtr()->SetW(W); + Range1D_t Q2l = kps->Q2Lim_W(); + if (Q2l.min<0 || Q2l.max<0) + return 0.0; + double Q2 = Q2l.min+(Q2l.max-Q2l.min)*xin[1]; + fInteraction->KinePtr()->SetQ2(Q2); + double xsec = fModel->XSec(fInteraction, kPSWQ2fE)*(fWmax-fWmin)*(Q2l.max-Q2l.min); + return xsec/(1E-38 * units::cm2); +} +ROOT::Math::IBaseFunctionMultiDim * + genie::utils::gsl::d2XSecRESFast_dWQ2_E::Clone() const +{ + return + new genie::utils::gsl::d2XSecRESFast_dWQ2_E(fModel,fInteraction); +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.h b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.h new file mode 100644 index 0000000000..23845074d5 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/DMRESXSecWithCacheFast.h @@ -0,0 +1,106 @@ +//____________________________________________________________________________ +/*! + +\class genie::DMRESXSecWithCacheFast + +\brief Class that caches resonance neutrinoproduction cross sections on free + nucleons according to the Rein-Sehgal model. This significantly speeds + the cross section calculation for multiple nuclear targets (eg at the + spline construction phase). This class integrates cross sections faster, + than DMRESXSecWithCache because of integration area transformation. + +\ref D.Rein and L.M.Sehgal, Neutrino Excitation of Baryon Resonances + and Single Pion Production, Ann.Phys.133, 79 (1981) + +\author Igor Kakorin + Joint Institute for Nuclear Research + + based on code of + Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + + adjusted for Dark Matter by + Zachary W. Orr + Colorado State University + +\created March 01, 2017 + +\cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _DM_RES_XSEC_WITH_CACHE_FAST_H_ +#define _DM_RES_XSEC_WITH_CACHE_FAST_H_ + +#include +#include + +#include "Physics/XSectionIntegration/XSecIntegratorI.h" +#include "Framework/ParticleData//BaryonResList.h" +#include "Framework/ParticleData/BaryonResonance.h" +#include "Framework/Utils/Range1.h" +#include "Framework/Interaction/KPhaseSpace.h" + +namespace genie { + +class DMRESXSecWithCacheFast : public XSecIntegratorI { + +protected: + DMRESXSecWithCacheFast(); + DMRESXSecWithCacheFast(string name); + DMRESXSecWithCacheFast(string name, string config); + virtual ~DMRESXSecWithCacheFast(); + + // Don't implement the XSecIntegratorI interface - leave it for the concrete + // subclasses. Just define utility methods and data + void CacheResExcitationXSec (const Interaction * interaction) const; + string CacheBranchName(Resonance_t r, InteractionType_t it, int nu, int nuc) const; + + bool fUsingDisResJoin; + double fWcut; + double fEMax; + + mutable const XSecAlgorithmI * fSingleResXSecModel; + BaryonResList fResList; +}; + + +class XSecAlgorithmI; +class Interaction; + +namespace utils { +namespace gsl { + + + +//..................................................................................... +// +// genie::utils::gsl::d2XSecRESFast_dWQ2_E +// A 2-D cross section function: d2xsec/dWdQ2 = f(W,Q2)|(fixed E) +// +class d2XSecRESFast_dWQ2_E: public ROOT::Math::IBaseFunctionMultiDim +{ +public: + d2XSecRESFast_dWQ2_E(const XSecAlgorithmI * m, const Interaction * i); + ~d2XSecRESFast_dWQ2_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 fWmin; + double fWmax; + bool isfWcutLessfWmin; + KPhaseSpace * kps; +}; + +} // gsl namespace +} // utils namespace +} // genie namespace + +#endif // _DM_RES_XSEC_WITH_CACHE_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/FKRDM.cxx b/src/Physics/BoostedDarkMatter/XSection/FKRDM.cxx new file mode 100644 index 0000000000..f107cb64d0 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/FKRDM.cxx @@ -0,0 +1,70 @@ +//____________________________________________________________________________ +/* +Dark Matter model parameters calculated in the same form as +the Feynmann-Kislinger-Ravndall (FKR) baryon excitation model parameters. + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + +#include + +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Messenger/Messenger.h" +#include "Physics/BoostedDarkMatter/XSection/FKRDM.h" + +using std::endl; + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +namespace genie +{ + ostream & operator << (ostream & stream, const FKRDM & parameters) + { + parameters.Print(stream); + return stream; + } +} +//____________________________________________________________________________ +FKRDM::FKRDM() +{ + this->Reset(); +} +//____________________________________________________________________________ +FKRDM::~FKRDM() +{ + +} +//____________________________________________________________________________ +void FKRDM::Print(ostream & stream) const +{ + stream << endl; + stream << " lamda = " << Lamda << endl; + stream << " Tv = " << Tv << endl; + stream << " Rv = " << Rv << endl; + stream << " S = " << S << endl; + stream << " Ta = " << Ta << endl; + stream << " Ra = " << Ra << endl; + stream << " Bs = " << Bs << endl; + stream << " Cs = " << Cs << endl; + stream << " Bz = " << Bz << endl; + stream << " Cz = " << Cz << endl; +} +//____________________________________________________________________________ +void FKRDM::Reset(void) +{ + Lamda = 0.0; + Tv = 0.0; + Rv = 0.0; + S = 0.0; + Ta = 0.0; + Ra = 0.0; + Bs = 0.0; + Cs = 0.0; + Bz = 0.0; + Cz = 0.0; +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/FKRDM.h b/src/Physics/BoostedDarkMatter/XSection/FKRDM.h new file mode 100644 index 0000000000..7e661199b7 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/FKRDM.h @@ -0,0 +1,51 @@ +//____________________________________________________________________________ +/* +Dark Matter model parameters calculated in the same form as +the Feynmann-Kislinger-Ravndall (FKR) baryon excitation model parameters. (header) + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + +#ifndef _FKRDM_H_ +#define _FKRDM_H_ + +#include + +using std::ostream; + +namespace genie { + +class FKRDM; +ostream & operator<< (ostream & stream, const FKRDM & parameters); + +class FKRDM { + +public: + + friend ostream & operator<< (ostream & stream, const FKRDM & parameters); + + //Unchanged variables + double Lamda; + double Tv; + double Rv; + double S; + double Ta; + double Ra; + + //New variables for DM case that differ from the FKR calculation + double Bs; + double Cs; + double Bz; + double Cz; + + void Reset (void); + void Print (ostream & stream) const; + + FKRDM(); + ~FKRDM(); +}; + +} // genie namespace + +#endif // _FKRDM_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/LinkDef.h b/src/Physics/BoostedDarkMatter/XSection/LinkDef.h index 35371cf40d..48c8ba854a 100644 --- a/src/Physics/BoostedDarkMatter/XSection/LinkDef.h +++ b/src/Physics/BoostedDarkMatter/XSection/LinkDef.h @@ -14,4 +14,14 @@ #pragma link C++ class genie::DMDISXSec; #pragma link C++ class genie::DMElectronXSec; +#pragma link C++ class genie::DMRESPXSec; +#pragma link C++ class genie::DMRESXSec; +#pragma link C++ class genie::DMRESXSecFast; +#pragma link C++ class genie::DMRESXSecWithCacheFast; +#pragma link C++ class genie::FKRDM; +#pragma link C++ class genie::RSHelicityAmplModelDMn; +#pragma link C++ class genie::RSHelicityAmplModelDMp; +#pragma link C++ class genie::RSHelicityAmplModelDMI; +#pragma link C++ class genie::RSHelicityAmplDM; + #endif diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.cxx b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.cxx new file mode 100644 index 0000000000..1e208f416c --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.cxx @@ -0,0 +1,68 @@ +//____________________________________________________________________________ +/* +Dark Matter Calculations of |f|^2 terms for the cross section, based on the +procedure of FKR. + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h" + +using namespace genie; +using std::endl; + +//____________________________________________________________________________ +namespace genie { + ostream & operator<< (ostream & stream, const RSHelicityAmplDM & hamp) + { + hamp.Print(stream); + return stream; + } +} +//____________________________________________________________________________ +RSHelicityAmplDM::RSHelicityAmplDM() +{ + this->Init(); +} +//____________________________________________________________________________ +RSHelicityAmplDM::RSHelicityAmplDM(const RSHelicityAmplDM & hamp) +{ + fMinus1 = hamp.AmpMinus1(); + fPlus1 = hamp.AmpPlus1(); + fMinus3 = hamp.AmpMinus3(); + fPlus3 = hamp.AmpPlus3(); + f0Minus = hamp.Amp0Minus(); + f0Plus = hamp.Amp0Plus(); + //New terms specific to DM + fz0Minus = hamp.Ampz0Minus(); + fz0Plus = hamp.Ampz0Plus(); +} +//____________________________________________________________________________ +void RSHelicityAmplDM::Print(ostream & stream) const +{ + stream << endl; + stream << " f(-1) = " << fMinus1 << endl; + stream << " f(+1) = " << fPlus1 << endl; + stream << " f(-3) = " << fMinus3 << endl; + stream << " f(+3) = " << fPlus3 << endl; + stream << " f(0-) = " << f0Minus << endl; + stream << " f(0+) = " << f0Plus << endl; + //New terms specific to DM + stream << " fz(0-) = " << fz0Minus << endl; + stream << " fz(0+) = " << fz0Plus << endl; +} +//____________________________________________________________________________ +void RSHelicityAmplDM::Init(void) +{ + fMinus1 = 0.0; + fPlus1 = 0.0; + fMinus3 = 0.0; + fPlus3 = 0.0; + f0Minus = 0.0; + f0Plus = 0.0; + //New terms specific to DM + fz0Minus = 0.0; + fz0Plus = 0.0; +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h new file mode 100644 index 0000000000..6fd34b220d --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h @@ -0,0 +1,106 @@ +//____________________________________________________________________________ +/* +Dark Matter Calculations of |f|^2 terms for the cross section, based on the +procedure of FKR. (header) + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + + + + +//____________________________________________________________________________ +/*! + +\class genie::RSHelicityAmpl + +\brief A class holding the Rein-Sehgal's helicity amplitudes. + + This class is using the \b Strategy Pattern. \n + It can accept requests to calculate itself, for a given interaction, + that it then delegates to the algorithmic object, implementing the + RSHelicityAmplModelI interface, that it finds attached to itself. + +\author Costas Andreopoulos + University of Liverpool & STFC Rutherford Appleton Laboratory + +\created May 03, 2004 + +\cpright Copyright (c) 2003-2020, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org +*/ +//____________________________________________________________________________ + +#ifndef _RS_HELICITY_AMPL_DM_H_ +#define _RS_HELICITY_AMPL_DM_H_ + +#include + +#include + +#include "Framework/Interaction/Interaction.h" +#include "Physics/BoostedDarkMatter/XSection/FKRDM.h" + +using std::ostream; + +namespace genie { + +class RSHelicityAmplDM; +ostream & operator<< (ostream & stream, const RSHelicityAmplDM & hamp); + +class RSHelicityAmplDM { +//Only two possibilites +friend class RSHelicityAmplModelDMp; +friend class RSHelicityAmplModelDMn; + +public: + + RSHelicityAmplDM(); + RSHelicityAmplDM(const RSHelicityAmplDM & hamp); + ~RSHelicityAmplDM() { } + + //! return helicity amplitude + double AmpMinus1 (void) const { return fMinus1; } /* f(-1) */ + double AmpPlus1 (void) const { return fPlus1; } /* f(+1) */ + double AmpMinus3 (void) const { return fMinus3; } /* f(-3) */ + double AmpPlus3 (void) const { return fPlus3; } /* f(+3) */ + double Amp0Minus (void) const { return f0Minus; } /* f(0-) */ + double Amp0Plus (void) const { return f0Plus; } /* f(0+) */ + //New terms specific to DM + double Ampz0Minus (void) const { return fz0Minus; } /* fz(0-) */ + double Ampz0Plus (void) const { return fz0Plus; } /* fz(0+) */ + + //! return |helicity amplitude|^2 + double Amp2Minus1 (void) const { return TMath::Power(fMinus1, 2.); } /* |f(-1)|^2 */ + double Amp2Plus1 (void) const { return TMath::Power(fPlus1, 2.); } /* |f(+1)|^2 */ + double Amp2Minus3 (void) const { return TMath::Power(fMinus3, 2.); } /* |f(-3)|^2 */ + double Amp2Plus3 (void) const { return TMath::Power(fPlus3, 2.); } /* |f(+3)|^2 */ + double Amp20Minus (void) const { return TMath::Power(f0Minus, 2.); } /* |f(0-)|^2 */ + double Amp20Plus (void) const { return TMath::Power(f0Plus, 2.); } /* |f(0+)|^2 */ + //New terms specific to DM + double Ampz20Minus (void) const { return TMath::Power(fz0Minus, 2.); } /* |fz(0-)|^2 */ + double Ampz20Plus (void) const { return TMath::Power(fz0Plus, 2.); } /* |fz(0+)|^2 */ + + friend ostream & operator<< (ostream & stream, const RSHelicityAmplDM & hamp); + + void Print(ostream & stream) const; + +private: + + void Init(void); + + double fMinus1; + double fPlus1; + double fMinus3; + double fPlus3; + double f0Minus; + double f0Plus; + //New terms specific to DM + double fz0Minus; + double fz0Plus; +}; + +} // genie namespace + +#endif // _RS_HELICITY_AMPL_DM_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.cxx b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.cxx new file mode 100644 index 0000000000..0e9492e125 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.cxx @@ -0,0 +1,59 @@ +//____________________________________________________________________________ +/* +RSHelicityAmplModelDMI: the interface for the Dark Matter helicity Amplitudes + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h" //Specific to DM + +using namespace genie; + +//____________________________________________________________________________ +RSHelicityAmplModelDMI::RSHelicityAmplModelDMI() : +Algorithm() +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMI::RSHelicityAmplModelDMI(string name) : +Algorithm(name) +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMI::RSHelicityAmplModelDMI(string name, string config) : +Algorithm(name, config) +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMI::~RSHelicityAmplModelDMI() +{ + +} +//____________________________________________________________________________ + +//___________added for quark charges +//____________________________________________________________________________ +void RSHelicityAmplModelDMI::Configure(string config) +{ + Algorithm::Configure(config); + this->LoadConfig(); +} +//____________________________________________________________________________ +void RSHelicityAmplModelDMI::LoadConfig(void) +{ + + // quark couplings to mediator + double QuL, QuR, QdL, QdR; + this->GetParam( "UpLeftCharge", QuL ) ; + this->GetParam( "UpRightCharge", QuR ) ; + this->GetParam( "DownLeftCharge", QdL ) ; + this->GetParam( "DownRightCharge", QdR ) ; + fQuV = 0.5*(QuL + QuR); + fQuA = 0.5*(- QuL + QuR); + fQdV = 0.5*(QdL + QdR); + fQdA = 0.5*(- QdL + QdR); +} +//_____________ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h new file mode 100644 index 0000000000..c7c2c82c2d --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h @@ -0,0 +1,50 @@ +//____________________________________________________________________________ +/* +RSHelicityAmplModelDMI: the interface for the Dark Matter helicity Amplitudes (header) + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + +#ifndef _REIN_SEHGAL_HELICITY_AMPL_MODEL_DM_I_H_ +#define _REIN_SEHGAL_HELICITY_AMPL_MODEL_DM_I_H_ + +#include "Framework/Algorithm/Algorithm.h" +#include "Framework/ParticleData/BaryonResonance.h" +#include "Physics/BoostedDarkMatter/XSection/FKRDM.h" +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h" + +namespace genie { + +class RSHelicityAmplModelDMI : public Algorithm +{ +public: +//______________________ added for quark charges + void Configure(string config); +//______________________ + + virtual ~RSHelicityAmplModelDMI(); + + virtual const RSHelicityAmplDM & Compute(Resonance_t res, const FKRDM & fkrdm) const = 0; + + + +protected: + RSHelicityAmplModelDMI(); + RSHelicityAmplModelDMI(string name); + RSHelicityAmplModelDMI(string name, string config); + + //______________added for quark charges + void LoadConfig(void); + + double fQuV; + double fQuA; + double fQdV; + double fQdA; + //________________ + +}; + +} // namespace + +#endif // _REIN_SEHGAL_HELICITY_AMPL_MODEL_DM_I_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.cxx b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.cxx new file mode 100644 index 0000000000..f351a4c8df --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.cxx @@ -0,0 +1,312 @@ +//____________________________________________________________________________ +/* +Dark Matter Resonant Scattering Helicity Amplitudes for neutron Interactions + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + +#include "Framework/Algorithm/AlgConfigPool.h" +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Messenger/Messenger.h" +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.h" +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplDM.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +RSHelicityAmplModelDMn::RSHelicityAmplModelDMn() : +RSHelicityAmplModelDMI("genie::RSHelicityAmplModelDMn") +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMn::RSHelicityAmplModelDMn(string config) : +RSHelicityAmplModelDMI("genie::RSHelicityAmplModelDMn", config) +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMn::~RSHelicityAmplModelDMn() +{ + +} +//____________________________________________________________________________ +const RSHelicityAmplDM & + RSHelicityAmplModelDMn::Compute( + Resonance_t res, const FKRDM & fkrdm) const +{ + switch(res) { + + case (kP33_1232) : + { + fAmpl.fMinus3 = kSqrt6 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = kSqrt2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus1 = kSqrt2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus3 = kSqrt6 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.f0Plus = 2.0 * kSqrt2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.f0Minus = 2.0 * kSqrt2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.fz0Plus = 2.0 * kSqrt2 * fkrdm.Cz * (fQuA - fQdA); + fAmpl.fz0Minus = 2.0 * kSqrt2 * fkrdm.Cz * (fQuA - fQdA); + + break; + } + case (kS11_1535) : + { + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQuA + 5 * fQdA) - fkrdm.Rv * (fQuV + 5 * fQdV)) + kSqrt3 * fkrdm.Ta * (fQdA - fQuA) + kSqrt3 * fkrdm.Tv * (fQuV - fQdV); + fAmpl.fPlus1 = -k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQuA + 5 * fQdA) + fkrdm.Rv * (fQuV + 5 * fQdV)) + kSqrt3 * fkrdm.Ta * (fQuA - fQdA) + kSqrt3 * fkrdm.Tv * (fQuV - fQdV); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV - fQuV) - 3 * fkrdm.Bs * (fQuA + 5 * fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQuA + 5 * fQdA)); + fAmpl.f0Minus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV - fQdV) - 3 * fkrdm.Bs * (fQuA + 5 * fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQuA + 5 * fQdA)); + fAmpl.fz0Plus = -k1_Sqrt6 * ((fQuA + 5 * fQdA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz)); + fAmpl.fz0Minus = -k1_Sqrt6 * ((fQuA + 5 * fQdA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz)); + + break; + } + case (kD13_1520) : + { + fAmpl.fMinus3 = k3_Sqrt2 * (fkrdm.Tv * (fQuV - fQdV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.fMinus1 = 0.5 * k1_Sqrt3 * (2.0 * fkrdm.Lamda * (fkrdm.Rv * (5 * fQdV + fQuV) - fkrdm.Ra * (5 * fQdA + fQuA)) + 3 * kSqrt2 * (fkrdm.Tv * (fQuV - fQdV) + fkrdm.Ta * (fQdA - fQuA))); + fAmpl.fPlus1 = 0.5 * k1_Sqrt3 * (-2.0 * fkrdm.Lamda * (fkrdm.Ra * fQuA + fkrdm.Rv * (5 * fQdV + fQuV)) + fQdA * (3 * kSqrt2 * fkrdm.Ta - 10 * fkrdm.Lamda * fkrdm.Ra) + 3 * kSqrt2 * fkrdm.Tv * (fQdV - fQuV) - 3 * kSqrt2 * fQuA * fkrdm.Ta); + fAmpl.fPlus3 = k3_Sqrt2 * (fkrdm.Tv * (fQdV - fQuV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.f0Plus = -fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQdV - fQuV) + fkrdm.Cs * (fQuA + 5 * fQdA)); + fAmpl.f0Minus = fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQuV - fQdV) + fkrdm.Cs * (fQuA + 5 * fQdA)); + fAmpl.fz0Plus = -fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQuA + 5 * fQdA); + fAmpl.fz0Minus = fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQuA + 5 * fQdA); + + break; + } + case (kS11_1650) : + { + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt6 * fkrdm.Lamda *(fkrdm.Rv * (2.0 * fQuV + fQdV) - fkrdm.Ra * (2.0 * fQuA + fQdA)); + fAmpl.fPlus1 = k1_Sqrt6* fkrdm.Lamda *(fkrdm.Rv * (2.0 * fQuV + fQdV) + fkrdm.Ra * (2.0 * fQuA + fQdA)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = -kSqrt2_3 * (2.0 * fQuA + fQdA) * (3 * fkrdm.Bs - fkrdm.Lamda * fkrdm.Cs); + fAmpl.f0Minus = -kSqrt2_3 * (2.0 * fQuA + fQdA) * (3 * fkrdm.Bs - fkrdm.Lamda * fkrdm.Cs); + fAmpl.fz0Plus = -kSqrt2_3 * (2.0 * fQuA + fQdA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + fAmpl.fz0Minus = -kSqrt2_3 * (2.0 * fQuA + fQdA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + + break; + } + case (kD13_1700) : + { + fAmpl.fMinus3 = (1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQuV + fQdV) - fkrdm.Ra * (2.0 * fQuA + fQdA)); + fAmpl.fMinus1 = k1_Sqrt30 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQuV + fQdV) - fkrdm.Ra * (2.0 * fQuA + fQdA)); + fAmpl.fPlus1 = -k1_Sqrt30 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQuV + fQdV) + fkrdm.Ra * (2.0 * fQuA + fQdA)); + fAmpl.fPlus3 = -(1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQuV + fQdV) + fkrdm.Ra * (2.0 * fQuA + fQdA)); + fAmpl.f0Plus = (kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQuA + fQdA); + fAmpl.f0Minus = -(kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQuA + fQdA); + fAmpl.fz0Plus = (kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQuA + fQdA); + fAmpl.fz0Minus = -(kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQuA + fQdA); + + break; + } + case (kD15_1675) : + { + fAmpl.fMinus3 = kSqrt3_5 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQuA + fQdA) - fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.fMinus1 = kSqrt3_10 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQuA + fQdA) - fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.fPlus1 = -kSqrt3_10 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQuA + fQdA) + fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.fPlus3 = -kSqrt3_5 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQuA + fQdA) + fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.f0Plus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQuA + fQdA); + fAmpl.f0Minus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQuA + fQdA); + fAmpl.fz0Plus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQuA + fQdA); + fAmpl.fz0Minus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQuA + fQdA); + + break; + } + case (kS31_1620) : + { + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)) + kSqrt3 * fkrdm.Ta * (fQuA - fQdA) + kSqrt3 * fkrdm.Tv * (fQdV - fQuV); + fAmpl.fPlus1 = k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)) + kSqrt3 * fkrdm.Ta * (fQdA - fQuA) + kSqrt3 * fkrdm.Tv * (fQdV - fQuV); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV - fQdV) + 3 * fkrdm.Bs * (fQuA - fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - fQuA)); + fAmpl.f0Minus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV - fQuV) + 3 * fkrdm.Bs * (fQuA - fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - fQuA)); + fAmpl.fz0Plus = -k1_Sqrt6 * (fQdA - fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + fAmpl.fz0Minus = -k1_Sqrt6 * (fQdA - fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + + break; + } + case (kD33_1700) : + { + fAmpl.fMinus3 = -k3_Sqrt2 * (fkrdm.Tv * (fQuV - fQdV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.fMinus1 = 0.5 * k1_Sqrt3 * ((fQuA - fQdA) * (2.0 * fkrdm.Lamda * fkrdm.Ra + 3 * kSqrt2 * fkrdm.Ta) + (fQdV - fQuV) * (2.0 * fkrdm.Lamda * fkrdm.Rv + 3 * kSqrt2 * fkrdm.Tv)); + fAmpl.fPlus1 = 0.5 * k1_Sqrt3 * ((fQuA - fQdA) * (2.0 * fkrdm.Lamda * fkrdm.Ra + 3 * kSqrt2 * fkrdm.Ta) - (fQdV - fQuV) * (2.0 * fkrdm.Lamda * fkrdm.Rv + 3 * kSqrt2 * fkrdm.Tv)); + fAmpl.fPlus3 = -k3_Sqrt2 * (fkrdm.Tv * (fQdV - fQuV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.f0Plus = fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQdV - fQuV) + fkrdm.Cs * (fQuA - fQdA)); + fAmpl.f0Minus = fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQdV - fQuV) + fkrdm.Cs * (fQdA - fQuA)); + fAmpl.fz0Plus = fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQuA - fQdA); + fAmpl.fz0Minus = fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQdA - fQuA); + + break; + } + case (kP11_1440) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = 0.5 * k1_Sqrt3 * L2 * (fkrdm.Ra * (fQuA - 4 * fQdA) - fkrdm.Rv * (fQuV - 4 * fQdV)); + fAmpl.fPlus1 = 0.5 * k1_Sqrt3 * L2 * (fkrdm.Ra * (fQuA - 4 * fQdA) + fkrdm.Rv * (fQuV - 4 * fQdV)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = (fkrdm.Lamda/2.0) * k1_Sqrt3 * (-3 * fkrdm.Lamda * fkrdm.S * (fQuV + 2.0 * fQdV) - 2.0 * fkrdm.Bs * (fQuA - 4 * fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQuA - 4 * fQdA)); + fAmpl.f0Minus = (fkrdm.Lamda/2.0) * k1_Sqrt3 * (-3 * fkrdm.Lamda * fkrdm.S * (fQuV + 2.0 * fQdV) + 2.0 * fkrdm.Bs * (fQuA - 4 * fQdA) + fkrdm.Lamda * fkrdm.Cs * (4 * fQdA - fQuA)); + fAmpl.fz0Plus = (fkrdm.Lamda/2.0) * k1_Sqrt3 * ((fQuA - 4 * fQdA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz)); + fAmpl.fz0Minus = -(fkrdm.Lamda/2.0) * k1_Sqrt3 * ((fQuA - 4 * fQdA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz)); + + break; + } + case (kP33_1600) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = k1_Sqrt2 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fMinus1 = k1_Sqrt6 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus1 = k1_Sqrt6 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus3 = k1_Sqrt2 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.f0Plus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs); + fAmpl.f0Minus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs); + fAmpl.fz0Plus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + fAmpl.fz0Minus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + + break; + } + case (kP13_1720) : + { + fAmpl.fMinus3 = (1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Ta * (fQuA + 2.0 * fQdA) - fkrdm.Tv * (fQuV + 2.0 * fQdV)); + fAmpl.fMinus1 = k1_Sqrt15 * (fkrdm.Lamda/2.0) * (2.0 * fkrdm.Lamda * (fkrdm.Ra * (fQuA - 4 * fQdA) - fkrdm.Rv * (fQuV - 4 * fQdV)) + 9 * kSqrt2 * (fkrdm.Tv * (fQuV + 2.0 * fQdV) - fkrdm.Ta * (fQuA + 2.0 * fQdA))); + fAmpl.fPlus1 = k1_Sqrt15 * (fkrdm.Lamda/2.0) * (9 * kSqrt2 * (fkrdm.Tv * (fQuV + 2.0 * fQdV) + fkrdm.Ta * (fQuA + 2.0 * fQdA)) - 2.0 * fkrdm.Lamda * (fkrdm.Ra * (fQuA - 4 * fQdA) + fkrdm.Rv * (fQuV - 4 * fQdV))); + fAmpl.fPlus3 = -(1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Ta * (fQuA + 2.0 * fQdA) + fkrdm.Tv * (fQuV + 2.0 * fQdV)); + fAmpl.f0Plus = fkrdm.Lamda * k1_Sqrt15 * (-3 * fkrdm.Lamda * fkrdm.S * (fQuV + 2.0 * fQdV) - 5 * fkrdm.Bs * (fQuA - 4 * fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQuA - 4 * fQdA)); + fAmpl.f0Minus = fkrdm.Lamda * k1_Sqrt15 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV + 2.0 * fQdV) - 5 * fkrdm.Bs * (fQuA - 4 * fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQuA - 4 * fQdA)); + fAmpl.fz0Plus = fkrdm.Lamda * k1_Sqrt15 * ((4 * fQdA - fQuA) * (5 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz)); + fAmpl.fz0Minus = fkrdm.Lamda * k1_Sqrt15 * ((4 * fQdA - fQuA) * (5 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz)); + + break; + } + case (kF15_1680) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 3 * fkrdm.Lamda * kSqrt2_5 * (fkrdm.Tv * (fQuV + 2.0 * fQdV) - fkrdm.Ta * (fQuA + 2.0 * fQdA)); + fAmpl.fMinus1 = k1_Sqrt5 * (fkrdm.Lamda/2.0) * (kSqrt2 * fkrdm.Lamda * (fkrdm.Rv * (fQuV - 4 * fQdV) - fkrdm.Ra * (fQuA - 4 * fQdA)) + 6 * (fkrdm.Tv * (fQuV + 2.0 * fQdV) - fkrdm.Ta * (fQuA + 2.0 * fQdA))); + fAmpl.fPlus1 = k1_Sqrt5 * (fkrdm.Lamda/2.0) * (-kSqrt2 * fkrdm.Lamda * (fkrdm.Rv * (fQuV - 4 * fQdV) + fkrdm.Ra * (fQuA - 4 * fQdA)) - 6 * (fkrdm.Tv * (fQuV + 2.0 * fQdV) + fkrdm.Ta * (fQuA + 2.0 * fQdA))); + fAmpl.fPlus3 = -3 * fkrdm.Lamda * kSqrt2_5 * (fkrdm.Tv * (fQuV + 2.0 * fQdV) + fkrdm.Ta * (fQuA + 2.0 * fQdA)); + fAmpl.f0Plus = (1.0/kSqrt10) * L2 * (3 * fkrdm.S * (fQuV + 2.0 * fQdV) - fkrdm.Cs * (fQuA - 4 * fQdA)); + fAmpl.f0Minus = (1.0/kSqrt10) * L2 * (3 * fkrdm.S * (fQuV + 2.0 * fQdV) + fkrdm.Cs * (fQuA - 4 * fQdA)); + fAmpl.fz0Plus = -(1.0/kSqrt10) * L2 * fkrdm.Cz * (fQuA - 4 * fQdA); + fAmpl.fz0Minus = (1.0/kSqrt10) * L2 * fkrdm.Cz * (fQuA - 4 * fQdA); + + break; + } + case (kP31_1910) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = -k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.f0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.fz0Plus = -k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + fAmpl.fz0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + + break; + } + case (kP33_1920) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = k1_Sqrt5 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus3 = k1_Sqrt5 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.f0Plus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.f0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.fz0Plus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + fAmpl.fz0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + + break; + } + case (kF35_1905) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = -3 * L2 * (kSqrt2 * k1_Sqrt35) * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = L2 * k1_Sqrt35 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus1 = L2 * k1_Sqrt35 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus3 = -3 * L2 * (kSqrt2 * k1_Sqrt35) * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.f0Plus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cs * (fQdA - fQuA); + fAmpl.f0Minus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.fz0Plus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cz * (fQdA - fQuA); + fAmpl.fz0Minus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cz * (fQuA - fQdA); + + break; + } + case (kF37_1950) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = kSqrt2_7 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = kSqrt6_35 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus1 = kSqrt6_35 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus3 = kSqrt2_7 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.f0Plus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.f0Minus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.fz0Plus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cz * (fQuA - fQdA); + fAmpl.fz0Minus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cz * (fQuA - fQdA); + + break; + } + case (kP11_1710) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = 0.5 * k1_Sqrt6 * L2 * (fkrdm.Ra * (fQuA + 5 * fQdA) - fkrdm.Rv * (fQuV + 5 * fQdV)); + fAmpl.fPlus1 = 0.5 * k1_Sqrt6 * L2 * (fkrdm.Ra * (fQuA + 5 * fQdA) + fkrdm.Rv * (fQuV + 5 * fQdV)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = fkrdm.Lamda * 0.5 * k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV - fQuV) + (fQuA + 5 * fQdA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs)); + fAmpl.f0Minus = -fkrdm.Lamda * 0.5 * k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV - fQdV) + (fQuA + 5 * fQdA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs)); + fAmpl.fz0Plus = fkrdm.Lamda * 0.5 * k1_Sqrt6 * (fQuA + 5 * fQdA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + fAmpl.fz0Minus = -fkrdm.Lamda * 0.5 * k1_Sqrt6 * (fQuA + 5 * fQdA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + + break; + } + case (kF17_1970) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = (k1_Sqrt2 * k1_Sqrt7) * L2 * (fkrdm.Ra * (2.0 * fQuA + fQdA) - fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.fMinus1 = (kSqrt3_2 * k1_Sqrt35) * L2 * (fkrdm.Ra * (2.0 * fQuA + fQdA) - fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.fPlus1 = -(kSqrt3_2 * k1_Sqrt35) * L2 * (fkrdm.Ra * (2.0 * fQuA + fQdA) + fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.fPlus3 = -(k1_Sqrt2 * k1_Sqrt7) * L2 * (fkrdm.Ra * (2.0 * fQuA + fQdA) + fkrdm.Rv * (2.0 * fQuV + fQdV)); + fAmpl.f0Plus = -kSqrt6_35 * L2 * fkrdm.Cs * (2.0 * fQuA + fQdA); + fAmpl.f0Minus = -kSqrt6_35 * L2 * fkrdm.Cs * (2.0 * fQuA + fQdA); + fAmpl.fz0Plus = -kSqrt6_35 * L2 * fkrdm.Cz * (2.0 * fQuA + fQdA); + fAmpl.fz0Minus = -kSqrt6_35 * L2 * fkrdm.Cz * (2.0 * fQuA + fQdA); + + break; + } + + default: + { + LOG("DMHAmpl", pWARN) << "*** UNRECOGNIZED RESONANCE!"; + fAmpl.fMinus1 = 0.; + fAmpl.fPlus1 = 0.; + fAmpl.fMinus3 = 0.; + fAmpl.fPlus3 = 0.; + fAmpl.f0Minus = 0.; + fAmpl.f0Plus = 0.; + break; + } + + }//switch + + return fAmpl; +} +//____________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.h b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.h new file mode 100644 index 0000000000..fdd1c7f6d4 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMn.h @@ -0,0 +1,31 @@ +//____________________________________________________________________________ +/* +Dark Matter Resonant Scattering Amplitude Calculation for neutron Interactions (header) + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + + +#ifndef _HELICITY_AMPL_MODEL_DM_N_H_ +#define _HELICITY_AMPL_MODEL_DM_N_H_ + +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h" + +namespace genie { + +class RSHelicityAmplModelDMn : public RSHelicityAmplModelDMI { + +public: + RSHelicityAmplModelDMn(); + RSHelicityAmplModelDMn(string config); + virtual ~RSHelicityAmplModelDMn(); + + const RSHelicityAmplDM & Compute(Resonance_t res, const FKRDM & fkrdm) const; + +private: + mutable RSHelicityAmplDM fAmpl; +}; + +} // genie namespace +#endif // _HELICITY_AMPL_MODEL_DM_N_H_ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.cxx b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.cxx new file mode 100644 index 0000000000..09af5615e5 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.cxx @@ -0,0 +1,314 @@ +//____________________________________________________________________________ +/* +Dark Matter Resonant Scattering Helicity Amplitudes for proton Interactions + + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + +#include "Framework/ParticleData/BaryonResUtils.h" +#include "Framework/Conventions/Constants.h" +#include "Framework/Messenger/Messenger.h" +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.h" + +using namespace genie; +using namespace genie::constants; + +//____________________________________________________________________________ +RSHelicityAmplModelDMp::RSHelicityAmplModelDMp() : +RSHelicityAmplModelDMI("genie::RSHelicityAmplModelDMp") +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMp::RSHelicityAmplModelDMp(string config) : +RSHelicityAmplModelDMI("genie::RSHelicityAmplModelDMp", config) +{ + +} +//____________________________________________________________________________ +RSHelicityAmplModelDMp::~RSHelicityAmplModelDMp() +{ + +} +//____________________________________________________________________________ +const RSHelicityAmplDM & + RSHelicityAmplModelDMp::Compute( + Resonance_t res, const FKRDM & fkrdm) const +{ + switch(res) { + + + case (kP33_1232) : + { + fAmpl.fMinus3 = kSqrt6 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = kSqrt2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus1 = kSqrt2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus3 = kSqrt6 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.f0Plus = 2.0 * kSqrt2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.f0Minus = 2.0 * kSqrt2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.fz0Plus = 2.0 * kSqrt2 * fkrdm.Cz * (fQuA - fQdA); + fAmpl.fz0Minus = 2.0 * kSqrt2 * fkrdm.Cz * (fQuA - fQdA); + + break; + } + case (kS11_1535) : + { + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQdA + 5 * fQuA) - fkrdm.Rv * (fQdV + 5 * fQuV)) + kSqrt3 * fkrdm.Ta * (fQuA - fQdA) + kSqrt3 * fkrdm.Tv * (fQdV - fQuV); + fAmpl.fPlus1 = -k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQdA + 5 * fQuA) + fkrdm.Rv * (fQdV + 5 * fQuV)) + kSqrt3 * fkrdm.Ta * (fQdA - fQuA) + kSqrt3 * fkrdm.Tv * (fQdV - fQuV); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV - fQdV) - 3 * fkrdm.Bs * (fQdA + 5 * fQuA) + fkrdm.Lamda * fkrdm.Cs * (fQdA + 5 * fQuA)); + fAmpl.f0Minus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV - fQuV) - 3 * fkrdm.Bs * (fQdA + 5 * fQuA) + fkrdm.Lamda * fkrdm.Cs * (fQdA + 5 * fQuA)); + fAmpl.fz0Plus = -k1_Sqrt6 * ((fQdA + 5 * fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz)); + fAmpl.fz0Minus = -k1_Sqrt6 * ((fQdA + 5 * fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz)); + + break; + } + case (kD13_1520) : + { + fAmpl.fMinus3 = -k3_Sqrt2 * (fkrdm.Tv * (fQuV - fQdV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.fMinus1 = 0.5 * k1_Sqrt3 * (2.0 * fkrdm.Lamda * (fkrdm.Rv * (fQdV + 5 * fQuV) - 5 * fkrdm.Ra * fQuA) - fQdA * (2.0 * fkrdm.Lamda * fkrdm.Ra + 3 * kSqrt2 * fkrdm.Ta) + 3 * kSqrt2 * fkrdm.Tv * (fQdV - fQuV) + 3 * kSqrt2 * fQuA * fkrdm.Ta); + fAmpl.fPlus1 = 0.5 * k1_Sqrt3 * (-2.0 * fkrdm.Lamda * (fkrdm.Rv * (fQdV + 5 * fQuV) + 5 * fkrdm.Ra * fQuA) - fQdA * (2.0 * fkrdm.Lamda * fkrdm.Ra + 3 * kSqrt2 * fkrdm.Ta) + 3 * kSqrt2 * fkrdm.Tv * (fQuV - fQdV) + 3 * kSqrt2 * fQuA * fkrdm.Ta); + fAmpl.fPlus3 = -k3_Sqrt2 * (fkrdm.Tv * (fQdV - fQuV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.f0Plus = -fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQuV - fQdV) + fkrdm.Cs * (fQdA + 5 * fQuA)); + fAmpl.f0Minus = fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQdV - fQuV) + fkrdm.Cs * (fQdA + 5 * fQuA)); + fAmpl.fz0Plus = -fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQdA + 5 * fQuA); + fAmpl.fz0Minus = fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQdA + 5 * fQuA); + + break; + } + case (kS11_1650) : + { + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt6 * fkrdm.Lamda *(fkrdm.Rv * (2.0 * fQdV + fQuV) - fkrdm.Ra * (2.0 * fQdA + fQuA)); + fAmpl.fPlus1 = k1_Sqrt6* fkrdm.Lamda *(fkrdm.Rv * (2.0 * fQdV + fQuV) + fkrdm.Ra * (2.0 * fQdA + fQuA)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = -kSqrt2_3 * (2.0 * fQdA + fQuA) * (3 * fkrdm.Bs - fkrdm.Lamda * fkrdm.Cs); + fAmpl.f0Minus = -kSqrt2_3 * (2.0 * fQdA + fQuA) * (3 * fkrdm.Bs - fkrdm.Lamda * fkrdm.Cs); + fAmpl.fz0Plus = -kSqrt2_3 * (2.0 * fQdA + fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + fAmpl.fz0Minus = -kSqrt2_3 * (2.0 * fQdA + fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + + break; + } + case (kD13_1700) : + { + fAmpl.fMinus3 = (1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQdV + fQuV) - fkrdm.Ra * (2.0 * fQdA + fQuA)); + fAmpl.fMinus1 = k1_Sqrt30 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQdV + fQuV) - fkrdm.Ra * (2.0 * fQdA + fQuA)); + fAmpl.fPlus1 = -k1_Sqrt30 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQdV + fQuV) + fkrdm.Ra * (2.0 * fQdA + fQuA)); + fAmpl.fPlus3 = -(1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Rv * (2.0 * fQdV + fQuV) + fkrdm.Ra * (2.0 * fQdA + fQuA)); + fAmpl.f0Plus = (kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQdA + fQuA); + fAmpl.f0Minus = -(kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQdA + fQuA); + fAmpl.fz0Plus = (kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQdA + fQuA); + fAmpl.fz0Minus = -(kSqrt2 * k1_Sqrt15) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQdA + fQuA); + + break; + } + case (kD15_1675) : + { + fAmpl.fMinus3 = kSqrt3_5 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQdA + fQuA) - fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.fMinus1 = kSqrt3_10 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQdA + fQuA) - fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.fPlus1 = -kSqrt3_10 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQdA + fQuA) + fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.fPlus3 = -kSqrt3_5 * fkrdm.Lamda * (fkrdm.Ra * (2.0 * fQdA + fQuA) + fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.f0Plus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQdA + fQuA); + fAmpl.f0Minus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cs * (2.0 * fQdA + fQuA); + fAmpl.fz0Plus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQdA + fQuA); + fAmpl.fz0Minus = -(kSqrt3 * kSqrt2_5) * fkrdm.Lamda * fkrdm.Cz * (2.0 * fQdA + fQuA); + + break; + } + case (kS31_1620) : + { + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)) + kSqrt3 * fkrdm.Ta * (fQuA - fQdA) + kSqrt3 * fkrdm.Tv * (fQdV - fQuV); + fAmpl.fPlus1 = k1_Sqrt6 * fkrdm.Lamda * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)) + kSqrt3 * fkrdm.Ta * (fQdA - fQuA) + kSqrt3 * fkrdm.Tv * (fQdV - fQuV); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV - fQdV) + 3 * fkrdm.Bs * (fQuA - fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - fQuA)); + fAmpl.f0Minus = k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV - fQuV) + 3 * fkrdm.Bs * (fQuA - fQdA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - fQuA)); + fAmpl.fz0Plus = -k1_Sqrt6 * (fQdA - fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + fAmpl.fz0Minus = -k1_Sqrt6 * (fQdA - fQuA) * (3 * fkrdm.Bz - fkrdm.Lamda * fkrdm.Cz); + + break; + } + case (kD33_1700) : + { + fAmpl.fMinus3 = -k3_Sqrt2 * (fkrdm.Tv * (fQuV - fQdV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.fMinus1 = 0.5 * k1_Sqrt3 * ((fQuA - fQdA) * (2.0 * fkrdm.Lamda * fkrdm.Ra + 3 * kSqrt2 * fkrdm.Ta) + (fQdV - fQuV) * (2.0 * fkrdm.Lamda * fkrdm.Rv + 3 * kSqrt2 * fkrdm.Tv)); + fAmpl.fPlus1 = 0.5 * k1_Sqrt3 * ((fQuA - fQdA) * (2.0 * fkrdm.Lamda * fkrdm.Ra + 3 * kSqrt2 * fkrdm.Ta) - (fQdV - fQuV) * (2.0 * fkrdm.Lamda * fkrdm.Rv + 3 * kSqrt2 * fkrdm.Tv)); + fAmpl.fPlus3 = -k3_Sqrt2 * (fkrdm.Tv * (fQdV - fQuV) + fkrdm.Ta * (fQdA - fQuA)); + fAmpl.f0Plus = fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQdV - fQuV) + fkrdm.Cs * (fQuA - fQdA)); + fAmpl.f0Minus = fkrdm.Lamda * k1_Sqrt3 * (3 * fkrdm.S * (fQdV - fQuV) + fkrdm.Cs * (fQdA - fQuA)); + fAmpl.fz0Plus = fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQuA - fQdA); + fAmpl.fz0Minus = fkrdm.Lamda * fkrdm.Cz * k1_Sqrt3 * (fQdA - fQuA); + + break; + } + case (kP11_1440) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = 0.5 * k1_Sqrt3 * L2 * (fkrdm.Ra * (fQdA - 4 * fQuA) - fkrdm.Rv * (fQdV - 4 * fQuV)); + fAmpl.fPlus1 = 0.5 * k1_Sqrt3 * L2 * (fkrdm.Ra * (fQdA - 4 * fQuA) + fkrdm.Rv * (fQdV - 4 * fQuV)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = (fkrdm.Lamda/2.0) * k1_Sqrt3 * (-3 * fkrdm.Lamda * fkrdm.S * (fQdV + 2.0 * fQuV) - 2.0 * fkrdm.Bs * (fQdA - 4 * fQuA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - 4 * fQuA)); + fAmpl.f0Minus = -(fkrdm.Lamda/2.0) * k1_Sqrt3 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV + 2.0 * fQuV) - 2.0 * fkrdm.Bs * (fQdA - 4 * fQuA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - 4 * fQuA)); + fAmpl.fz0Plus = (fkrdm.Lamda/2.0) * k1_Sqrt3 * ((fQdA - 4 * fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz)); + fAmpl.fz0Minus = -(fkrdm.Lamda/2.0) * k1_Sqrt3 * ((fQdA - 4 * fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz)); + + break; + } + case (kP33_1600) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = k1_Sqrt2 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fMinus1 = k1_Sqrt6 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus1 = k1_Sqrt6 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus3 = k1_Sqrt2 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.f0Plus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs); + fAmpl.f0Minus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs); + fAmpl.fz0Plus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + fAmpl.fz0Minus = kSqrt2_3 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + + break; + } + case (kP13_1720) : + { + fAmpl.fMinus3 = (1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Ta * (fQdA + 2.0 * fQuA) - fkrdm.Tv * (fQdV + 2.0 * fQuV)); + fAmpl.fMinus1 = k1_Sqrt15 * (fkrdm.Lamda/2.0) * (2.0 * fkrdm.Lamda * (fkrdm.Ra * (fQdA - 4 * fQuA) - fkrdm.Rv * (fQdV - 4 * fQuV)) + 9 * kSqrt2 * (fkrdm.Tv * (fQdV + 2.0 * fQuV) - fkrdm.Ta * (fQdA + 2.0 * fQuA))); + fAmpl.fPlus1 = k1_Sqrt15 * (fkrdm.Lamda/2.0) * (9 * kSqrt2 * (fkrdm.Tv * (fQdV + 2.0 * fQuV) + fkrdm.Ta * (fQdA + 2.0 * fQuA)) - 2.0 * fkrdm.Lamda * (fkrdm.Ra * (fQdA - 4 * fQuA) + fkrdm.Rv * (fQdV - 4 * fQuV))); + fAmpl.fPlus3 = -(1.0/kSqrt10) * 3 * fkrdm.Lamda * (fkrdm.Ta * (fQdA + 2.0 * fQuA) + fkrdm.Tv * (fQdV + 2.0 * fQuV)); + fAmpl.f0Plus = fkrdm.Lamda * k1_Sqrt15 * (-3 * fkrdm.Lamda * fkrdm.S * (fQdV + 2.0 * fQuV) - 5 * fkrdm.Bs * (fQdA - 4 * fQuA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - 4 * fQuA)); + fAmpl.f0Minus = fkrdm.Lamda * k1_Sqrt15 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV + 2.0 * fQuV) - 5 * fkrdm.Bs * (fQdA - 4 * fQuA) + fkrdm.Lamda * fkrdm.Cs * (fQdA - 4 * fQuA)); + fAmpl.fz0Plus = fkrdm.Lamda * k1_Sqrt15 * ((fQdA - 4 * fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz)); + fAmpl.fz0Minus = fkrdm.Lamda * k1_Sqrt15 * ((fQdA - 4 * fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz)); + + break; + } + case (kF15_1680) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 3 * fkrdm.Lamda * kSqrt2_5 * (fkrdm.Tv * (fQdV + 2.0 * fQuV) - fkrdm.Ta * (fQdA + 2.0 * fQuA)); + fAmpl.fMinus1 = k1_Sqrt5 * (fkrdm.Lamda/2.0) * (kSqrt2 * fkrdm.Lamda * (fkrdm.Rv * (fQdV - 4 * fQuV) - fkrdm.Ra * (fQdA - 4 * fQuA)) + 6 * (fkrdm.Tv * (fQdV + 2.0 * fQuV) - fkrdm.Ta * (fQdA + 2.0 * fQuA))); + fAmpl.fPlus1 = k1_Sqrt5 * (fkrdm.Lamda/2.0) * (-kSqrt2 * fkrdm.Lamda * (fkrdm.Rv * (fQdV - 4 * fQuV) + fkrdm.Ra * (fQdA - 4 * fQuA)) - 6 * (fkrdm.Tv * (fQdV + 2.0 * fQuV) + fkrdm.Ta * (fQdA + 2.0 * fQuA))); + fAmpl.fPlus3 = -3 * fkrdm.Lamda * kSqrt2_5 * (fkrdm.Tv * (fQdV + 2.0 * fQuV) + fkrdm.Ta * (fQdA + 2.0 * fQuA)); + fAmpl.f0Plus = (1.0/kSqrt10) * L2 * (3 * fkrdm.S * (fQdV + 2.0 * fQuV) - fkrdm.Cs * (fQdA - 4 * fQuA)); + fAmpl.f0Minus = (1.0/kSqrt10) * L2 * (3 * fkrdm.S * (fQdV + 2.0 * fQuV) + fkrdm.Cs * (fQdA - 4 * fQuA)); + fAmpl.fz0Plus = -(1.0/kSqrt10) * L2 * fkrdm.Cz * (fQdA - 4 * fQuA); + fAmpl.fz0Minus = (1.0/kSqrt10) * L2 * fkrdm.Cz * (fQdA - 4 * fQuA); + + break; + } + case (kP31_1910) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = -k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.f0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.fz0Plus = -k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + fAmpl.fz0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + + break; + } + case (kP33_1920) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = k1_Sqrt5 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus1 = k1_Sqrt15 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus3 = k1_Sqrt5 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.f0Plus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.f0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cs - 5 * fkrdm.Bs); + fAmpl.fz0Plus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + fAmpl.fz0Minus = k1_Sqrt15 * 2.0 * fkrdm.Lamda * (fQdA - fQuA) * (fkrdm.Lamda * fkrdm.Cz - 5 * fkrdm.Bz); + + break; + } + case (kF35_1905) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = -3 * L2 * (kSqrt2 * k1_Sqrt35) * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = L2 * k1_Sqrt35 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.fPlus1 = L2 * k1_Sqrt35 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus3 = -3 * L2 * (kSqrt2 * k1_Sqrt35) * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQdV - fQuV)); + fAmpl.f0Plus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cs * (fQdA - fQuA); + fAmpl.f0Minus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.fz0Plus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cz * (fQdA - fQuA); + fAmpl.fz0Minus = 2.0 * L2 * k1_Sqrt35 * fkrdm.Cz * (fQuA - fQdA); + + break; + } + case (kF37_1950) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = kSqrt2_7 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fMinus1 = kSqrt6_35 * L2 * (fkrdm.Ra * (fQdA - fQuA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus1 = kSqrt6_35 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.fPlus3 = kSqrt2_7 * L2 * (fkrdm.Ra * (fQuA - fQdA) + fkrdm.Rv * (fQuV - fQdV)); + fAmpl.f0Plus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.f0Minus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cs * (fQuA - fQdA); + fAmpl.fz0Plus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cz * (fQuA - fQdA); + fAmpl.fz0Minus = 2.0 * kSqrt6_35 * L2 * fkrdm.Cz * (fQuA - fQdA); + + break; + } + case (kP11_1710) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = 0; + fAmpl.fMinus1 = 0.5 * k1_Sqrt6 * L2 * (fkrdm.Ra * (fQdA + 5 * fQuA) - fkrdm.Rv * (fQdV + 5 * fQuV)); + fAmpl.fPlus1 = 0.5 * k1_Sqrt6 * L2 * (fkrdm.Ra * (fQdA + 5 * fQuA) + fkrdm.Rv * (fQdV + 5 * fQuV)); + fAmpl.fPlus3 = 0; + fAmpl.f0Plus = fkrdm.Lamda * 0.5 * k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQuV - fQdV) + (fQdA + 5 * fQuA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs)); + fAmpl.f0Minus = -fkrdm.Lamda * 0.5 * k1_Sqrt6 * (3 * fkrdm.Lamda * fkrdm.S * (fQdV - fQuV) + (fQdA + 5 * fQuA) * (fkrdm.Lamda * fkrdm.Cs - 2.0 * fkrdm.Bs)); + fAmpl.fz0Plus = fkrdm.Lamda * 0.5 * k1_Sqrt6 * (fQdA + 5 * fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + fAmpl.fz0Minus = -fkrdm.Lamda * 0.5 * k1_Sqrt6 * (fQdA + 5 * fQuA) * (fkrdm.Lamda * fkrdm.Cz - 2.0 * fkrdm.Bz); + + break; + } + case (kF17_1970) : + { + double L2 = TMath::Power(fkrdm.Lamda, 2); + + fAmpl.fMinus3 = (k1_Sqrt2 * k1_Sqrt7) * L2 * (fkrdm.Ra * (2.0 * fQdA + fQuA) - fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.fMinus1 = (kSqrt3_2 * k1_Sqrt35) * L2 * (fkrdm.Ra * (2.0 * fQdA + fQuA) - fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.fPlus1 = -(kSqrt3_2 * k1_Sqrt35) * L2 * (fkrdm.Ra * (2.0 * fQdA + fQuA) + fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.fPlus3 = -(k1_Sqrt2 * k1_Sqrt7) * L2 * (fkrdm.Ra * (2.0 * fQdA + fQuA) + fkrdm.Rv * (2.0 * fQdV + fQuV)); + fAmpl.f0Plus = -kSqrt6_35 * L2 * fkrdm.Cs * (2.0 * fQdA + fQuA); + fAmpl.f0Minus = -kSqrt6_35 * L2 * fkrdm.Cs * (2.0 * fQdA + fQuA); + fAmpl.fz0Plus = -kSqrt6_35 * L2 * fkrdm.Cz * (2.0 * fQdA + fQuA); + fAmpl.fz0Minus = -kSqrt6_35 * L2 * fkrdm.Cz * (2.0 * fQdA + fQuA); + + break; + } + default: + { + LOG("DMHAmpl", pWARN) << "*** UNRECOGNIZED RESONANCE!"; + fAmpl.fMinus1 = 0.; + fAmpl.fPlus1 = 0.; + fAmpl.fMinus3 = 0.; + fAmpl.fPlus3 = 0.; + fAmpl.f0Minus = 0.; + fAmpl.f0Plus = 0.; + //New terms specific to DM + fAmpl.fz0Minus = 0.; + fAmpl.fz0Plus = 0.; + break; + } + + }//switch + + return fAmpl; + } + //___________________________________________________________________________ diff --git a/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.h b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.h new file mode 100644 index 0000000000..39127a7f91 --- /dev/null +++ b/src/Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMp.h @@ -0,0 +1,32 @@ +//____________________________________________________________________________ +/* +Dark Matter Resonant Scattering Amplitude Calculation for proton Interactions (header) + +Zachary W. Orr, Colorado State University +*/ +//____________________________________________________________________________ + + +#ifndef _HELICITY_AMPL_MODEL_DM_P_H_ +#define _HELICITY_AMPL_MODEL_DM_P_H_ + +#include "Physics/BoostedDarkMatter/XSection/RSHelicityAmplModelDMI.h" + +namespace genie { + +class RSHelicityAmplModelDMp : public RSHelicityAmplModelDMI { + +public: + RSHelicityAmplModelDMp(); + RSHelicityAmplModelDMp(string config); + virtual ~RSHelicityAmplModelDMp(); + + + const RSHelicityAmplDM & Compute(Resonance_t res, const FKRDM & fkrdm) const; + +private: + mutable RSHelicityAmplDM fAmpl; +}; + +} // genie namespace +#endif // _HELICITY_AMPL_MODEL_DM_N_H_ diff --git a/src/Tools/Flux/GAtmoFlux.cxx b/src/Tools/Flux/GAtmoFlux.cxx index c1809d2455..9e37a69a0a 100644 --- a/src/Tools/Flux/GAtmoFlux.cxx +++ b/src/Tools/Flux/GAtmoFlux.cxx @@ -63,6 +63,11 @@ bool GAtmoFlux::GenerateNext(void) double wght = this->Weight(); bool accept = (E<=Emax && E>=Emin && wght>0); + + LOG("Flux", pDEBUG) << "Neutrino w/ energy : " << E; + LOG("Flux", pDEBUG) << "Max : " << Emax; + LOG("Flux", pDEBUG) << "Min : " << Emin; + LOG("Flux", pDEBUG) << "Weight : " << wght; if(accept) return true; } return false;