ExampleEEWW: eeww_polarized.sin

# $Id: eeww_polarized.sin 3070 2011-03-28 08:09:25Z jr_reuter $
# Use WHIZARD's initial and final state polarization facilities to
# examine the angular distributions for all 36 helicity combinations
# in e+, e- -> W+, W-
########################################################################
#
# Copyright (C) 1999-2011 by 
#     Wolfgang Kilian <kilian@physik.uni-siegen.de>
#     Thorsten Ohl <ohl@physik.uni-wuerzburg.de>
#     Juergen Reuter <juergen.reuter@desy.de>
#     Christian Speckner <Christian.Speckner@physik.uni-freiburg.de>
#
# WHIZARD is free software; you can redistribute it and/or modify it
# under the terms of the GNU General Public License as published by 
# the Free Software Foundation; either version 2, or (at your option)
# any later version.
#
# WHIZARD is distributed in the hope that it will be useful, but
# WITHOUT ANY WARRANTY; without even the implied warranty of
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the 
# GNU General Public License for more details.
#
# You should have received a copy of the GNU General Public License
# along with this program; if not, write to the Free Software
# Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
#
########################################################################


# Define and compile the process
model = SM
process eeww = "e+", "e-" => "W+", "W-"
compile

# Define the beams and event count (for simulation)
sqrts = 200 GeV
beams = "e+", "e-"
n_events = 100000

# Tell WHIZARD to retain the W polarization in the generated events
polarized "W+", "W-"

# Loop over e+ / e- helicity
scan int hel_ep = (-1, 1) {
  scan int hel_em = (-1, 1) {
    # Apply the helicity setup
    beam_polarization = diagonal_density (hel_ep:1), diagonal_density (hel_em:1)

    # (Re)calculate the integral for each initial state polarization
    integrate (eeww)

    # Define histograms for all different final state helicies. We
    # have to automatically generate tags from the initial state
    # helicites
    string $hist_mm = sprintf "cta%i_%i_-1_-1" (hel_ep, hel_em)
    histogram $hist_mm (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i -1 -1" (hel_ep, hel_em)
    }
    string $hist_ml = sprintf "cta%i_%i_-1_0" (hel_ep, hel_em)
    histogram $hist_ml (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i -1  0" (hel_ep, hel_em)
    }
    string $hist_mp = sprintf "cta%i_%i_-1_1" (hel_ep, hel_em)
    histogram $hist_mp (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i -1  1" (hel_ep, hel_em)
    }
    string $hist_lm = sprintf "cta%i_%i_0_-1" (hel_ep, hel_em)
    histogram $hist_lm (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i  0 -1" (hel_ep, hel_em)
    }
    string $hist_ll = sprintf "cta%i_%i_0_0" (hel_ep, hel_em)
    histogram $hist_ll (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i  0  0" (hel_ep, hel_em)
    }
    string $hist_lp = sprintf "cta%i_%i_0_1" (hel_ep, hel_em)
    histogram $hist_lp (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i  0  1" (hel_ep, hel_em)
    }
    string $hist_pm = sprintf "cta%i_%i_1_-1" (hel_ep, hel_em)
    histogram $hist_pm (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i  1 -1" (hel_ep, hel_em)
    }
    string $hist_pl = sprintf "cta%i_%i_1_0" (hel_ep, hel_em)
    histogram $hist_pl (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i  1  0" (hel_ep, hel_em)
    }
    string $hist_pp = sprintf "cta%i_%i_1_1" (hel_ep, hel_em)
    histogram $hist_pp (-1, 1, 0.1) {
      $title = sprintf "polarization: %i %i  1  1" (hel_ep, hel_em)
    }

    # The analysis setup fills the different histograms
    analysis =
      if (all Hel == -1 ["W+"] and all Hel == -1 ["W-"] ) then 
        record $hist_mm (eval cos (Theta) ["W+"])
      endif; if (all Hel == -1 ["W+"] and all Hel ==  0 ["W-"] ) then
        record $hist_ml (eval cos (Theta) ["W+"])
      endif; if (all Hel == -1 ["W+"] and all Hel ==  1 ["W-"] ) then
        record $hist_mp (eval cos (Theta) ["W+"])
      endif; if (all Hel ==  0 ["W+"] and all Hel == -1 ["W-"] ) then
        record $hist_lm (eval cos (Theta) ["W+"])
      endif; if (all Hel ==  0 ["W+"] and all Hel ==  0 ["W-"] ) then
        record $hist_ll (eval cos (Theta) ["W+"])
      endif; if (all Hel ==  0 ["W+"] and all Hel ==  1 ["W-"] ) then
        record $hist_lp (eval cos (Theta) ["W+"])
      endif; if (all Hel ==  1 ["W+"] and all Hel == -1 ["W-"] ) then
        record $hist_pm (eval cos (Theta) ["W+"])
      endif; if (all Hel ==  1 ["W+"] and all Hel ==  0 ["W-"] ) then
        record $hist_pl (eval cos (Theta) ["W+"])
      endif; if (all Hel ==  1 ["W+"] and all Hel ==  1 ["W-"] ) then
        record $hist_pp (eval cos (Theta) ["W+"])
      endif

    # Generate events (keep the helicity information). As this takes
    # some time, we request status information.
    simulate (eeww) {
      ?polarized_events = true
      checkpoint = 10000
    }
  }
}

# Write out the generated histograms in eeww_polarized.ps / eeww_polarized.pdf
compile_analysis { $out_file = "eeww_polarized.dat" }