HOWTO use hdgeant4 as a beam Monte Carlo generator

From GlueXWiki
Revision as of 14:30, 25 May 2021 by Jonesrt (Talk | contribs) (→‎GENBEAM event generation rates)

Jump to: navigation, search

Internal beam generator in hdgeant

In HOWTO to generate electromagnetic background the use of the BEAM card in control.in to activate the built-in coherent photon beam generator was described. One use case for this internal beam generator is to model beam-related scattering background that overlays the signals from hadronic reactions produced by the beam interacting in the GlueX target. Now that we have large samples of random and PS triggers to provide a more accurate representation of these backgrounds, that functionality is no longer widely used for physics simulations. However, there are other use cases for the built-in coherent beam generator that remain relevant, including the following.

  1. simulation of the PS acceptance
  2. simulation of the beam transverse intensity and polarization profile
  3. simulation of the TPOL trigger, acceptance, and analyzing power
  4. simulation of Compton and Bethe Heitler reactions in the GlueX target
  5. any other process whose acceptance might be sensitive to the transverse beam position and momentum profile

In order to support the full range of possible applications listed above, I decided to introduce a new feature to control.in, in addition to the BEAM card. This new feature is enabled by the GENBEAM card. When the GENBEAM feature is active, hdgeant switches from functioning primarily as a particle track simulator that follows particles and induces interactions as they traverse the setup geometry to a Monte Carlo generator that generates output records consisting not of detector hits, but either a single beam particle or an interaction vertex with final state positions and momenta to be simulated later. In this way, GENBEAM transforms hdgeant from a simulation to an event generator, with the one proviso that the results of the generator are dependent in a significant way on the beamline geometry and the specific beam properties from the conditions database for the run number specified in the control.in file.

GENBEAM 'precol'

This is the simplest GENBEAM option. All photons above the minimum energy specified in the BEAM card are included in the output, with one photon per event. There are no detector hits in the output hddm file. The only content in the output events is the <reaction> tag, as illustrated in the example below. 

    <reaction type="0" weight="0">
      <beam type="Photon">
        <momentum E="8.8443" px="-2.36241e-05" py="-0.000160787" pz="8.8443" />
        <polarization Px="2.98534e-17" Py="0.487544" Pz="8.86341e-06" />
      </beam>
      <vertex>
        <product decayVertex="0" id="1" mech="0" parentid="0" pdgtype="11" type="Electron">
          <momentum E="12.0003" px="-0.000104861" py="2.44356e-05" pz="12.0003" />
        </product>
        <origin t="-78.2157" vx="0.0370969" vy="0.0463277" vz="-2400" />
      </vertex>
      <random seed1="691877422" seed2="281623319" seed3="709975947" seed4="912931183" />
    </reaction>

The meaning of this reaction tag is as follows. The leading <beam> tag tells the energy, momentum, and polarization of the beam photon that was generated for this event. The <vertex> tag following that gives the 4-momentum of the radiating electron that produced the photon, together with projection of the electron momentum vector forward onto the reference plane at z=-24m, just upstream of the entrance to the primary collimator. The impact position of the beam photon on the reference plane is not given explicitly, but it can be computed by linearly propagating the electron ray from the reference plane back to the radiator at z=-100m, and then forward again along the photon momentum direction. The random seeds in <random> are saved so the event can be generated again by replay.

GENBEAM 'postcol'

This is the same as the 'precol' GENBEAM option, except that only beam photons above the minimum energy specified in the BEAM card that make it all the way through the beam collimation system to a second reference plane at z=-18m are included in the output. The event count specified in the TRIG card is the total number of radiated photons, so to get 1000 events in the output file with GENBEAM 'postcol', the control.in file might need to have TRIG 4000. As with GENBEAM 'precol', there are no detector hits in the output hddm file, the only content in the output events is the <reaction> tag. The content and meaning of the <reaction> tag is the same for GENBEAM 'precol' and GENBEAM 'postcol', the only difference is that 'precol' includes all photons produced at the radiator in the output file, while 'postcol' only includes those that make it through the secondary collimator without being absorbed or scattering.

Monte Carlo beam generation in hdgeant4

