HOWTO generate Bethe-Heitler pairs in the GlueX target

From GlueXWiki
Jump to: navigation, search

The generation and simulation of Bethe-Heitler events in the GlueX spectrometer are both performed using the same hdgeant4 simulation program, by the following three-step procedure.

  1. Make sure that your hdgeant4 build has the dirac++ library enable if it is not the case do the following:
    • Download and build the Dirac++ library from https://github.com/rjones30/Diracxx.git
    • Define a new environment variable DIRACXX_HOME pointing to the directory you just installed above and do not forget to add DIRACXX_HOME in LD_LIBRARY_PATH
      • export DIRACXX_HOME=path to /Diracxx
      • export LD_LIBRARY_PATH=$DIRACXX_HOME:$LD_LIBRARY_PATH
    • Go into your HDGeant4 build directory and do "make clean; make" to rebuild hdgeant4 including Diracxx support
  2. A sample of beam photons from the GlueX radiator is generated over a limited spectral range, typically just the energy range of the tagger microscope and above. The coherent bremsstrahlung photon beam generator built into hdgeant4 is used for this. It gives a physically realistic photon transverse beam intensity profile incident on the GlueX target, with the proper energy and polarization spectrum. The correlation between energy and transverse position that is characteristic of coherent bremsstrahlung is correctly simulated. Using the following special configuration of control.in, every beam photon entering the liquid hydrogen target can be forced to undergo the Bethe Heitler process at a random location along the length of the target. The process is treated in a fashion fully consistent with tree-level QED, taking into account the polarization of the photon and both space-like and time-like form factors of the proton target. Internal radiative corrections are not currently included, but external radiation is automatically taken into account by the Geant4 tracking library. Check the following lines in control.in to make sure they are properly configured, then generate your Bethe-Heitler sample using the simple command "hdgeant4".
    • INFIle 'xxx' - make sure this line is commented out
    • TRIG 1000000 - set this to the number of Bethe-Heitler events you want to generate
    • RUNG 51500 - or whatever run you want to simulate
    • BEAM 12. 9.0 8.0 76.00 0.005 4.e-9 50.e-6 0.5e-3 - check the comments in control.in for the meaning of these fields
    • GENBEAM 'BHgen' 2.156 - very important, make sure this line is there and not commented out, and replace 2.156 with the minimum pair invariant mass that you want to generate, or leave it blank to generate all pair masses.
    • OUTFile 'my_BHgen.hddm' - or whatever you want to name the output file
    • RNDM - if you are generating your sample in multiple batch jobs, make this unique for each job
    • BGRATE - make sure this line is commented out
    • TRAJ - make sure this line is commented out to prevent explosive growth of your output file
    • CUTS - make sure this line is commented out to prevent the pairs/triplets from being pruned
    • ... - the remaining lines don't matter much, leave the defaults
  3. The sample of pairs that were generated by GENBEAM 'BHgen' in step 1 now needs to be pruned. First of all, there are a significant number of pairs with a low energy electron or positron that is outside the GlueX acceptance. Second, the events are generated with unequal weights, and you probably want to implement an accept-reject filter to equalize the weights before you simulate them. Third, you probably want to look them over before investing in all of the cpu cycles that are needed to simulate and analyze your sample. Write your own ROOT macro to apply minimal cuts to satisfy the experimental trigger, implement the accept-reject filter, and make various plots to make sure they look like they should. You will write your filtered event sample to a new file, here named my_bhgen_filtered.hddm. Below you will find a couple of simple examples of ROOT macros to illustrate how to write your own filter.
  4. The filtered sample of generated pairs now needs to be simulated. For this, a new control.in configuration is needed. Once these edits are in place, simulate your generated sample with the command 'hdgeant4'. This will take a while, so you may want to subdivide your sample and run on the OSG.
    • INFIle 'my_bhgen_filtered.hddm' - make sure this line is present and enabled
    • TRIG 1000000000 - set this to greater than the number of events in your input file
    • RUNG 51500 - or whatever run you want to simulate
    • BEAM 12. 9.0 8.0 76.00 0.005 4.e-9 50.e-6 0.5e-3 - check the comments in control.in for the meaning of these fields
    • BGRATE - if you want to simulate out-of-time tags, see associated comments in control.in for how to set this
    • BGGATE - ditto
    • GENBEAM 'BHgen' 2.156 - very important, make sure this line is NOT there or IS commented out
    • OUTFile 'my_bhgen_sim.hddm' - or whatever you want to name the output file
    • RNDM - if you are generating your sample in multiple batch jobs, make this unique for each job
    • TRAJ - make sure this line is commented out to prevent explosive growth of your output file
    • CUTS - make sure this line is commented out to prevent the pairs/triplets from being pruned
    • ... - the remaining lines don't matter much, leave the defaults

