=======================================================
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