The hdgeant4 simulation supports the 'precol' and 'postcol' GENBEAM options in control.in the same as hdgeant, but in addition it also provides two other GENBEAM options, 'postconv' and 'BHgen'. When 'postconv' is active, the built-in photon beam generator configured according to the standard BEAM options is used to generate coherent bremsstrahlung photons, which are propagated down the beamline through the collimators to the position of the TPOL converter. All photons that stop before reaching the TPOL converter are dropped from the output record, but those that reach the converter are forced to undergo conversion at a random depth within the converter target.

GENBEAM 'postconv'

The QED differential cross section with full polarization dependence is applied in the Monte Carlo simulation of the photon pair conversion process, which takes place either coherently from an atom in the converter (coherent pair production) or else incoherently from a single electron (triplet production). An example of the latter is shown in the following output event produced by hdgeant4 with GENBEAM 'postconv' activated in the control.in file. 

    <reaction type="0" weight="0">
      <beam type="Photon">
        <momentum E="10.835" px="-0.000177278" py="-5.44413e-05" pz="10.835" />
        <polarization Px="0.0011168" Py="0.00412377" Pz="3.89926e-08" />
      </beam>
      <vertex>
        <product decayVertex="0" id="1" mech="0" parentid="0" pdgtype="11" type="Electron">
          <momentum E="11.9997" px="-0.000117054" py="1.90218e-05" pz="11.9997" />
        </product>
        <origin t="-74.2076" vx="0.137271" vy="0.039639" vz="-2400" />
      </vertex>
      <vertex>
        <product decayVertex="1" id="2" mech="1447972675" parentid="1" pdgtype="-11" type="Positron">
          <momentum E="4.13102" px="0.00145518" py="0.000154387" pz="4.13102" />
        </product>
        <product decayVertex="1" id="3" mech="1447972675" parentid="1" pdgtype="11" type="Electron">
          <momentum E="6.70397" px="-0.00179934" py="-0.000355582" pz="6.70397" />
        </product>
        <product decayVertex="1" id="4" mech="1447972675" parentid="1" pdgtype="11" type="Electron">
          <momentum E="0.000559377" px="0.000166887" py="0.000146754" pz="4.89378e-05" />
        </product>
        <origin t="-44.7902" vx="0.0726302" vy="-0.0150264" vz="-1518.09" />
      </vertex>
      <random seed1="1955466792" seed2="75289136" seed3="709975949" seed4="912931185" />
    </reaction>

The leading <beam> photon and radiator <vertex> elements are the same as described above under GENBEAM 'precol'. The second <vertex> element describes a triplet conversion that takes place inside the TPOL converter. The first electron listed is the forward-going one of the pair, while the second is the recoiling atomic electron from the triplet process. In the case of a coherent pair production event, the second electron <product> tag is omitted from the second <vertex>. The weight parameter in the <reaction> tag for GENBEAM 'postconv' simulations is always set to zero and can be ignored; these events are already weighted by the differential cross section at the generation stage, and do not need to be assigned different weights at simulation time. Instructions for how to use the GENBEAM 'postconv' option to estimate the trigger efficiency and analyzing power of the TPOL are given in HOWTO perform a triplet polarimeter simulation.

GENBEAM 'BHgen'