Beam Polarization

By default, the internal photon beam generator in hdgeant/hdgeant4 is configured for a diamond radiator in the 90o orientation. The degree of linear polarization in the generated beam photons varies with photon energy, as dictated by the differential cross section for coherent bremsstrahlung. However, this may not be what you want if you are using the internal beam generator to produce Bethe-Heitler pairs, and using the simulation of these pairs to estimate the analyzing power of this process for a particular trigger and set of final state selection cuts. To measure the analyzing power, one wants to simulate photons with 100% polarization over the relevant region of the energy spectrum, and be able to do that for different polarization orientations. This capability has been added to hdgeant4 through a special command /PhotonBeam/enableFixedPolarization. Here is a sample session that illustrates how this command might be used to generate 10M events with 100% polarization along the 45o orientation, followed by another 10M events with 100% polarization along the -45o direction.

$ cat special_run.mac
/Photonbeam/enableFixedPolarization 1.0 45
/run/beamOn 100000000
/Photonbeam/enableFixedPolarization 1.0 135
/run/beamOn 100000000
exit
$ hdgeant4 -t4 special_run.mac

Analysis

Analysis of the simulation results proceeds in the usual way, passing them through mcsmear to hd_root to your analysis TSelector. Instructions for all of that are beyond the scope of this document.

Sample plots

The total Bethe-Heitler cross section vs photon energy is obtained by plotting the average event weight in each energy bin. The plot on the left shows this result before the accept-reject filter, and on the right after the filter. After the filter, all of the events have the same weight so the filter effectively converts the weight into a convenient marker of the total Bethe-Heitler cross section at the specific photon energy for that event. The average cross section over the range 3-12 GeV is 13613 +/- 5 microbarns. The error bars shown in the plots are statistical. Note that the error bars show the modulation of the beam spectrum intensity across the coherent edge, but the estimated cross section is smooth over the coherent peak, showing that the sampled pair and incident flux spectra are proportional. These plots were produced from a run of 100M beam photons, about 1/4 of which make it through the collimator. The weighted plot contains 24M events (left) whereas the one on the right after accept-reject contains only 150k events.

BHtotalXS weighted.png BHtotalXS unweightedC.png

Next I repeat the same experiment, but with a minimum pair mass of 1 GeV. Both of the plots below show weighted events, the left showing 1 GeV cut on a sample that was generated with no minimum pair mass in the generator, and the right showing an uncut sample generated with a 1 GeV minimum pair mass cut in the generator. The left plot was made starting from a sample of 100M beam photons, 24M Bethe-Heitler pairs. Notice how large the statistical error is, even for such a large sample, once a cut on the invariant mass of the pair is applied. The plot on the right is made starting from a similar size sample, but with the 1 GeV mass cut applied at the generator level. Comparison of these two plots of the same observable illustrates the power of applying the mass cut at the generator stage for improving the statistical accuracy at high pair mass. The mean partial cross section in the right-hand plot is 0.0073 +/- 0.0003 ub.

BHpartialXS masscut.png BHpartialXS uncutdB.png

The next two plots show the same total Bethe-Heitler cross section for pairs > 1GeV/c2, with increasing sample size.

BHpartialXS x500.png BHpartialXS x5000.png

The plots below show the transverse and longitudinal distributions of the Bethe Heitler pair production vertex, showing that it has the expected shape and size corresponding to the overlap region between the beam and the liquid hydrogen target.

BHvspot trans.png BHvspot long.png

Example macros

  • accept.py - implements an accept-reject algorithm to equalize the weights of all events for a given photon energy.
  • masscut.py - pass over the results of a BHgen generator run and histogram the weights vs photon energy of all events with a pair invariant mass above some threshold value.
  • yields.py - pass over the results of a BHgen generator run and histogram the weights vs photon energy of all events.
  • vspot.py - pass over the results of a BHgen generator run and histogram the transverse and longitudinal distribution of the generated Bethe-Heitler vertices.