Skip to content

Latest commit

 

History

History
317 lines (260 loc) · 15.1 KB

dev_bergersehgal.adoc

File metadata and controls

317 lines (260 loc) · 15.1 KB

Implementing New Models

Implementing a new model in GENIE is relatively straightforward, factoring out the complexity of the actual physics and considering only the software architecture. Models in GENIE are broken down into:

  1. A series of event record visitors. Each visitor carries out once link in the computation chain required to connect the initial state (the projectile and the target) and the final state. In principle all of these steps may be combined into one class, but it is better to break them into modular (even reusable) steps.

  2. One of the links in the event record visitor chain should compute differential cross sections. It is additionally necessary to write a cross sections module that will integrate the differential cross section to produce a total cross section when generating splines. Before the differential cross section code is available, you can cheat and write the spline by hand!

  3. Modifications to the configuration XML that provide quantitative values for parameters and specify whether the model should be included in a simulation and how to chain the event record visitors together.

  4. Machinery to hook into ROOT, both for interactive use inside CINT and because ROOT is deeply integrated into GENIE (including, for example, class loading).

Because it is difficult to explain all these steps in an abstract way, we will present case studies and hope they are sufficient to get developers started.

Implementing the Berger-Sehgal Model

The first and best step is to examine the most similar existing models. (If you have nothing similar, then just read any model carefully.) Two very useful exercies are:

  1. Step through a calculation of the integrated cross section for spline generation. Go to [Generating-a-Spline] and create a spline file for only the process you wish to study. Step through the spline generation with gdb using gdb -tui gmkspl and set the arguments to whatever is required to generate your study channel only. This will give you a good idea as to how the numerical integrators are called. (You needn’t be concerned with the actual methods of integration though unless that is a focus of your work.)

  2. Step through the computation of the differential cross section. Stepping through the spline generation carefully will also bring you through the codes that do these steps but you may also wish to just run them on their own, again watching in gdb wtih gdb -tui gevgen and setting the arguments appropriately in gdb.

In this example, the Berger-Sehgal coherent pion model is very similar to the Rein- Sehgal model. Both models have identical initial and final states and both exploit the PCAC hypothesis.

If we check the structure of Rein-Sehgal, we find a directory under src (unfortunately misspelled) that contains the coherent pion code (and other things as well, since Rein and Sehgal produced several important models together).

By examing the configuation XML discussed [Configuration-XML], we realize that the Rein-Sehgal coherent pion model is called for both the COH-NC and COH-CC event lists by default in GENIE 2.8.0. By examining the list of event record visitors for those events, we can find the names of several other classes that look to play an important role in the calculation. Running find on the GENIE src directory, we can find the important classes quickly:

$ ls -l ReinSeghal/ReinSeghalCOHPiPXSec.*
-rw-r--r-- 1 perdue e938   8333 Aug 23  2013 ReinSeghal/ReinSeghalCOHPiPXSec.cxx
-rw-r--r-- 1 perdue e938   2370 Aug 23  2013 ReinSeghal/ReinSeghalCOHPiPXSec.h

