HOPL

An annotated example
2.2.1 The BEGIN action
The BEGIN action is executed before the first event is processed. Typically it defines some variables and allocates histograms. There can be several BEGIN actions in a hepawk script, they will be concatenated.
This BEGIN action starts by printing some information on the Monte Carlo run which will be processed. The builtin function printf is similar to the function with the same name from the standard C library. See section A.8.2 for details.
# multiphoton.hepawk
BEGIN
{
printf ("\nKRONOS, %s:\n", REV);
printf ("Run: %d, Date: %s\n\n", RUN, DATE);
The builtin variables REV, RUN, and DATE are preset to the version and revison date of the Monte Carlo, the number of the run, and the date of the run respectively. This information is taken from /hepevt/. See section A.7 for more automatic variables.
Next we introduce variables as symbolic names for kinematical cuts, detector geometry, and histogramming parameters.
x_min = 0.001; x_max = 0.1; Q2_min = 100; # kinematical cuts
E_min = 0.5; E_max = 40; # histogramming range
n_max = 6; E_cut_max = 3; E_cut_step = 0.1; # multiplicity cuts
E_cut_bins = E_cut_max / E_cut_step;
th_min_lumi = 0; th_max_lumi = 0.5e-3; # luminosity monitors
th_max_main = PI - 0.100; # main detector (fwd)
th_min_main = 0.150; # (bwd)
th_sep_jet = 0.300; # separation from jet
th_sep_em = 0.150; # sep. from electron
Finally we allocate counters and histograms for the various event types and
observables we want to study. See section A.8.4 for more details on the interface
to the histogramming package.
incut_0 = 0; # counting 0 photon evts
incut_1 = 0; # counting 1 photon evts
h_E = book1 (0, "dsigma/dE", E_max, 0, E_max);
incut_2 = 0; # counting 2 photon evts
h_E1 = book1 (0, "dsigma/dE1", E_max, 0, E_max);
h_E2 = book1 (0, "dsigma/dE2", E_max, 0, E_max);
# Multiplicities
h_m = book2 (0, "Multiplicity", n_max, - 0.5, n_max - 0.5,
E_cut_bins, 0, E_cut_max);
h_E1_E2 = book2 (0, "dsigma/dE1dE2", # Lego plot
E_max, 0, E_max, E_max, 0, E_max);
}

