Changes between Initial Version and Version 1 of SysCalc

Timestamp:

Nov 22, 2013, 2:56:34 PM (12 years ago)

Author:

Olivier Mattelaer

Comment:

--

SysCalc

@@ +1 @@
+=======================================================
+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