$ ls -l Coherent/*.cxx | grep -v El
-rw-r--r-- 1 perdue e938  8547 Oct  8 13:53 Coherent/COHHadronicSystemGenerator.cxx
-rw-r--r-- 1 perdue e938  4082 Sep 30 10:25 Coherent/COHInteractionListGenerator.cxx
-rw-r--r-- 1 perdue e938 25105 Oct  8 13:53 Coherent/COHKinematicsGenerator.cxx
-rw-r--r-- 1 perdue e938  4697 Oct  8 13:53 Coherent/COHPrimaryLeptonGenerator.cxx

In this case it is simple to mimic the structure of Rein-Sehgal and add a directory to src to hold the new model:

$ ls -l BergerSehgal/*.cxx
-rw-r--r-- 1 perdue e938 8373 Sep 30 10:26 BergerSehgal/BergerSehgalCOHPiPXSec.cxx

During event generation we go through classes in the src/Coherent directory and find clauses like:

  if (fXSecModel->Id().Name() == "genie::ReinSeghalCOHPiPXSec") {
       CalculateHadronicSystem_ReinSeghal(evrec);
  } else if ((fXSecModel->Id().Name() == "genie::AlvarezRusoCOHXSec")) {
       CalculateHadronicSystem_AlvarezRuso(evrec);
  }

Therefore, we need to add the lines:

  } else if ((fXSecModel->Id().Name() == "genie::BergerSehgalCOHPiXSec")) {
       CalculateHadronicSystem_BergerSehgal(evrec);

along with the appropriate function implementations everywhere appropriate (there are a half-dozen or so required locations).

Additionally, while generating the cross section splines, we stepped through a class in src/CrossSections named COHXSec that needs to adopt this look-up mechanism so it can call the appropriate model when producing the contribution to the total cross section.

Once we have implemented the code, we need to update the configuration XML. We update config/master_config.xml by adding the line

<config alg="genie::BergerSehgalCOHPiPXSec"> BergerSehgalCOHPiPXSec.xml </config>

to the <!-- CONFIGURATION FOR XSEC ALGORITHMS -→ block. The BergerSehgalCOHPiPXSec.xml does not exist by default, we need to create it as our next step. In this case we can essentially copy and rename the appropriate Rein-Sehgal analog, but if one is attempting to do something very non-standard, this step may require more thought.

Next, we need to update config/UserPhysicsOptions.xml if we want to use this model. In this case, we just swap out Rein-Sehgal for Berger-Sehgal by commenting out the RS model and adding the lines:

<param type="alg" name="XSecModel@genie::EventGenerator/COH-CC">
  genie::BergerSehgalCOHPiPXSec/Default
</param>
<param type="alg" name="XSecModel@genie::EventGenerator/COH-NC">
  genie::BergerSehgalCOHPiPXSec/Default
</param>

In cases where we want mutliple models present, we should probably add a controller class that contains the logic for choosing.

We also need to be sure the LinkDef.h file in our new package has the right class names defined. In this case, we can just copy and edit the Rein-Sehgal file (trimming unneeded lines for the resonance model, etc.).

To get our new module to build, we add BergerSehgal to the top-level Makefile. We also need to update the libs variable in src/scripts/genie-config.

Finally, on should be sure before running the code that the model is represented in the cross section spline. One may need to edit a file by hand for quick-and-dirty testing, but we will need to eventually be able to run gmkspl to produce the splines correctly. This may require you to update and/or add an integrator, so it is reasonable to put that step off at the very beginning and use a ``by hand'' file with fake numbers and a monochromatic neutrino beam for early-stage testing.

Implementation Notes for the LANL STA Model

We need to decide where to put the physics. For now, we choose

$GENIE/src/Physics/Multinucleon

Generally, we begin by defining a new possible interaction list generator. When GENIE first runs, it populates a list of possible interaction types, and will draw from this list using a cross section weighted probability. The interactions generated by the list generator are defined by the Interaction class defined in $GENIE/src/Framework/Interaction.

In this case, we need to extend the EScatteringType enum (it is generally safe to extend enums by appending - do not prepend or insert or you lose the ability to read back older versions of the files with current software while preserving the semantic meaning of the enum). We will add an enum type to ScatteringType.h in the $GENIE/src/Framework/Interaction folder. In this case, we choose kScQE1B2BInterfere, although that choice may not be permanent.

Returning to the InteractionList generator, we make a new class, here in $GENIE/src/Physics/Multinucleon/EventGen (the choice of the Multinucleon folder is pseudo-arbitrary here - we should put physics models where they may be easily found, and we may indeed choose to create an entirely new folder, etc.). We call the new class STAInteractionListGenerator and create header and implementation files. One key requirement is the file will use the named constructor pattern to create new Interaction objects with the right data. We need to make these constructors in the Interaction class found in src/Framework/Interaction. We edit the header, for example, to include the lines:

  static Interaction * STACC     (int tgt, int probe, double E=0);
  static Interaction * STANC     (int tgt, int probe, double E=0);
  static Interaction * STAEM     (int tgt, int probe, double E=0);

as public methods. Next, we define those functions in the implementation file. Here we are assuming we will want to have different interaction configurations for the three different possible scattering types, but this is not necessarily the case - we may collapse to one constructor depending on the needs of the model. Right now, this configuration requires an STAInteractionListGenerator.xml file in the $GENIE/config directory similar to the MECInteractionListGenerator.xml - with three different algorithm configurations for CC, NC, and EM interactions and boolean switches for the different modes.

While we’re working on these items, we’ll add an IsSTA() method to $GENIE/src/Framework/Interaction/ProcessInfo. This abbreviation ('STA') may not stick - there is likely a better option.

In order to get our STAInteractionListGenerator class to build, we must add a line to $GENIE/src/Physics/Multinucleon/LinkDef.h. The GENIE build system uses this file to create ROOT dictionaries able to read the class. The line we need to add is:

  #pragma link C++ class genie::STAInteractionListGenerator;

Next, we need to consider what the physics process will look like at runtime. If we examine $GENIE/config/EventGeneratorListAssembler.xml, we find the definition of the --event-generator-list flags we may pass to apps like gevgen. For example, we see definitions like:

  <param_set name="CCQE+CCMEC">
     <param type="int" name="NGenerators">  2                                  </param>
     <param type="alg" name="Generator-0">  genie::EventGenerator/MEC-CC       </param>
     <param type="alg" name="Generator-1">  genie::EventGenerator/QEL-CC       </param>
  </param_set>

  <param_set name="EMMEC">
     <param type="int" name="NGenerators">   1                                  </param>
     <param type="alg" name="Generator-0">   genie::EventGenerator/MEC-EM       </param>
  </param_set>

We will need STA lists defined in this style. In the neutrino-scattering portion of the file, we add:

  <param_set name="CCSTA">
     <param type="int" name="NGenerators">  1                                  </param>
     <param type="alg" name="Generator-0">  genie::EventGenerator/STA-CC       </param>
  </param_set>

  <param_set name="NCSTA">
     <param type="int" name="NGenerators">  1                                  </param>
     <param type="alg" name="Generator-0">  genie::EventGenerator/STA-NC       </param>
  </param_set>

  <param_set name="WeakSTA">
     <param type="int" name="NGenerators">   2                                  </param>
     <param type="alg" name="Generator-0">   genie::EventGenerator/STA-CC       </param>
     <param type="alg" name="Generator-1">   genie::EventGenerator/STA-NC       </param>
  </param_set>

And, in the electron-scattering portion of the file, we add:

  <param_set name="EMSTA">
     <param type="int" name="NGenerators">   1                                  </param>
     <param type="alg" name="Generator-0">   genie::EventGenerator/STA-EM       </param>
  </param_set>

Now, we need to define the lists, STA-EM, STA-CC, and STA-NC in the file $GENIE/config/EventGenerator.xml:

  <param_set name="STA-CC">
     <param type="string" name="VldContext"> </param>
     <param type="int"    name="NModules">   4                                                   </param>
     <param type="alg"    name="Module-0">   genie::InitialStateAppender/Default                 </param>
     <param type="alg"    name="Module-2">   genie::STAGenerator/Default                         </param>
     <param type="alg"    name="Module-3">   genie::HadronTransporter/Default                    </param>
     <param type="alg"    name="Module-4">   genie::UnstableParticleDecayer/AfterHadronTransport </param>
     <param type="alg"    name="ILstGen">    genie::STAInteractionListGenerator/CC-Default       </param>
  </param_set>

we also add the NC and EM analogs of this list. Note that we have dropped VertexGenerator because the physics code we are using here will incorporate an implicit ground state that will be used to seed the interaction in {x,y,z} space.

A new class appeared in the event generator XML: STAGenerator. We will create a stub class based in $GENIE/src/Physics/Multinucleon/EventGen with that name, and declare the class in the local LinkDef.h. We also need to add an STAGenerator.xml file to the $GENIE/config folder. We follow the MECGenerator.xml example - note that we drop the Common entry here (referring to an entry in CommonParams.xml).

In STAGenerator::ProcessEventRecord(), we add a check that the cross section algorithm is PastoreCarlsonSTAPXSec2018. We must add this class to $GENIE/src/Physics/Multinucleon/XSection. We must also add an integrator class to the same folder, here STAXSec. Of course, these classes must be added to the local LinkDef.h.

We must also add a configuration XML, PastoreCarlsonSTAPXSec2018.xml to $GENIE/config, along with a STAXSec.xml.

Finally, we must add all the XML files we have created to the list of files in master_config.xml.