2.2.2 Kinematics
The calculation of the basic kinematical variables should be almost self explaining. This action will be executed for all events, because no selection condition is given.
The basic data types in hepawk version 1.6 are real , vector , and particle. The lexical convention is to use identifiers starting with $ for vectors and identifiers starting with @ for particles. Components of these data structures are addressed by trailing selections, such as :p for the four momentum of a particle. Therefore the first line in the following action calculates the Mandelstam variable s from the four momenta of the two beam particles @B1 and @B2. See section A.2 for more details on the available data types.
{
S = (@B1:p + @B2:p)^2; # Center of mass energy
Next we loop over all outgoing electrons and store the momentum in $el.
for (@p in ELECTRONS) # Collect the outgoing electron.
$el = @p:p;
This method of obtaining the outgoing electron momentum is actually quite naive and uses the fact that the Monte Carlo KRONOS will generate one and only one electron. A more sophisticated example would have to use some physical criteria for distinguishing the scattered electron from electrons created in hadronic decays. An example would be to select the electron which is furthest away from any hadronic jet. Similar code selecting the hardest photon can be found in section 2.2.3.
In the next step the electron momentum is used to calculate the (electronic) momentum transfer and the usual Bjorken variables x, y.
$q = @B1:p - $el; # (Electronic) momentum transfer
Q2 = - $q^2;
x = Q2/(2*$q*@B2:p); # Bjorken variables
y = Q2/(x*S);
Finally the same naive method is used to determine the single outgoing quark momentum.
for (@p in QUARKS) # Collect the outgoing quark (= jet).
$jet = @p:p;
}
2.2.3 Event processing
After having determined the Bjorken variables in the last section, we can now use them for kinematical cuts. The next action will only be executed for those events that satisfy xmin. = x = xmax. and Q2 = Q2 min.. Note that (unlike other programming languages) comparisons can be combined to intervalls in a natural way. See section A.5 for more details on comparisons and logical operators.
(x_min <= x <= x_max) && (Q2_min <= Q2)
{
The first histogram will represent the total cross section inside the kinematical cuts as a function of the photon multiplicity and the infrared cut off in the photon energy. hepawk version 1.6 does not have a for loop similar to C (for is at the moment only used for looping over particles), thus we implement the loop over the cut offs with a while loop.
E_cut = 0.0;
while (E_cut <= E_cut_max)
{
For each cut off energy we now loop over all photons and select those that satisfy the energy cut and the angular cuts.
i = 0; # multiplicity counter
for (@p in PHOTONS)
if (@p:p:E >= E_cut)
{
A photon is counted if it falls either in the forward or backward luminosity monitor or in the main detector. Furthermore it has to be sufficiently separated from the outgoing electron and quark jet.
theta = angle (@p:p, @B1:p);
if ((th_min_lumi <= theta <= th_max_lumi # bwd
|| th_min_lumi <= PI - theta <= th_max_lumi # fwd
|| th_min_main <= theta <= th_max_main) # main
&& angle (@p:p, $quark) >= th_sep_jet
&& angle (@p:p, $el) >= th_sep_em)
i++;
}
Before reiterating the while loop, we fill the bin in the histogram corresponding to the measured multiplicity and the applied energy cut.
fill (h_m, i, E_cut + 0.01);
E_cut = E_cut + E_cut_step;
}
The other histograms record the energy spectrum of the hardest and second hardest photon in a multiphoton event. This is implemented by looping over all existing photons satisfying a certain energy cut and satisfying the same angular and separation cuts as above.
n = 0; $p_hard = $p_soft = $NULL;
for (@p in PHOTONS)
{
if (@p:p:E >= E_min)
{
theta = angle (@p:p, @B1:p);
if ((th_min_lumi <= theta <= th_max_lumi # bw lumi
|| th_min_lumi <= PI - theta <= th_max_lumi # fw lumi
|| th_min_main <= theta <= th_max_main) # main
&& (angle (@p:p, $quark) >= th_sep_jet)
&& (angle (@p:p, $el) >= th_sep_em))
{
If the current photon is harder than the second hardest photon seen in the
event so far, update the variables holding the momenta of these two photons.
if (@p:p:E > $p_hard:E) # the hardest photon yet
{
$p_soft = $p_hard; $p_hard = @p:p;
}
else if (@p:p:E > $p_soft:E) # the second hardest
$p_soft = @p:p;
n++;
}
}
}
Once we have decided which are the hardest photons in this event, we can
fill the corresponding histograms. Furthermore a counter is incremented which
will later be used to determine the accumulated cross section.
if (n == 0)
incut_0++; # count this event
else if (n == 1)
{
incut_1++; # count this event
fill (h_E, $p_hard:E);
}
else
{
incut_2++;
fill (h_E1, $p_hard:E); fill (h_E2, $p_soft:E);
fill (h_E1_E2, $p_hard:E, $p_soft:E);
}
}

2.2.4 The END action
The END action will be executed after the last event has been processed. Typically it is used to print results and to write histograms to a permanent storage medium. Multiple instances of the END action will be concatenated.
This END action starts by printing the integrated cross sections for zero, one, and multi photon events. These are constructed from the counters incremented earlier, the total number of scanned events, and the total cross section. The latter two variables are provided by the Monte Carlo in /hepevt/.
END
{
printf ("null photon events: %d, cross section: %gmb\n",
incut_0, (incut_0/NEVENT) * XSECT);
printf ("one photon events: %d, cross section: %gmb\n",
incut_1, (incut_1/NEVENT) * XSECT);
printf ("two photon events: %d, cross section: %gmb\n",
incut_2, (incut_2/NEVENT) * XSECT);
Now a line printer style plot of the histograms is printed out.
printf ("\nHISTOGRAMS:\n");
printf ("***********\n\n");
plot ();
And finally the contents of the histograms is normalized to picobarns per
bin and written to the file "PAW".
scale (1.0e9 * XSECT/NEVENT); # normalize to picobarns
save ("PAW");
printf ("\ndone.\n");
}

hepawk(ID:6208/hep001)