The 'BHgen' option is identical in many ways to the 'postconv' option described above, except that in this case the conversion event takes place inside the liquid hydrogen target inside the GlueX solenoid, instead of the TPOL conversion target. The QED differential cross section with full polarization dependence is applied in the Monte Carlo simulation of the photon pair conversion process, which takes place either coherently from a hydrogen atom in the target (coherent pair production) or else incoherently from a single electron (triplet production). An example of the former is shown in the following output event produced by hdgeant4 with GENBEAM 'BHgen' activated in the control.in file. The example shows a case of coherent pair production where the energy transfer to the target is large enough to produce a recoil proton track. Even though the scattering is from the entire atomic charge distribution, the probability that the electron remains bound to the scattered proton in the final state is negligible, resulting in the full recoil momentum being assigned to the final-state proton in the simulation. Equivalently, screening by the electron at the corresponding value for Q2 is negligible, so coherent scattering from the atom at high Q2 is essentially incoherent scattering from the proton. Hadronic excitations of the target, ie. interferences between the Bethe-Heitler and the DVCS process, are not taken into account in this QED calculation. 

    <reaction type="221" weight="0.000802113">
      <beam type="Photon">
        <momentum E="9.1424" px="-6.1159e-05" py="-0.000116298" pz="9.1424" />
        <polarization Px="-0.0389073" Py="0.0241197" Pz="4.65457e-08" />
      </beam>
      <vertex>
        <product decayVertex="0" id="1" mech="0" parentid="0" pdgtype="11" type="Electron">
          <momentum E="12.0002" px="3.66898e-05" py="3.34213e-05" pz="12.0002" />
        </product>
        <origin t="-102.264" vx="-0.0201677" vy="0.000241349" vz="-2400" />
      </vertex>
      <vertex>
        <product decayVertex="1" id="2" mech="1447972675" parentid="1" pdgtype="-11" type="Positron">
          <momentum E="0.637599" px="0.501584" py="-0.391354" pz="0.0422823" />
        </product>
        <product decayVertex="1" id="3" mech="1447972675" parentid="1" pdgtype="11" type="Electron">
          <momentum E="7.5973" px="-0.0888835" py="0.0834119" pz="7.59632" />
        </product>
        <product decayVertex="1" id="4" mech="1447972675" parentid="1" pdgtype="2212" type="Proton">
          <momentum E="1.84578" px="-0.412761" py="0.307825" pz="1.5038" />
        </product>
        <origin t="-20.3886" vx="-0.110665" vy="-0.148827" vz="54.5567" />
      </vertex>
      <random seed1="972847607" seed2="497991566" seed3="709976015" seed4="912931251" />
    </reaction>

The big difference with GENBEAM 'BHgen' versus the previous GENBEAM modes is seen in the appearance of the non-zero value for the "weight" parameter, as seen in the above <reaction> tag. This weight must be carried with each event in the generated sample, and used as a weight factor in all histograms that are produced in the analysis of the samples. Below is a set of recurring questions that arise in the meaning and use of these weighted samples.

