What is SysCalc

==Running the code

./sys_calc inputfile configfile outfile

1. inputfile is a valid input file (typically a lhe file generated be MG5_aMC with the option use_syst=T
2. configfile is a file as described below
3. outputfile is the place where to write the output file.

Standalone Interface

parameter	meaning	options
scalefact	various value used for the variation of the factorization and renormalization scale
scalecorrelation	Define how the correlation between the factorization and renormalization scale has to be handle	-1: all combination
		-2 only correlated ones
		a list of number of type (i*N+j) will select this list with i: index of the factorization and j the index of the renormalization scale (both starts at zero)
alpsfact	scale variation for the first emission in MLM matching
matchscale	QCUT variation of Pythia	QCUT should be set to the minimal value (default in madevent) to have this parameter meaningful.
PDF	pdf variation

This is an example of configuration file

# Central scale factors
scalefact:
0.5 1 2
# Scale correlation
# Special value -1: all combination (N**2)
# Special value -2: only correlated variation
# Otherwise list of index N*fac_index + ren_index #index starts at 0
scalecorrelation:
0
# \alpha_s emission scale factors
alpsfact:
0.5 1 2
# matching scales
matchscale:
32 50 80
# PDF sets and number of members (optional)
PDF:
CT10.LHgrid 53
MSTW2008nlo68cl.LHgrid

MadEvent Interface

As all external module, SysCalc can be called either in the main flow or in a a-posteriori way via the command "syscalc" in the command interface (./bin/madevent or launch -i from the MG5 interface)

In both case, the information from the run are used from the run_card and require to have use_syst on the run_card set on T. All the parameter correspond to those describe above.

Concerning the PDF since the multiple line is not supported in the run_card, you can separate the various PDF line information by '&&'

WARNING: If you run with use_syst=T with matching the output of pythia6 will not not have performed the (complete) veto of the events. This allows SysCalc to compute the final veto and assign it has a 0/1 weight. But it means that you can not use the output of pythia6 without running SysCalc. Also all the plot should then take care of those weights which is not currently supported by MadAnalysis.

Output

at parton level

• LHEF version 3 format containing event by event weight

at pythia level

• a second file containing the event by event weight.

Note that SysCalc prints also the sum of the weight on the screen [corresponding to the cross-section and/or of the effective number of events depending how the weight of the original file are defined].

Technical Details

=======================================================
Description of variables stored for systematics studies
=======================================================
By: Johan Alwall, 28/3/2012

-------------------------------------------------------------------
Turn on systematics info with the flag use_syst in the run_card.dat
Note that systematics only works with matching (ickkw set to 1)
-------------------------------------------------------------------

Parameters that can be varied after-the-fact (without need to rerun
Pythia+detector simulation):

- Central renormalization scale
- Central factorization scale
- PDF choice
- Emission renormalization scale factor
- PDF reweighting scale factor (not available at present)
- QCUT scale

=============================================
Variational parameters in reweight.f:
=============================================

*********************************************
Central scale (ren scale):
*********************************************
line 702: scale (asref = alphas(scale))

Event weight given by:
----------------------
alpha_s weight: alphas(scale)^N where 
                N=#(QCD vertices) - #(emission alpha_s vertices below)

*********************************************
Emission alpha_s reweighting:
*********************************************
line 912: qnow = sqrt(q2now)

Event weight given by:
----------------------
alphas(alpsfact*qnow)

*********************************************
PDF reweighting:
*********************************************
line 873  (initial pdf): ipdgcl(idacl(n,i)), xnow(j), q2now (pdgini,xini, q2ini)
line 1070 (cont.   pdf): ipdgcl(idacl(n,i)), xnow(j), q2now (pdgint,xint, q2int)
                         etc. (for both sides (1,2))

Event weight given by:
----------------------
initial pdf(pdgini,xini,q2ini)
*pdf(pdgint,xint,q2int)/pdf(pdgint,xint,q2ini)
... etc.
Note: Central fact scale variation corresponds to reweighting the last 
      scale only on each side

=============================================
Variation of QCUT in ME2pythia.f:
=============================================

failing criteria in parentheses. line numbers approximate.

SHOWERKT: QCUT
line 999: PTSORT(1) (in lhe file) (< QCUT)
line 1012 (non-highest mult): shower kt (> QCUT)
line 1025 (highest mult): shower kt (> PTSORT(1))

kT-MLM: YCUT=QCUT**2
line 1090: NJETS (< NLJETS)
           actually Y(NLJETS) (< YCUT)
line 1107 (non-highest mult): Y(NLJETS+1) (> YCUT)
if highest mult case: YCUT=PTSORT(1)
line 1133: Y(NN) (> YCUT)
line 1145: If not clustered, fail
line 1176: Y(2) (> YCUT)

So, just need three number for systematic variation of QCUT:
SMIN. Fail if < QCUT: 
   For SHOWERKT: PTSORT(1)
   For kT-MLM: Y(NLJETS)
SCOMP. Comparison number: 
   For highest mult: max(QCUT,PTSORT(1))
   Otherwise QCUT
   Perhaps use minimum safe QCUT or 0 instead of QCUT for systematics studies
   In any case, use max(QCUTcurr,comparison number) for arbitrary QCUTcurr.
SMAX. Fail if > comparison number:
   For SHOWERKT: shower kt
   For kt-MLM: max(Y(NLJETS+1),Y(NN),Y(2))

Note that some events will always fail - I suggest to simply ignore
those (as well as requiring minimum safe QCUT = xqcut for SHOWERKT and
max(xqcut+10,xqcut*1.3) for kT-MLM).

Event weight given by:
----------------------
1 if QCUT < SMIN and SMAX < max(QCUT, SCOMP)
otherwise 0


===============================================
Each line in the syst.dat file has the entries:
===============================================

<mgrwt event="event_num">
<rscale>n_qcd ren_scale</rscale>
<asrwt>n_alpsem alpsem_scale(1) ... alpsem_scale(n_alpsem)</asrwt>
<pdfrwt beam="1">n_pdfrw1 pdf_pdg_code1(1) ... pdf_pdg_code1(n_pdgrw1) \
 pdf_x1(1) ... pdf_x1(n_pdfrw1) pdf_q1(1) ... pdf_q1(n_pdfrw1)</pdfwgt>
<pdfrwt beam="2">n_pdfrw2 pdf_pdg_code2(1) ... pdf_pdg_code2(n_pdgrw2) \
 pdf_x2(1) ... pdf_x2(n_pdfrw2) pdf_q2(1) ... pdf_q2(n_pdfrw2)</pdfwgt>
<totfact>total_reweight_factor</totfact>
<matchscale>SMIN SCOMP SMAX</matchscale>
</mgrwt>

Total event weight for event event_num given by:
------------------------------------------------
alpha_s(scalefact*ren_scale)^(n_qcd) *  # central ren scale
alpha_s(alpsfact*alpsem_scale(1)) *              # emission ren scale
alpha_s(alpsfact*alpsem_scale(2)) *              # emission ren scale
...
pdf(pdf_pdg_code1(1),pdf_x1(1),pdf_q1(1))*       # initial state pdf
pdf(pdf_pdg_code1(2),pdf_x1(2),pdf_q1(2))/
pdf(pdf_pdg_code1(2),pdf_x1(2),pdf_q1(2))*       # pdf reweighting
pdf(pdf_pdg_code1(3),pdf_x1(3),pdf_q1(3))/
pdf(pdf_pdg_code1(3),pdf_x1(3),pdf_q1(2))*       # pdf reweighting
 ...
pdf(pdf_pdg_code1(n_pdfrw1),pdf_x1(n_pdfrw1),scalefact*pdf_q1(n_pdfrw1))/
pdf(pdf_pdg_code1(n_pdfrw1),pdf_x1(n_pdfrw1),pdf_q1(n_pdfrw1-1))*
# Note the central scale reweighting by scalefact above
# (if n_pdfrw1 = 1, need to reweight the initial state pdf scale)
# Also note that no scale should be larger than the last one (including
# scalefact) for that beam.
#
# Now the same thing for all pdfs in beam 2 (n_pdfrw2)
pdf(pdf_pdg_code2(1),pdf_x2(1),pdf_q2(1))*       # initial state pdf
 ...
pdf(pdf_pdg_code2(n_pdfrw2),pdf_x2(n_pdfrw2),scalefact*pdf_q2(n_pdfrw2))/
pdf(pdf_pdg_code2(n_pdfrw2),pdf_x2(n_pdfrw2),pdf_q2(n_pdfrw2-1))*
/ total_reweight_factor                          # corr. factor from MG run
* 0 if (QCUT > SMIN or SMAX > max(QCUT, SCOMP)), otherwise 1}}}