diff --git a/src/Apps/gNNBarOscEvGen.cxx b/src/Apps/gNNBarOscEvGen.cxx index 37a896c47..81e81a85b 100644 --- a/src/Apps/gNNBarOscEvGen.cxx +++ b/src/Apps/gNNBarOscEvGen.cxx @@ -1,18 +1,9 @@ -//________________________________________________________________________________________ -/*! - -\program gevgen_nosc - -\brief A GENIE-based n-nbar oscillation event generation application. - - *** Synopsis : - - gevgen_nosc [-h] - [-r run#] +//________________________________________________________________________________________ /*! \program gevgen_nosc \brief A GENIE-based n-nbar oscillation event generation application. *** Synopsis : gevgen_nosc [-h] + [-r run#] -n n_of_events [-m decay_mode] -g geometry - [-L geometry_length_units] + [-L geometry_length_units] [-D geometry_density_units] [-t geometry_top_volume_name] [-o output_event_file_prefix] @@ -25,17 +16,17 @@ [] Denotes an optional argument - -h + -h Prints out the gevgen_nosc syntax and exits. - -r + -r Specifies the MC run number [default: 1000]. - -n + -n Specifies how many events to generate. - -m + -m Nucleon decay mode ID: --------------------------------------------------------- - ID | Decay Mode - | + ID | Decay Mode + | --------------------------------------------------------- 0 | Random decay mode 1 | p + nbar --> \pi^{+} + \pi^{0} @@ -56,31 +47,31 @@ 16 | n + nbar --> 2\pi^{+} + 2\pi^{-} + 2\pi^{0} --------------------------------------------------------- - -g + -g Input 'geometry'. This option can be used to specify any of: 1 > A ROOT file containing a ROOT/GEANT geometry description - [Examples] - - To use the master volume from the ROOT geometry stored + [Examples] + - To use the master volume from the ROOT geometry stored in the laguna-lbno.root file, type: '-g /some/path/laguna-lbno.root' 2 > A mix of target materials, each with its corresponding weight, typed as a comma-separated list of nuclear PDG codes (in the std PDG2006 convention: 10LZZZAAAI) with the weight fractions in brackets, eg code1[fraction1],code2[fraction2],... - If that option is used (no detailed input geometry description) + If that option is used (no detailed input geometry description) then the interaction vertices are distributed in the detector by the detector MC. - [Examples] + [Examples] - To use a target mix of 88.9% O16 and 11.1% Hydrogen type: '-g 1000080160[0.889],1000010010[0.111]' - -L + -L Input geometry length units, eg 'm', 'cm', 'mm', ... [default: 'mm'] - -D - Input geometry density units, eg 'g_cm3', 'clhep_def_density_unit',... + -D + Input geometry density units, eg 'g_cm3', 'clhep_def_density_unit',... [default: 'g_cm3'] - -t + -t Input 'top volume' for event generation. The option be used to force event generation in given sub-detector. [default: the 'master volume' of the input geometry] @@ -97,11 +88,11 @@ except the ones explicitly turned on. Vice versa, if the very first character is a `-', GENIE will keep all volumes except the ones explicitly turned off (feature contributed by J.Holeczek). - -o - Sets the prefix of the output event file. - The output filename is built as: + -o + Sets the prefix of the output event file. + The output filename is built as: [prefix].[run_number].[event_tree_format].[file_format] - The default output filename is: + The default output filename is: gntp.[run_number].ghep.root This cmd line arguments lets you override 'gntp' --seed @@ -109,23 +100,28 @@ \author Costas Andreopoulos University of Liverpool & STFC Rutherford Appleton Laboratory - + \created November 03, 2011 - + +\author Linyan Wan + Fermilab + +\created April 18, 2024 + \cpright Copyright (c) 2003-2023, The GENIE Collaboration For the full text of the license visit http://copyright.genie-mc.org - + */ //_________________________________________________________________________________________ #include #include -#include +#include #include #include -#include +#include #include "Framework/Algorithm/AlgFactory.h" #include "Framework/EventGen/EventRecord.h" @@ -135,6 +131,7 @@ #include "Framework/Messenger/Messenger.h" #include "Framework/Ntuple/NtpWriter.h" #include "Physics/NNBarOscillation/NNBarOscMode.h" +#include "Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.h" #include "Physics/NNBarOscillation/NNBarOscUtils.h" #include "Framework/Numerical/RandomGen.h" #include "Framework/ParticleData/PDGCodes.h" @@ -156,14 +153,14 @@ using namespace genie; // function prototypes void GetCommandLineArgs (int argc, char ** argv); void PrintSyntax (void); -int SelectAnnihilationMode (int pdg_code); +int SelectAnnihilationMode (int pdg_code, const EventRecordVisitorI* mcgen); int SelectInitState (void); const EventRecordVisitorI * NeutronOscGenerator(void); // string kDefOptGeomLUnits = "mm"; // default geometry length units string kDefOptGeomDUnits = "g_cm3"; // default geometry density units -NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // default event tree format +NtpMCFormat_t kDefOptNtpFormat = kNFGHEP; // default event tree format string kDefOptEvFilePrefix = "gntp"; // @@ -174,9 +171,9 @@ string gOptEvFilePrefix = kDefOptEvFilePrefix; // event file prefix bool gOptUsingRootGeom = false; // using root geom or target mix? map gOptTgtMix; // target mix (tgt pdg -> wght frac) / if not using detailed root geom string gOptRootGeom; // input ROOT file with realistic detector geometry -string gOptRootGeomTopVol = ""; // input geometry top event generation volume -double gOptGeomLUnits = 0; // input geometry length units -double gOptGeomDUnits = 0; // input geometry density units +string gOptRootGeomTopVol = ""; // input geometry top event generation volume +double gOptGeomLUnits = 0; // input geometry length units +double gOptGeomDUnits = 0; // input geometry density units long int gOptRanSeed = -1; // random number seed //_________________________________________________________________________________________ @@ -215,11 +212,11 @@ int main(int argc, char ** argv) EventRecord * event = new EventRecord; int target = SelectInitState(); - int decay = SelectAnnihilationMode(target); + int decay = SelectAnnihilationMode(target, mcgen); Interaction * interaction = Interaction::NOsc(target,decay); event->AttachSummary(interaction); - // Simulate decay + // Simulate decay mcgen->ProcessEventRecord(event); LOG("gevgen_nnbar_osc", pINFO) @@ -241,7 +238,7 @@ int main(int argc, char ** argv) return 0; } //_________________________________________________________________________________________ -int SelectAnnihilationMode(int pdg_code) +int SelectAnnihilationMode(int pdg_code, const EventRecordVisitorI* mcgen) { // if the mode is set to 'random' (the default), pick one at random! if (gOptDecayMode == kNORandom) { @@ -263,15 +260,16 @@ int SelectAnnihilationMode(int pdg_code) double proton_frac = ((double)n_protons) / ((double)n_nucleons); double neutron_frac = 1 - proton_frac; - // set branching ratios, taken from bubble chamber data - const int n_modes = 16; - double br [n_modes] = { 0.010, 0.080, 0.100, 0.220, - 0.360, 0.160, 0.070, 0.020, - 0.015, 0.065, 0.110, 0.280, - 0.070, 0.240, 0.100, 0.100 }; + // set branching ratios, taken from SK2021 + const auto* osc_gen = dynamic_cast< + const NNBarOscPrimaryVtxGenerator* >( mcgen ); + const int np_modes = osc_gen->GetNPModes(); + const int nn_modes = osc_gen->GetNNModes(); + const int n_modes = np_modes + nn_modes; + auto br = osc_gen->GetBRs(); for (int i = 0; i < n_modes; i++) { - if (i < 7) + if (i < np_modes) br[i] *= proton_frac; else br[i] *= neutron_frac; @@ -340,7 +338,7 @@ void GetCommandLineArgs(int argc, char ** argv) { LOG("gevgen_nnbar_osc", pINFO) << "Parsing command line arguments"; - // Common run options. + // Common run options. RunOpt::Instance()->ReadFromCommandLine(argc,argv); // Parse run options for this app @@ -366,11 +364,11 @@ void GetCommandLineArgs(int argc, char ** argv) // number of events if( parser.OptionExists('n') ) { - LOG("gevgen_nnbar_osc", pDEBUG) + LOG("gevgen_nnbar_osc", pDEBUG) << "Reading number of events to generate"; gOptNev = parser.ArgAsInt('n'); } else { - LOG("gevgen_nnbar_osc", pFATAL) + LOG("gevgen_nnbar_osc", pFATAL) << "You need to specify the number of events"; PrintSyntax(); exit(0); @@ -379,14 +377,14 @@ void GetCommandLineArgs(int argc, char ** argv) // decay mode int mode = 0; if( parser.OptionExists('m') ) { - LOG("gevgen_nnbar_osc", pDEBUG) + LOG("gevgen_nnbar_osc", pDEBUG) << "Reading annihilation mode"; mode = parser.ArgAsInt('m'); } gOptDecayMode = (NNBarOscMode_t) mode; bool valid_mode = utils::nnbar_osc::IsValidMode(gOptDecayMode); if(!valid_mode) { - LOG("gevgen_nnbar_osc", pFATAL) + LOG("gevgen_nnbar_osc", pFATAL) << "You need to specify a valid annihilation mode"; PrintSyntax(); exit(0); @@ -403,14 +401,14 @@ void GetCommandLineArgs(int argc, char ** argv) geom = parser.ArgAsString('g'); // is it a ROOT file that contains a ROOT geometry? - bool accessible_geom_file = + bool accessible_geom_file = ! (gSystem->AccessPathName(geom.c_str())); if (accessible_geom_file) { gOptRootGeom = geom; gOptUsingRootGeom = true; - } + } } else { - LOG("gevgen_nnbar_osc", pFATAL) + LOG("gevgen_nnbar_osc", pFATAL) << "No geometry option specified - Exiting"; PrintSyntax(); exit(1); @@ -421,7 +419,7 @@ void GetCommandLineArgs(int argc, char ** argv) // legth units: if( parser.OptionExists('L') ) { - LOG("gevgen_nnbar_osc", pDEBUG) + LOG("gevgen_nnbar_osc", pDEBUG) << "Checking for input geometry length units"; lunits = parser.ArgAsString('L'); } else { @@ -430,24 +428,24 @@ void GetCommandLineArgs(int argc, char ** argv) } // -L // density units: if( parser.OptionExists('D') ) { - LOG("gevgen_nnbar_osc", pDEBUG) + LOG("gevgen_nnbar_osc", pDEBUG) << "Checking for input geometry density units"; dunits = parser.ArgAsString('D'); } else { LOG("gevgen_nnbar_osc", pDEBUG) << "Using default geometry density units"; dunits = kDefOptGeomDUnits; - } // -D + } // -D gOptGeomLUnits = utils::units::UnitFromString(lunits); gOptGeomDUnits = utils::units::UnitFromString(dunits); - // check whether an event generation volume name has been + // check whether an event generation volume name has been // specified -- default is the 'top volume' if( parser.OptionExists('t') ) { LOG("gevgen_nnbar_osc", pDEBUG) << "Checking for input volume name"; gOptRootGeomTopVol = parser.ArgAsString('t'); } else { LOG("gevgen_nnbar_osc", pDEBUG) << "Using the "; - } // -t + } // -t } // using root geom? @@ -462,18 +460,18 @@ void GetCommandLineArgs(int argc, char ** argv) if(tgtmix.size()==1) { int pdg = atoi(tgtmix[0].c_str()); double wgt = 1.0; - gOptTgtMix.insert(map::value_type(pdg, wgt)); + gOptTgtMix.insert(map::value_type(pdg, wgt)); } else { vector::const_iterator tgtmix_iter = tgtmix.begin(); for( ; tgtmix_iter != tgtmix.end(); ++tgtmix_iter) { string tgt_with_wgt = *tgtmix_iter; string::size_type open_bracket = tgt_with_wgt.find("["); string::size_type close_bracket = tgt_with_wgt.find("]"); - if (open_bracket ==string::npos || - close_bracket==string::npos) + if (open_bracket ==string::npos || + close_bracket==string::npos) { - LOG("gevgen_nnbar_osc", pFATAL) - << "You made an error in specifying the target mix"; + LOG("gevgen_nnbar_osc", pFATAL) + << "You made an error in specifying the target mix"; PrintSyntax(); exit(1); } @@ -483,7 +481,7 @@ void GetCommandLineArgs(int argc, char ** argv) string::size_type jend = close_bracket; int pdg = atoi(tgt_with_wgt.substr(ibeg,iend-ibeg).c_str()); double wgt = atof(tgt_with_wgt.substr(jbeg,jend-jbeg).c_str()); - LOG("gevgen_nnbar_osc", pDEBUG) + LOG("gevgen_nnbar_osc", pDEBUG) << "Adding to target mix: pdg = " << pdg << ", wgt = " << wgt; gOptTgtMix.insert(map::value_type(pdg, wgt)); @@ -542,7 +540,7 @@ void GetCommandLineArgs(int argc, char ** argv) << "\n\n" << utils::print::PrintFramedMesg("gevgen_nosc job configuration"); - LOG("gevgen_nnbar_osc", pNOTICE) + LOG("gevgen_nnbar_osc", pNOTICE) << "\n @@ Run number: " << gOptRunNu << "\n @@ Random number seed: " << gOptRanSeed << "\n @@ Decay channel $ " << utils::nnbar_osc::AsString(gOptDecayMode) @@ -553,7 +551,7 @@ void GetCommandLineArgs(int argc, char ** argv) // Temporary warnings... // if(gOptUsingRootGeom) { - LOG("gevgen_nnbar_osc", pWARN) + LOG("gevgen_nnbar_osc", pWARN) << "\n ** ROOT geometries not supported yet in neutron oscillation mode" << "\n ** (they will be in the very near future)" << "\n ** Please specify a `target mix' instead."; @@ -564,7 +562,7 @@ void GetCommandLineArgs(int argc, char ** argv) //_________________________________________________________________________________________ void PrintSyntax(void) { - LOG("gevgen_nnbar_osc", pFATAL) + LOG("gevgen_nnbar_osc", pFATAL) << "\n **Syntax**" << "\n gevgen_nnbarosc [-h] " << "\n [-r run#]" diff --git a/src/Physics/NNBarOscillation/NNBarOscMode.h b/src/Physics/NNBarOscillation/NNBarOscMode.h index 1dd512ebc..54ddc2f8a 100644 --- a/src/Physics/NNBarOscillation/NNBarOscMode.h +++ b/src/Physics/NNBarOscillation/NNBarOscMode.h @@ -1,18 +1,24 @@ //____________________________________________________________________________ /*! -\class genie::NNBarOscMode + \class genie::NNBarOscMode -\brief Enumeration of neutron oscillation annihilation modes. + \brief Enumeration of neutron oscillation annihilation modes. -\author Jeremy Hewes, Georgia Karagiorgi - University of Manchester + \author Jeremy Hewes, Georgia Karagiorgi + University of Manchester -\created November, 2016 + \created November, 2016 -\cpright Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org -*/ + \author Linyan Wan + Fermilab + + \updated April, 2024 + + + \cpright Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org + */ //____________________________________________________________________________ #ifndef _N_NBAR_OSC_MODE_H_ @@ -24,30 +30,58 @@ namespace genie { - typedef enum ENNBarOscMode { - - // i just replaced all the nucleon decay modes with nnbar modes -j - - kNONull = -1, - kNORandom, // Will select a random decay mode -j - kNOpto1pip1pi0, // p + nbar --> \pi^{+} + \pi^{0} - kNOpto1pip2pi0, // p + nbar --> \pi^{+} + 2\pi^{0} - kNOpto1pip3pi0, // p + nbar --> \pi^{+} + 3\pi^{0} - kNOpto2pip1pim1pi0, // p + nbar --> 2\pi^{+} + \pi^{-} + \pi^{0} - kNOpto2pip1pim2pi0, // p + nbar --> 2\pi^{+} + \pi^{-} + 2\pi^{0} - kNOpto2pip1pim2o, // p + nbar --> 2\pi^{+} + \pi^{-} + 2\omega^{0} - kNOpto3pip2pim1pi0, // p + nbar --> 3\pi^{+} + 2\pi^{-} + \pi^{0} - kNOnto1pip1pim, // n + nbar --> \pi^{+} + \pi^{-} - kNOnto2pi0, // n + nbar --> 2\pi^{0} - kNOnto1pip1pim1pi0, // n + nbar --> \pi^{+} + \pi^{-} + \pi^{0} - kNOnto1pip1pim2pi0, // n + nbar --> \pi^{+} + \pi^{-} + 2\pi^{0} - kNOnto1pip1pim3pi0, // n + nbar --> \pi^{+} + \pi^{-} + 3\pi^{0} - kNOnto2pip2pim, // n + nbar --> 2\pi^{+} + 2\pi^{-} - kNOnto2pip2pim1pi0, // n + nbar --> 2\pi^{+} + 2\pi^{-} + \pi^{0} - kNOnto1pip1pim1o, // n + nbar --> \pi^{+} + \pi^{-} + \omega^{0} - kNOnto2pip2pim2pi0 // n + nbar --> 2\pi^{+} + 2\pi^{-} + 2\pi^{0} - - } NNBarOscMode_t; + typedef enum ENNBarOscMode { + + // i just replaced all the nucleon decay modes with nnbar modes -j + kNONull = -1, + kNORandom, // Will select a random decay mode -j + kNOpto1pip1pi0, // p + nbar --> \pi^{+} + \pi^{0} + kNOpto1pip2pi0, // p + nbar --> \pi^{+} + 2\pi^{0} + kNOpto1pip3pi0, // p + nbar --> \pi^{+} + 3\pi^{0} + kNOpto1pip4pi0, // p + nbar --> \pi^{+} + 4\pi^{0} + kNOpto2pip1pim, // p + nbar --> 2\pi^{+} + \pi^{-} + kNOpto2pip1pim1pi0, // p + nbar --> 2\pi^{+} + \pi^{-} + \pi^{0} + kNOpto2pip1pim2pi0, // p + nbar --> 2\pi^{+} + \pi^{-} + 2\pi^{0} + kNOpto2pip1pim3pi0, // p + nbar --> 2\pi^{+} + \pi^{-} + 3\pi^{0} + kNOpto3pip2pim, // p + nbar --> 3\pi^{+} + 2\pi^{-} + kNOpto3pip2pim1pi0, // p + nbar --> 3\pi^{+} + 2\pi^{-} + \pi^{0} + kNOpto1pip1pi01o, // p + nbar --> 1\pi^{+} + \pi^{0} + 1\omega^{0} + kNOpto2pip1pim1o, // p + nbar --> 2\pi^{+} + \pi^{-} + 1\omega^{0} + + kNOnto2pi0, // n + nbar --> 2\pi^{0} + kNOnto3pi0, // n + nbar --> 3\pi^{0} + kNOnto4pi0, // n + nbar --> 4\pi^{0} + kNOnto5pi0, // n + nbar --> 5\pi^{0} + kNOnto7pi0, // n + nbar --> 7\pi^{0} + kNOnto1pip1pim, // n + nbar --> \pi^{+} + \pi^{-} + kNOnto1pip1pim1pi0, // n + nbar --> \pi^{+} + \pi^{-} + \pi^{0} + kNOnto1pip1pim2pi0, // n + nbar --> \pi^{+} + \pi^{-} + 2\pi^{0} + kNOnto1pip1pim3pi0, // n + nbar --> \pi^{+} + \pi^{-} + 3\pi^{0} + kNOnto1pip1pim4pi0, // n + nbar --> \pi^{+} + \pi^{-} + 4\pi^{0} + kNOnto1pip1pim5pi0, // n + nbar --> \pi^{+} + \pi^{-} + 5\pi^{0} + kNOnto2pip2pim, // n + nbar --> 2\pi^{+} + 2\pi^{-} + kNOnto2pip2pim1pi0, // n + nbar --> 2\pi^{+} + 2\pi^{-} + \pi^{0} + kNOnto2pip2pim2pi0, // n + nbar --> 2\pi^{+} + 2\pi^{-} + 2\pi^{0} + kNOnto2pip2pim3pi0, // n + nbar --> 2\pi^{+} + 2\pi^{-} + 3\pi^{0} + kNOnto3pip3pim, // n + nbar --> 3\pi^{+} + 3\pi^{-} + kNOnto3pip3pim1pi0, // n + nbar --> 3\pi^{+} + 3\pi^{-} + 1\pi^{0} + kNOnto1rho1pi0, // n + nbar --> \rho^{0} + \pi^{0} + kNOnto1rhopm1pimp, // n + nbar --> \rho^{+-} + \pi^{-+} + kNOnto2o, // n + nbar --> 2\omega^{0} + kNOnto1rho1o, // n + nbar --> \rho^{0} + \omega^{0} + kNOnto2pi01o, // n + nbar --> 2\pi^{0} + \omega^{0} + kNOnto1pip1pim1o, // n + nbar --> \pi^{+} + \pi^{-} + \omega^{0} + kNOnto1eta1o, // n + nbar --> \eta^{0} + \omega^{0} + kNOnto1pip1pim1eta, // n + nbar --> \pi^{+} + \pi^{-} + \eta^{0} + kNOnto1Kp1Km, // n + nbar --> \K^{+} + \K^{-} + kNOnto1Kp1Km1o, // n + nbar --> K^{+} + K^{-} + \omega^{0} + kNOnto2Ks1o, // n + nbar --> 2K^{s} + \omega^{0} + kNOnto1pipm1Kmp1Kl, // n + nbar --> \pi^{+-} + K^{-+} + K^{l} + kNOnto1pipm1Kmp1Ks, // n + nbar --> \pi^{+-} + K^{-+} + K^{s} + kNOnto1pipm1Kmp1pi0K0,// n + nbar --> \pi^{+-} + K^{-+} + \pi^{0} + K^{0} + kNOnto1pip1pim1Ks1Kl // n + nbar --> \pi^{+} + \pi^{-} + K^{s} + K^{l} + + } NNBarOscMode_t; } #endif diff --git a/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.cxx b/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.cxx index e2cd34bd5..4d141afaa 100644 --- a/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.cxx +++ b/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.cxx @@ -1,11 +1,14 @@ //____________________________________________________________________________ /* - Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org - Jeremy Hewes, Georgia Karagiorgi - University of Manchester -*/ + Jeremy Hewes, Georgia Karagiorgi + University of Manchester + + Linyan Wan + Fermilab + */ //____________________________________________________________________________ #include "Framework/Algorithm/AlgFactory.h" @@ -33,14 +36,14 @@ using namespace genie; //____________________________________________________________________________ NNBarOscPrimaryVtxGenerator::NNBarOscPrimaryVtxGenerator() : -EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator") + EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator") { } //____________________________________________________________________________ NNBarOscPrimaryVtxGenerator::NNBarOscPrimaryVtxGenerator( - string config) : -EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator",config) + string config) : + EventRecordVisitorI("genie::NNBarOscPrimaryVtxGenerator",config) { } @@ -59,8 +62,8 @@ void NNBarOscPrimaryVtxGenerator::ProcessEventRecord(GHepRecord * event) const // spit out that info -j LOG("NNBarOsc", pNOTICE) - << "Simulating decay " << genie::utils::nnbar_osc::AsString(fCurrDecayMode) - << " for an initial state with code: " << fCurrInitStatePdg; + << "Simulating decay " << genie::utils::nnbar_osc::AsString(fCurrDecayMode) + << " for an initial state with code: " << fCurrInitStatePdg; // check if nucleon is bound fNucleonIsBound = (pdg::IonPdgCodeToA(fCurrInitStatePdg) > 1); @@ -76,15 +79,15 @@ void NNBarOscPrimaryVtxGenerator::ProcessEventRecord(GHepRecord * event) const void NNBarOscPrimaryVtxGenerator::AddInitialState(GHepRecord * event) const { -// add initial state to event record + // add initial state to event record -// event record looks like this: -// id: 0, nucleus (kIStInitialState) -// | -// |---> { id: 1, neutron (kIStDecayedState) -// { id: 2, nucleon (kIStDecayedState) -// { id: 3, remnant nucleus (kIStStableFinalState) -// + // event record looks like this: + // id: 0, nucleus (kIStInitialState) + // | + // |---> { id: 1, neutron (kIStDecayedState) + // { id: 2, nucleon (kIStDecayedState) + // { id: 3, remnant nucleus (kIStStableFinalState) + // TLorentzVector v4(0,0,0,0); @@ -123,12 +126,12 @@ void NNBarOscPrimaryVtxGenerator::AddInitialState(GHepRecord * event) const } //____________________________________________________________________________ void NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition( - GHepRecord * event) const + GHepRecord * event) const { GHepParticle * initial_nucleus = event->Particle(0); int A = initial_nucleus->A(); if(A <= 2) { - return; + return; } RandomGen * rnd = RandomGen::Instance(); @@ -138,14 +141,14 @@ void NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition( double R = R0 * TMath::Power(dA, 1./3.); LOG("NNBarOsc", pINFO) - << "Generating vertex according to a realistic nuclear density profile"; + << "Generating vertex according to a realistic nuclear density profile"; // get inputs to the rejection method double ymax = -1; double rmax = 3*R; double dr = R/40.; for(double r = 0; r < rmax; r+=dr) { - ymax = TMath::Max(ymax, r*r * utils::nuclear::Density(r,A)); + ymax = TMath::Max(ymax, r*r * utils::nuclear::Density(r,A)); } ymax *= 1.2; @@ -153,38 +156,38 @@ void NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition( TLorentzVector vtx(0,0,0,0); unsigned int iter = 0; while(1) { - iter++; - - // throw an exception if it hasn't find a solution after many attempts - if(iter > controls::kRjMaxIterations) { - LOG("NNBarOsc", pWARN) - << "Couldn't generate a vertex position after " << iter << " iterations"; - genie::exceptions::EVGThreadException exception; - exception.SetReason("Couldn't generate vertex"); - exception.SwitchOnFastForward(); - throw exception; - } - - double r = rmax * rnd->RndFsi().Rndm(); - double t = ymax * rnd->RndFsi().Rndm(); - double y = r*r * utils::nuclear::Density(r,A); - if(y > ymax) { - LOG("NNBarOsc", pERROR) - << "y = " << y << " > ymax = " << ymax << " for r = " << r << ", A = " << A; - } - bool accept = (t < y); - if(accept) { - double phi = 2*constants::kPi * rnd->RndFsi().Rndm(); - double cosphi = TMath::Cos(phi); - double sinphi = TMath::Sin(phi); - double costheta = -1 + 2 * rnd->RndFsi().Rndm(); - double sintheta = TMath::Sqrt(1-costheta*costheta); - vtx.SetX(r*sintheta*cosphi); - vtx.SetY(r*sintheta*sinphi); - vtx.SetZ(r*costheta); - vtx.SetT(0.); - break; - } + iter++; + + // throw an exception if it hasn't find a solution after many attempts + if(iter > controls::kRjMaxIterations) { + LOG("NNBarOsc", pWARN) + << "Couldn't generate a vertex position after " << iter << " iterations"; + genie::exceptions::EVGThreadException exception; + exception.SetReason("Couldn't generate vertex"); + exception.SwitchOnFastForward(); + throw exception; + } + + double r = rmax * rnd->RndFsi().Rndm(); + double t = ymax * rnd->RndFsi().Rndm(); + double y = r*r * utils::nuclear::Density(r,A); + if(y > ymax) { + LOG("NNBarOsc", pERROR) + << "y = " << y << " > ymax = " << ymax << " for r = " << r << ", A = " << A; + } + bool accept = (t < y); + if(accept) { + double phi = 2*constants::kPi * rnd->RndFsi().Rndm(); + double cosphi = TMath::Cos(phi); + double sinphi = TMath::Sin(phi); + double costheta = -1 + 2 * rnd->RndFsi().Rndm(); + double sintheta = TMath::Sqrt(1-costheta*costheta); + vtx.SetX(r*sintheta*cosphi); + vtx.SetY(r*sintheta*sinphi); + vtx.SetZ(r*costheta); + vtx.SetT(0.); + break; + } } // while 1 // giving position to oscillating neutron @@ -199,12 +202,12 @@ void NNBarOscPrimaryVtxGenerator::GenerateOscillatingNeutronPosition( } //____________________________________________________________________________ void NNBarOscPrimaryVtxGenerator::GenerateFermiMomentum( - GHepRecord * event) const + GHepRecord * event) const { GHepParticle * initial_nucleus = event->Particle(0); int A = initial_nucleus->A(); if(A <= 2) { - return; + return; } GHepParticle * oscillating_neutron = event->Particle(1); @@ -219,8 +222,12 @@ void NNBarOscPrimaryVtxGenerator::GenerateFermiMomentum( // start with oscillating neutron tgt.SetHitNucPdg(kPdgNeutron); + //Changes Made Here + double annihilationPos = oscillating_neutron->X4()->Vect().Mag(); + tgt.SetHitNucPosition( annihilationPos ); + // generate nuclear model & fermi momentum - fNuclModel->GenerateNucleon(tgt); + fNuclModel->GenerateNucleon(tgt,annihilationPos); TVector3 p3 = fNuclModel->Momentum3(); double w = fNuclModel->RemovalEnergy(); @@ -232,14 +239,14 @@ void NNBarOscPrimaryVtxGenerator::GenerateFermiMomentum( oscillating_neutron->SetMomentum(p4); LOG("NNBarOsc", pINFO) - << "Generated neutron momentum: (" - << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), " - << "|p| = " << p3.Mag(); + << "Generated neutron momentum: (" + << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), " + << "|p| = " << p3.Mag(); // then rinse repeat for the annihilation nucleon tgt.SetHitNucPdg(annihilation_nucleon->Pdg()); // use nuclear model to generate fermi momentum - fNuclModel->GenerateNucleon(tgt); + fNuclModel->GenerateNucleon(tgt,annihilationPos); p3 = fNuclModel->Momentum3(); w = fNuclModel->RemovalEnergy(); // use mass & momentum to figure out energy @@ -250,9 +257,9 @@ void NNBarOscPrimaryVtxGenerator::GenerateFermiMomentum( annihilation_nucleon->SetMomentum(p4); LOG("NNBarOsc", pINFO) - << "Generated nucleon momentum: (" - << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), " - << "|p| = " << p3.Mag(); + << "Generated nucleon momentum: (" + << p3.Px() << ", " << p3.Py() << ", " << p3.Pz() << "), " + << "|p| = " << p3.Mag(); // now figure out momentum for the nuclear remnant p3 = -1 * (oscillating_neutron->GetP4()->Vect() + annihilation_nucleon->GetP4()->Vect()); @@ -265,7 +272,7 @@ void NNBarOscPrimaryVtxGenerator::GenerateFermiMomentum( } //____________________________________________________________________________ void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts( - GHepRecord * event) const + GHepRecord * event) const { LOG("NNBarOsc", pINFO) << "Generating decay..."; @@ -282,14 +289,14 @@ void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts( double * mass = new double[pdgv.size()]; double sum = 0; for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) { - int pdgc = *pdg_iter; - double m = PDGLibrary::Instance()->Find(pdgc)->Mass(); - mass[idx++] = m; - sum += m; + int pdgc = *pdg_iter; + double m = PDGLibrary::Instance()->Find(pdgc)->Mass(); + mass[idx++] = m; + sum += m; } LOG("NNBarOsc", pINFO) - << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum; + << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum; int initial_nucleus_id = 0; int oscillating_neutron_id = 1; int annihilation_nucleon_id = 2; @@ -319,7 +326,7 @@ void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts( delete p4_2; LOG("NNBarOsc", pINFO) - << "Decaying system p4 = " << utils::print::P4AsString(p4d); + << "Decaying system p4 = " << utils::print::P4AsString(p4d); // Set the decay bool permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass); @@ -327,120 +334,119 @@ void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts( // If the decay is not energetically allowed, select a new final state while(!permitted) { - LOG("NNBarOsc", pINFO) - << "Not enough energy to generate decay products! Selecting a new final state."; - - int mode = 0; - - int initial_nucleus_pdg = initial_nucleus->Pdg(); - - std::string nucleus_pdg = std::to_string(static_cast(initial_nucleus_pdg)); - if (nucleus_pdg.size() != 10) { - LOG("NNBarOsc", pERROR) - << "Expecting the nuclear PDG code to be a 10-digit integer, but it is " << nucleus_pdg << ", which is " - << nucleus_pdg.size() << " digits long. Drop me an email at jezhewes@gmail.com ; exiting..."; - exit(1); - } - - int n_nucleons = std::stoi(nucleus_pdg.substr(6,3)) - 1; - int n_protons = std::stoi(nucleus_pdg.substr(3,3)); - double proton_frac = ((double)n_protons) / ((double) n_nucleons); - double neutron_frac = 1 - proton_frac; - - // set branching ratios, taken from bubble chamber data - const int n_modes = 16; - double br [n_modes] = { 0.010, 0.080, 0.100, 0.220, - 0.360, 0.160, 0.070, 0.020, - 0.015, 0.065, 0.110, 0.280, - 0.070, 0.240, 0.100, 0.100 }; - - for (int i = 0; i < n_modes; i++) { - // for first 7 branching ratios, multiply by relative proton density - if (i < 7) - br[i] *= proton_frac; - // for next 9, multiply by relative neutron density - else - br[i] *= neutron_frac; - } - - // randomly generate a number between 1 and 0 - RandomGen * rnd = RandomGen::Instance(); - rnd->SetSeed(0); - double p = rnd->RndNum().Rndm(); - - // loop through all modes, figure out which one our random number corresponds to - double threshold = 0; - for (int j = 0; j < n_modes; j++) { - threshold += br[j]; - if (p < threshold) { - // once we've found our mode, stop looping - mode = j + 1; - break; - } - } - - // create new event record with new final state - // TODO - I don't think Jeremy meant to make a _new_ record here; it - // shadows the one passed in... - // EventRecord * event = new EventRecord; - Interaction * interaction = Interaction::NOsc((int)fCurrInitStatePdg, mode); - event->AttachSummary(interaction); - - fCurrDecayMode = (NNBarOscMode_t) interaction->ExclTag().DecayMode(); - - pdgv = genie::utils::nnbar_osc::DecayProductList(fCurrDecayMode); - LOG("NNBarOsc", pINFO) << "Decay product IDs: " << pdgv; - assert ( pdgv.size() > 1); - - // get the decay particles again - LOG("NNBarOsc", pINFO) << "Performing a phase space decay..."; - idx = 0; - delete [] mass; - mass = new double[pdgv.size()]; - sum = 0; - for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) { - int pdgc = *pdg_iter; - double m = PDGLibrary::Instance()->Find(pdgc)->Mass(); - mass[idx++] = m; - sum += m; - } - - LOG("NNBarOsc", pINFO) - << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum; - LOG("NNBarOsc", pINFO) - << "Decaying system p4 = " << utils::print::P4AsString(p4d); - - permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass); + LOG("NNBarOsc", pINFO) + << "Not enough energy to generate decay products! Selecting a new final state."; + + int mode = 0; + + int initial_nucleus_pdg = initial_nucleus->Pdg(); + + std::string nucleus_pdg = std::to_string(static_cast(initial_nucleus_pdg)); + if (nucleus_pdg.size() != 10) { + LOG("NNBarOsc", pERROR) + << "Expecting the nuclear PDG code to be a 10-digit integer, but it is " << nucleus_pdg << ", which is " + << nucleus_pdg.size() << " digits long. Drop me an email at jezhewes@gmail.com ; exiting..."; + exit(1); + } + + int n_nucleons = std::stoi(nucleus_pdg.substr(6,3)) - 1; + int n_protons = std::stoi(nucleus_pdg.substr(3,3)); + double proton_frac = ((double)n_protons) / ((double) n_nucleons); + double neutron_frac = 1 - proton_frac; + + // set branching ratios, taken from SK2021 + const int np_modes = fNPModes; + const int nn_modes = fNNModes; + const int n_modes = np_modes + nn_modes; + auto br = this->GetBRs(); + + for (int i = 0; i < n_modes; i++) { + // for first 7 branching ratios, multiply by relative proton density + if (i < np_modes) + br[i] *= proton_frac; + // for next 9, multiply by relative neutron density + else + br[i] *= neutron_frac; + } + + // randomly generate a number between 1 and 0 + RandomGen * rnd = RandomGen::Instance(); + rnd->SetSeed(0); + double p = rnd->RndNum().Rndm(); + + // loop through all modes, figure out which one our random number corresponds to + double threshold = 0; + for (int j = 0; j < n_modes; j++) { + threshold += br[j]; + if (p < threshold) { + // once we've found our mode, stop looping + mode = j + 1; + break; + } + } + + // create new event record with new final state + // TODO - I don't think Jeremy meant to make a _new_ record here; it + // shadows the one passed in... + // EventRecord * event = new EventRecord; + Interaction * interaction = Interaction::NOsc((int)fCurrInitStatePdg, mode); + event->AttachSummary(interaction); + + fCurrDecayMode = (NNBarOscMode_t) interaction->ExclTag().DecayMode(); + + pdgv = genie::utils::nnbar_osc::DecayProductList(fCurrDecayMode); + LOG("NNBarOsc", pINFO) << "Decay product IDs: " << pdgv; + assert ( pdgv.size() > 1); + + // get the decay particles again + LOG("NNBarOsc", pINFO) << "Performing a phase space decay..."; + idx = 0; + delete [] mass; + mass = new double[pdgv.size()]; + sum = 0; + for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) { + int pdgc = *pdg_iter; + double m = PDGLibrary::Instance()->Find(pdgc)->Mass(); + mass[idx++] = m; + sum += m; + } + + LOG("NNBarOsc", pINFO) + << "Decaying N = " << pdgv.size() << " particles / total mass = " << sum; + LOG("NNBarOsc", pINFO) + << "Decaying system p4 = " << utils::print::P4AsString(p4d); + + permitted = fPhaseSpaceGenerator.SetDecay(*p4d, pdgv.size(), mass); } if(!permitted) { - LOG("NNBarOsc", pERROR) - << " *** Phase space decay is not permitted \n" - << " Total particle mass = " << sum << "\n" - << " Decaying system p4 = " << utils::print::P4AsString(p4d); - // clean-up - delete [] mass; - delete p4d; - delete v4d; - // throw exception - genie::exceptions::EVGThreadException exception; - exception.SetReason("Decay not permitted kinematically"); - exception.SwitchOnFastForward(); - throw exception; + LOG("NNBarOsc", pERROR) + << " *** Phase space decay is not permitted \n" + << " Total particle mass = " << sum << "\n" + << " Decaying system p4 = " << utils::print::P4AsString(p4d); + // clean-up + delete [] mass; + delete p4d; + delete v4d; + // throw exception + genie::exceptions::EVGThreadException exception; + exception.SetReason("Decay not permitted kinematically"); + exception.SwitchOnFastForward(); + throw exception; } // Get the maximum weight //double wmax = fPhaseSpaceGenerator.GetWtMax(); double wmax = -1; for(int i=0; i<200; i++) { - double w = fPhaseSpaceGenerator.Generate(); - wmax = TMath::Max(wmax,w); + double w = fPhaseSpaceGenerator.Generate(); + wmax = TMath::Max(wmax,w); } assert(wmax>0); wmax *= 2; LOG("NNBarOsc", pNOTICE) - << "Max phase space gen. weight @ current hadronic system: " << wmax; + << "Max phase space gen. weight @ current hadronic system: " << wmax; // Generate an unweighted decay RandomGen * rnd = RandomGen::Instance(); @@ -449,34 +455,34 @@ void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts( unsigned int itry=0; while(!accept_decay) { - itry++; - - if(itry > controls::kMaxUnweightDecayIterations) { - // report, clean-up and return - LOG("NNBarOsc", pWARN) - << "Couldn't generate an unweighted phase space decay after " - << itry << " attempts"; - // clean up - delete [] mass; - delete p4d; - delete v4d; - // throw exception - genie::exceptions::EVGThreadException exception; - exception.SetReason("Couldn't select decay after N attempts"); - exception.SwitchOnFastForward(); - throw exception; - } - double w = fPhaseSpaceGenerator.Generate(); - if(w > wmax) { - LOG("NNBarOsc", pWARN) - << "Decay weight = " << w << " > max decay weight = " << wmax; - } - double gw = wmax * rnd->RndHadro().Rndm(); - accept_decay = (gw<=w); - - LOG("NNBarOsc", pINFO) - << "Decay weight = " << w << " / R = " << gw - << " - accepted: " << accept_decay; + itry++; + + if(itry > controls::kMaxUnweightDecayIterations) { + // report, clean-up and return + LOG("NNBarOsc", pWARN) + << "Couldn't generate an unweighted phase space decay after " + << itry << " attempts"; + // clean up + delete [] mass; + delete p4d; + delete v4d; + // throw exception + genie::exceptions::EVGThreadException exception; + exception.SetReason("Couldn't select decay after N attempts"); + exception.SwitchOnFastForward(); + throw exception; + } + double w = fPhaseSpaceGenerator.Generate(); + if(w > wmax) { + LOG("NNBarOsc", pWARN) + << "Decay weight = " << w << " > max decay weight = " << wmax; + } + double gw = wmax * rnd->RndHadro().Rndm(); + accept_decay = (gw<=w); + + LOG("NNBarOsc", pINFO) + << "Decay weight = " << w << " / R = " << gw + << " - accepted: " << accept_decay; } //!accept_decay @@ -484,13 +490,13 @@ void NNBarOscPrimaryVtxGenerator::GenerateDecayProducts( TLorentzVector v4(*v4d); int idp = 0; for(pdg_iter = pdgv.begin(); pdg_iter != pdgv.end(); ++pdg_iter) { - int pdgc = *pdg_iter; - TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp); - GHepStatus_t ist = - utils::nnbar_osc::DecayProductStatus(fNucleonIsBound, pdgc); - p4fin->Boost(boost); - event->AddParticle(pdgc, ist, oscillating_neutron_id,-1,-1,-1, *p4fin, v4); - idp++; + int pdgc = *pdg_iter; + TLorentzVector * p4fin = fPhaseSpaceGenerator.GetDecay(idp); + GHepStatus_t ist = + utils::nnbar_osc::DecayProductStatus(fNucleonIsBound, pdgc); + p4fin->Boost(boost); + event->AddParticle(pdgc, ist, oscillating_neutron_id,-1,-1,-1, *p4fin, v4); + idp++; } // Clean-up @@ -513,8 +519,8 @@ void NNBarOscPrimaryVtxGenerator::Configure(string config) //___________________________________________________________________________ void NNBarOscPrimaryVtxGenerator::LoadConfig(void) { -// AlgConfigPool * confp = AlgConfigPool::Instance(); -// const Registry * gc = confp->GlobalParameterList(); + // AlgConfigPool * confp = AlgConfigPool::Instance(); + // const Registry * gc = confp->GlobalParameterList(); fNuclModel = 0; diff --git a/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.h b/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.h index 8e9dd8a4d..995cb1ff3 100644 --- a/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.h +++ b/src/Physics/NNBarOscillation/NNBarOscPrimaryVtxGenerator.h @@ -13,13 +13,15 @@ \created November, 2016 \cpright Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + For the full text of the license visit http://copyright.genie-mc.org */ //____________________________________________________________________________ #ifndef _NNBAR_OSC_PRIMARY_VTX_GENERATOR_H_ #define _NNBAR_OSC_PRIMARY_VTX_GENERATOR_H_ +#include + #include #include #include @@ -61,6 +63,26 @@ class NNBarOscPrimaryVtxGenerator: public EventRecordVisitorI { mutable TGenPhaseSpace fPhaseSpaceGenerator; const NuclearModelI * fNuclModel; + + static constexpr int fNPModes = 12; + static constexpr int fNNModes = 32; + static constexpr std::array< double, fNPModes + fNNModes > fBR = { + 0.001, 0.007, 0.148, 0.014, 0.020, 0.170, 0.108, 0.301, + 0.055, 0.032, 0.020, 0.124, //pnbar + 0.001, 0.007, 0.003, 0.010, 0.001, 0.003, 0.016, 0.131, 0.112, 0.033, + 0.014, 0.060, 0.136, 0.157, 0.006, 0.022, 0.020, 0.018, 0.037, 0.035, + 0.024, 0.027, 0.071, 0.016, 0.017, 0.001, 0.002, 0.001, 0.003, 0.003, + 0.010, 0.003 + }; + +public: + + inline int GetNPModes() const { return fNPModes; } + inline int GetNNModes() const { return fNNModes; } + inline int GetNModes() const { return fNPModes + fNNModes; } + inline std::array< double, fNPModes + fNNModes > GetBRs() const + { return fBR; } + }; } // genie namespace diff --git a/src/Physics/NNBarOscillation/NNBarOscUtils.cxx b/src/Physics/NNBarOscillation/NNBarOscUtils.cxx index a90848346..5de48ce5d 100644 --- a/src/Physics/NNBarOscillation/NNBarOscUtils.cxx +++ b/src/Physics/NNBarOscillation/NNBarOscUtils.cxx @@ -1,11 +1,14 @@ //____________________________________________________________________________ /* - Copyright (c) 2003-2023, The GENIE Collaboration - For the full text of the license visit http://copyright.genie-mc.org + Copyright (c) 2003-2023, The GENIE Collaboration + For the full text of the license visit http://copyright.genie-mc.org - Jeremy Hewes, Georgia Karagiorgi - University of Manchester -*/ + Jeremy Hewes, Georgia Karagiorgi + University of Manchester + + Linyan Wan + Fermilab + */ //____________________________________________________________________________ #include "Framework/Messenger/Messenger.h" @@ -24,42 +27,100 @@ string genie::utils::nnbar_osc::AsString(NNBarOscMode_t ndm) { // this just maps the decay mode integers to string descriptors. replaced. -j switch(ndm) { - case (kNORandom) : return "Random mode"; - break; - case (kNOpto1pip1pi0) : return "p + nbar --> pi+ pi0"; - break; - case (kNOpto1pip2pi0) : return "p + nbar --> pi+ 2pi0"; - break; - case (kNOpto1pip3pi0) : return "p + nbar --> pi+ 3pi0"; - break; - case (kNOpto2pip1pim1pi0) : return "p + nbar --> 2pi+ pi- pi0"; - break; - case (kNOpto2pip1pim2pi0) : return "p + nbar --> 2pi+ pi- 2pi0"; - break; - case (kNOpto2pip1pim2o) : return "p + nbar --> 2pi+ pi- 2omega"; - break; - case (kNOpto3pip2pim1pi0) : return "p + nbar --> 3pi+ 2pi- pi0"; - break; - case (kNOnto1pip1pim) : return "n + nbar --> pi+ pi-"; - break; - case (kNOnto2pi0) : return "n + nbar --> 2pi0"; - break; - case (kNOnto1pip1pim1pi0) : return "n + nbar --> pi+ pi- pi0"; - break; - case (kNOnto1pip1pim2pi0) : return "n + nbar --> pi+ pi- 2pi0"; - break; - case (kNOnto1pip1pim3pi0) : return "n + nbar --> pi+ pi- 3pi0"; - break; - case (kNOnto2pip2pim) : return "n + nbar --> 2pi+ 2pi-"; - break; - case (kNOnto2pip2pim1pi0) : return "n + nbar --> 2pi+ 2pi- pi0"; - break; - case (kNOnto1pip1pim1o) : return "n + nbar --> pi+ pi- omega"; - break; - case (kNOnto2pip2pim2pi0) : return "n + nbar --> 2pi+ 2pi- 2pi0"; - break; - default : return "?"; - break; + case (kNORandom) : return "Random mode"; + break; + case (kNOpto1pip1pi0) : return "p + nbar --> pi+ pi0"; + break; + case (kNOpto1pip2pi0) : return "p + nbar --> pi+ 2pi0"; + break; + case (kNOpto1pip3pi0) : return "p + nbar --> pi+ 3pi0"; + break; + case (kNOpto1pip4pi0) : return "p + nbar --> pi+ 4pi0"; + break; + case (kNOpto2pip1pim) : return "p + nbar --> 2pi+ pi-"; + break; + case (kNOpto2pip1pim1pi0) : return "p + nbar --> 2pi+ pi- pi0"; + break; + case (kNOpto2pip1pim2pi0) : return "p + nbar --> 2pi+ pi- 2pi0"; + break; + case (kNOpto2pip1pim3pi0) : return "p + nbar --> 2pi+ pi- 3pi0"; + break; + case (kNOpto3pip2pim) : return "p + nbar --> 3pi+ 2pi"; + break; + case (kNOpto3pip2pim1pi0) : return "p + nbar --> 3pi+ 2pi- pi0"; + break; + case (kNOpto1pip1pi01o) : return "p + nbar --> pi+ pi0 omega"; + break; + case (kNOpto2pip1pim1o) : return "p + nbar --> 2pi+ pi- omega"; + break; + + case (kNOnto2pi0) : return "n + nbar --> 2pi0"; + break; + case (kNOnto3pi0) : return "n + nbar --> 3pi0"; + break; + case (kNOnto4pi0) : return "n + nbar --> 4pi0"; + break; + case (kNOnto5pi0) : return "n + nbar --> 5pi0"; + break; + case (kNOnto7pi0) : return "n + nbar --> 7pi0"; + break; + case (kNOnto1pip1pim) : return "n + nbar --> pi+ pi-"; + break; + case (kNOnto1pip1pim1pi0) : return "n + nbar --> pi+ pi- pi0"; + break; + case (kNOnto1pip1pim2pi0) : return "n + nbar --> pi+ pi- 2pi0"; + break; + case (kNOnto1pip1pim3pi0) : return "n + nbar --> pi+ pi- 3pi0"; + break; + case (kNOnto1pip1pim4pi0) : return "n + nbar --> pi+ pi- 4pi0"; + break; + case (kNOnto1pip1pim5pi0) : return "n + nbar --> pi+ pi- 5pi0"; + break; + case (kNOnto2pip2pim) : return "n + nbar --> 2pi+ 2pi-"; + break; + case (kNOnto2pip2pim1pi0) : return "n + nbar --> 2pi+ 2pi- pi0"; + break; + case (kNOnto2pip2pim2pi0) : return "n + nbar --> 2pi+ 2pi- 2pi0"; + break; + case (kNOnto2pip2pim3pi0) : return "n + nbar --> 2pi+ 2pi- 3pi0"; + break; + case (kNOnto3pip3pim) : return "n + nbar --> 3pi+ 3pi-"; + break; + case (kNOnto3pip3pim1pi0) : return "n + nbar --> 3pi+ 3pi- 1pi0"; + break; + case (kNOnto1rho1pi0) : return "n + nbar --> rho0 pi0"; + break; + case (kNOnto1rhopm1pimp) : return "n + nbar --> rho+- pi-+"; + break; + case (kNOnto2o) : return "n + nbar --> 2omega"; + break; + case (kNOnto1rho1o) : return "n + nbar --> rho0 omega0"; + break; + case (kNOnto2pi01o) : return "n + nbar --> 2pi0 omega0"; + break; + case (kNOnto1pip1pim1o) : return "n + nbar --> pi+ pi- omega0"; + break; + case (kNOnto1eta1o) : return "n + nbar --> eta0 omega0"; + break; + case (kNOnto1pip1pim1eta) : return "n + nbar --> pi+ pi- eta0"; + break; + case (kNOnto1Kp1Km) : return "n + nbar --> K+ K-"; + break; + case (kNOnto1Kp1Km1o) : return "n + nbar --> K+ K- omega0"; + break; + case (kNOnto2Ks1o) : return "n + nbar --> Ks Ks omega0"; + break; + case (kNOnto1pipm1Kmp1Kl) : return "n + nbar --> pi+- K-+ Kl"; + break; + case (kNOnto1pipm1Kmp1Ks) : return "n + nbar --> pi+- K-+ Ks"; + break; + case (kNOnto1pipm1Kmp1pi0K0): return "n + nbar --> pi+- K-+ pi0 K0"; + break; + case (kNOnto1pip1pim1Ks1Kl): return "n + nbar --> pi+ pi- Ks Kl"; + break; + + default : return "?"; + break; } return "??"; } @@ -68,28 +129,58 @@ bool genie::utils::nnbar_osc::IsValidMode(NNBarOscMode_t ndm) { // checks if a mode is valid. just straight replaced. -j switch(ndm) { - case (kNORandom) : - case (kNOpto1pip1pi0) : - case (kNOpto1pip2pi0) : - case (kNOpto1pip3pi0) : - case (kNOpto2pip1pim1pi0) : - case (kNOpto2pip1pim2pi0) : - case (kNOpto2pip1pim2o) : - case (kNOpto3pip2pim1pi0) : - case (kNOnto1pip1pim) : - case (kNOnto2pi0) : - case (kNOnto1pip1pim1pi0) : - case (kNOnto1pip1pim2pi0) : - case (kNOnto1pip1pim3pi0) : - case (kNOnto2pip2pim) : - case (kNOnto2pip2pim1pi0) : - case (kNOnto1pip1pim1o) : - case (kNOnto2pip2pim2pi0) : - return true; - break; - default : - return false; - break; + case (kNORandom) : + case (kNOpto1pip1pi0) : + case (kNOpto1pip2pi0) : + case (kNOpto1pip3pi0) : + case (kNOpto1pip4pi0) : + case (kNOpto2pip1pim) : + case (kNOpto2pip1pim1pi0) : + case (kNOpto2pip1pim2pi0) : + case (kNOpto2pip1pim3pi0) : + case (kNOpto3pip2pim) : + case (kNOpto3pip2pim1pi0) : + case (kNOpto1pip1pi01o) : + case (kNOpto2pip1pim1o) : + + case (kNOnto2pi0) : + case (kNOnto3pi0) : + case (kNOnto4pi0) : + case (kNOnto5pi0) : + case (kNOnto7pi0) : + case (kNOnto1pip1pim) : + case (kNOnto1pip1pim1pi0) : + case (kNOnto1pip1pim2pi0) : + case (kNOnto1pip1pim3pi0) : + case (kNOnto1pip1pim4pi0) : + case (kNOnto1pip1pim5pi0) : + case (kNOnto2pip2pim) : + case (kNOnto2pip2pim1pi0) : + case (kNOnto2pip2pim2pi0) : + case (kNOnto2pip2pim3pi0) : + case (kNOnto3pip3pim) : + case (kNOnto3pip3pim1pi0) : + case (kNOnto1rho1pi0) : + case (kNOnto1rhopm1pimp) : + case (kNOnto2o) : + case (kNOnto1rho1o) : + case (kNOnto2pi01o) : + case (kNOnto1pip1pim1o) : + case (kNOnto1eta1o) : + case (kNOnto1pip1pim1eta) : + case (kNOnto1Kp1Km) : + case (kNOnto1Kp1Km1o) : + case (kNOnto2Ks1o) : + case (kNOnto1pipm1Kmp1Kl) : + case (kNOnto1pipm1Kmp1Ks) : + case (kNOnto1pipm1Kmp1pi0K0): + case (kNOnto1pip1pim1Ks1Kl): + + return true; + break; + default : + return false; + break; } return false; } @@ -99,29 +190,58 @@ int genie::utils::nnbar_osc::AnnihilatingNucleonPdgCode(NNBarOscMode_t ndm) // name isn't really accurate any more. instead of decayed nucleon, function // returns what particle the oscillated neutron annihilated with -j switch(ndm) { - case (kNOpto1pip1pi0) : return kPdgProton; break; - case (kNOpto1pip2pi0) : return kPdgProton; break; - case (kNOpto1pip3pi0) : return kPdgProton; break; - case (kNOpto2pip1pim1pi0) : return kPdgProton; break; - case (kNOpto2pip1pim2pi0) : return kPdgProton; break; - case (kNOpto2pip1pim2o) : return kPdgProton; break; - case (kNOpto3pip2pim1pi0) : return kPdgProton; break; - case (kNOnto1pip1pim) : return kPdgNeutron; break; - case (kNOnto2pi0) : return kPdgNeutron; break; - case (kNOnto1pip1pim1pi0) : return kPdgNeutron; break; - case (kNOnto1pip1pim2pi0) : return kPdgNeutron; break; - case (kNOnto1pip1pim3pi0) : return kPdgNeutron; break; - case (kNOnto2pip2pim) : return kPdgNeutron; break; - case (kNOnto2pip2pim1pi0) : return kPdgNeutron; break; - case (kNOnto1pip1pim1o) : return kPdgNeutron; break; - case (kNOnto2pip2pim2pi0) : return kPdgNeutron; break; - default : return 0; break; + case (kNOpto1pip1pi0) : return kPdgProton; break; + case (kNOpto1pip2pi0) : return kPdgProton; break; + case (kNOpto1pip3pi0) : return kPdgProton; break; + case (kNOpto1pip4pi0) : return kPdgProton; break; + case (kNOpto2pip1pim) : return kPdgProton; break; + case (kNOpto2pip1pim1pi0) : return kPdgProton; break; + case (kNOpto2pip1pim2pi0) : return kPdgProton; break; + case (kNOpto2pip1pim3pi0) : return kPdgProton; break; + case (kNOpto3pip2pim) : return kPdgProton; break; + case (kNOpto3pip2pim1pi0) : return kPdgProton; break; + case (kNOpto1pip1pi01o) : return kPdgProton; break; + case (kNOpto2pip1pim1o) : return kPdgProton; break; + + case (kNOnto2pi0) : return kPdgNeutron; break; + case (kNOnto3pi0) : return kPdgNeutron; break; + case (kNOnto4pi0) : return kPdgNeutron; break; + case (kNOnto5pi0) : return kPdgNeutron; break; + case (kNOnto7pi0) : return kPdgNeutron; break; + case (kNOnto1pip1pim) : return kPdgNeutron; break; + case (kNOnto1pip1pim1pi0) : return kPdgNeutron; break; + case (kNOnto1pip1pim2pi0) : return kPdgNeutron; break; + case (kNOnto1pip1pim3pi0) : return kPdgNeutron; break; + case (kNOnto1pip1pim4pi0) : return kPdgNeutron; break; + case (kNOnto1pip1pim5pi0) : return kPdgNeutron; break; + case (kNOnto2pip2pim) : return kPdgNeutron; break; + case (kNOnto2pip2pim1pi0) : return kPdgNeutron; break; + case (kNOnto2pip2pim2pi0) : return kPdgNeutron; break; + case (kNOnto2pip2pim3pi0) : return kPdgNeutron; break; + case (kNOnto3pip3pim) : return kPdgNeutron; break; + case (kNOnto3pip3pim1pi0) : return kPdgNeutron; break; + case (kNOnto1rho1pi0) : return kPdgNeutron; break; + case (kNOnto1rhopm1pimp) : return kPdgNeutron; break; + case (kNOnto2o) : return kPdgNeutron; break; + case (kNOnto1rho1o) : return kPdgNeutron; break; + case (kNOnto2pi01o) : return kPdgNeutron; break; + case (kNOnto1pip1pim1o) : return kPdgNeutron; break; + case (kNOnto1eta1o) : return kPdgNeutron; break; + case (kNOnto1pip1pim1eta) : return kPdgNeutron; break; + case (kNOnto1Kp1Km) : return kPdgNeutron; break; + case (kNOnto1Kp1Km1o) : return kPdgNeutron; break; + case (kNOnto2Ks1o) : return kPdgNeutron; break; + case (kNOnto1pipm1Kmp1Kl) : return kPdgNeutron; break; + case (kNOnto1pipm1Kmp1Ks) : return kPdgNeutron; break; + case (kNOnto1pipm1Kmp1pi0K0):return kPdgNeutron; break; + case (kNOnto1pip1pim1Ks1Kl): return kPdgNeutron; break; + default : return 0; break; } return 0; } //____________________________________________________________________________ PDGCodeList genie::utils::nnbar_osc::DecayProductList( - NNBarOscMode_t ndm) + NNBarOscMode_t ndm) { // ok so i think this is the first function where a straight replacement // isn't gonna cut it. all the nucleon decay modes are two-body, but that is @@ -134,119 +254,295 @@ PDGCodeList genie::utils::nnbar_osc::DecayProductList( PDGCodeList decay_products(allow_duplicate); switch(ndm) { - case (kNOpto1pip1pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPi0); - break; - case (kNOpto1pip2pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - case (kNOpto1pip3pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - case (kNOpto2pip1pim1pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - break; - case (kNOpto2pip1pim2pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - case (kNOpto2pip1pim2o) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgomega); - decay_products.push_back(kPdgomega); - break; - case (kNOpto3pip2pim1pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - break; - case (kNOnto1pip1pim) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - break; - case (kNOnto2pi0) : - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - case (kNOnto1pip1pim1pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - break; - case (kNOnto1pip1pim2pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - case (kNOnto1pip1pim3pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - case (kNOnto2pip2pim) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPiM); - break; - case (kNOnto2pip2pim1pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - break; - case (kNOnto1pip1pim1o) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgomega); - break; - case (kNOnto2pip2pim2pi0) : - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiP); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPiM); - decay_products.push_back(kPdgPi0); - decay_products.push_back(kPdgPi0); - break; - default : - break; + case (kNOpto1pip1pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto1pip2pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto1pip3pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto1pip4pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto2pip1pim) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + break; + case (kNOpto2pip1pim1pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto2pip1pim2pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto2pip1pim3pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto3pip2pim) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + break; + case (kNOpto3pip2pim1pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + break; + case (kNOpto1pip1pi01o) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgomega); + break; + case (kNOpto2pip1pim1o) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgomega); + break; + + case (kNOnto2pi0) : + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto3pi0) : + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto4pi0) : + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto5pi0) : + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto7pi0) : + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1pip1pim) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + break; + case (kNOnto1pip1pim1pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1pip1pim2pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1pip1pim3pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1pip1pim4pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1pip1pim5pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto2pip2pim) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + break; + case (kNOnto2pip2pim1pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto2pip2pim2pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto2pip2pim3pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto3pip3pim) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + break; + case (kNOnto3pip3pim1pi0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1rho1pi0) : + decay_products.push_back(kPdgRho0); + decay_products.push_back(kPdgPi0); + break; + case (kNOnto1rhopm1pimp) : + decay_products.push_back(kPdgRhoP); + decay_products.push_back(kPdgPiM); + break; + case (kNOnto2o) : + decay_products.push_back(kPdgomega); + decay_products.push_back(kPdgomega); + break; + case (kNOnto1rho1o) : + decay_products.push_back(kPdgRho0); + decay_products.push_back(kPdgomega); + break; + case (kNOnto2pi01o) : + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgomega); + break; + case (kNOnto1pip1pim1o) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgomega); + break; + case (kNOnto1eta1o) : + decay_products.push_back(kPdgEta); + decay_products.push_back(kPdgomega); + break; + case (kNOnto1pip1pim1eta) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgEta); + break; + case (kNOnto1Kp1Km) : + decay_products.push_back(kPdgKP); + decay_products.push_back(kPdgKM); + break; + case (kNOnto1Kp1Km1o) : + decay_products.push_back(kPdgKP); + decay_products.push_back(kPdgKM); + decay_products.push_back(kPdgomega); + break; + case (kNOnto2Ks1o) : + decay_products.push_back(kPdgK0S); + decay_products.push_back(kPdgK0S); + decay_products.push_back(kPdgomega); + break; + case (kNOnto1pipm1Kmp1Kl) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgKM); + decay_products.push_back(kPdgK0L); + break; + case (kNOnto1pipm1Kmp1Ks) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgKM); + decay_products.push_back(kPdgK0S); + break; + case (kNOnto1pipm1Kmp1pi0K0) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgKM); + decay_products.push_back(kPdgPi0); + decay_products.push_back(kPdgK0); + break; + case (kNOnto1pip1pim1Ks1Kl) : + decay_products.push_back(kPdgPiP); + decay_products.push_back(kPdgPiM); + decay_products.push_back(kPdgK0S); + decay_products.push_back(kPdgK0L); + break; + default : + break; } return decay_products; } //____________________________________________________________________________ GHepStatus_t genie::utils::nnbar_osc::DecayProductStatus( - bool in_nucleus, int pdgc) + bool in_nucleus, int pdgc) { // took out all the irrelevant particles -j if(in_nucleus) { - if( pdgc == kPdgPi0 || - pdgc == kPdgPiM || - pdgc == kPdgPiP || - pdgc == kPdgomega) - { - return kIStHadronInTheNucleus; - } + if( pdgc == kPdgPi0 || + pdgc == kPdgPiM || + pdgc == kPdgPiP || + pdgc == kPdgRho0 || + pdgc == kPdgRhoP || + pdgc == kPdgEta || + pdgc == kPdgKP || + pdgc == kPdgKM || + pdgc == kPdgK0 || + pdgc == kPdgK0S || + pdgc == kPdgK0L || + pdgc == kPdgomega) + { + return kIStHadronInTheNucleus; + } } return kIStStableFinalState;