Why do I see fewer events in the hdgeant4 output after running with GENBEAM 'BHgen' that were requested in my control.in?
The 'BHgen' event generation process is broken down into several stages, beginning with sampling a photon from the coherent bremsstrahlung spectrum emitted from the radiator, followed by the propagation of the photon through the collimator and beamline down to the liquid hydrogen target, and terminating with the conversion of the beam photon in the target. Not all photons emitted from the radiator make it all the way to the liquid hydrogen target. Most of the losses occur at the primary collimator. For example, with a 12 GeV electron beam on a diamond radiator with the coherent edge at 9 GeV, only about 25% of the radiated photons coming from the radiator make it through the primary collimator. Only those that survive all the way to the GlueX target without being scattered or absorbed are able to produce output events in a GENBEAM 'BHgen' run. For all intents and purposes, the actual number of events in the output sample should be considered the generated sample size, not the value that was originally requested in the control.in TRIG card.
What is the meaning of events in the hdgeant4 output with weight="0"? Can I just drop them?
If you look in the output from a hdgeant4 run with GENBEAM 'BHgen', you see a substantial fraction of the output events are marked with weight="0". There is no use simulating and reconstructing these events with weight="0", because they will not contribute any weight to any final signals emerging from an analysis. However, there are reasons why they were included in the output sample, one technical and one practical. The technical reason why these zero-weight events are produced by the generator lies in the importance sampling scheme employed to generate final-state kinematics with a pre-weighted distribution that approximates the strongly peaking behavior of the physical process. Because these peaks lie close to the edges of the kinematically allowed phase space, the generated values sometimes appear in corners of the phase space plot that are disallowed by energy-momentum conservation. When this happens, the corresponding events are written out, but they are assigned zero weight. Why write them out at all? The answer to this question raises the practical use that these zero-weight events have in the interpretation of the weighted sample. The meaning of the event weight is that its average value over the output Monte Carlo sample gives the integrated differential cross section in microbarns for the pair conversion process. This is described in more detail in a following question. The practical reason these events are included in the output sample is that they are needed to compute the average event weight over the Monte Carlo sample, which has this very specific meaning in the context of GENBEAM 'BHgen' output. If you don't care about knowing the cross section and only want a weighted sample of simulated events, these zero-weight events can simply be ignored and discarded.
Why do all of my output events with weight > 0 have a recoil proton? Isn't production from the target electrons also included in the generator?
In principle, both coherent (scattering from the proton, with screening corrections at low Q2 from the electron cloud) and incoherent (scattering from the atomic electron) are simulated, weighted according to their respective cross sections from leading-order QED. However, depending on the threshold pair mass that you set in the third field of the line "GENBEAM 'BHgen' [threshold_pair_mass]" in your control.in, the incoherent (triplet) process may be excluded due to kinematics. The maximum pair mass that is allowed by energy-momentum conservation in the triplet process is [2 me k]1/2 for a beam photon energy k. This generally means that pairs from the triplet process are too low in invariant mass to make it into the acceptance of the forward calorimeters.
What is the physical meaning of the magnitude of the weight, or is it only a relative factor?
The scale of the event weights from GENBEAM 'BHgen' is such that its average value, in the limit of a large sample size, equals the integrated differential cross section in microbarns for the pair conversion process of a high energy photon incident on a hydrogen atom target. Keep in mind that the output from a GENBEAM 'BHgen' run generally contains a mixture of many beam photon energies at once, so the average weight must be computed separately for each bin in photon beam energy in order for it to make sense in this way. The reason for stating "integrated differential cross section" and not "total cross section" is that the domain of integration is flexible and can be chosen by the user. To get the total cross section as a function of beam energy, set the pair mass threshold to zero, as in "GENBEAM 'BHgen' 0.0" (or simply "GENBEAM 'BHgen'") and average the weights over all output events, binned by beam energy. Likewise, to get the partial cross section for e+e- pair production above a pair mass threshold of 1.2 GeV, repeat the same run with the pair mass threshold specified in control.in as "GENBEAM 'BHgen' 1.2". To select only the triplet process, the same average per beam energy bin should be computed, while setting the weights of events with a recoil proton to zero, and conversely for the Bethe Heitler process. Whatever selection is being performed to narrow the final states of interest, the same complete set of all sample events including those with weight="0" must be counted in computing the denominator of the average. As a consequence of these rules, the scale of the weight assigned to any given final state changes, depending on the pair mass threshold for the sample to which it belongs.
How can I convert this weighted event sample to an unweighted sample, where all events have the same weight?
If you simply carry the generator weights along with the event through the reconstruction and analysis stages, and include them as a factor in the weight assigned each time a bin in a counting histogram is incremented for this event, the proper statistical error propagation for weighted events will be automatically taken care of. However, there are some cases where it is more convenient to work with an unweighted sample where all events are weighted equally and Poisson counting statistics apply. Instructions and an example script called accept.py that can be used to convert a weighted sample into an equivalent unweighted one are provided near the end of the wiki page, HOWTO generate Bethe-Heitler pairs in the GlueX target.
How do I figure out how many generated events I need to end up with an unweighted sample of a desired size?
The plots at the end of HOWTO generate Bethe-Heitler pairs in the GlueX target provide a means to calculate the number of weighted events you will need to end up with a given number of unweighted events, under nominal operating conditions of the GlueX beam line. But the most reliable way to answer this question is to generate a small sample of weighted events, run them through accept.py and count the number of output events. To be sure you are operating in the large-N limit, repeat the test again with twice the event size and verify that the input:output ratio is stable. In general, it is advantageous to set the minimum pair mass threshold on the GENBEAM line in control.in as high as possible for the given reconstruction analysis cuts, as long as this can be done without losing events at the low-mass end of the spectrum.

GENBEAM event generation rates

The following table shows the number and rate of output events from GENBEAM runs under various options. The rate measurements were taken on a 3.2 GHz Intel Xeon processor.

Type Triggers Nout Nout/cpu/sec
precol 10,000 10,000 300
postcol 10,000 2,411 280
preconv 10,000 2,411 250
postconv 10,000 2,411 15
BHgen 10,000 2,411 160
Retrieved from "https://halldweb.jlab.org/wiki/index.php?title=HOWTO_use_hdgeant4_as_a_beam_Monte_Carlo_generator&oldid=105413"