REPORT
final version 6.09
1. SM Higgs boson phenomenology
1.2 Theoretical owerview
Here I'd like to sum up what I have learned on the Higgs mechanism and point out important (in my opinion) and interesting details.
Motivation
EW gauge bosons are known to be massive, but if we introduce their masses with brute force, the following problems occur:
 Lagrangian is not gauge invariant any more
 Theory becomes nonrenormalizable
 Theory lost its unitarity
How it works
That is why in the SM the more elegant mechanism is used, the mechanism of spontaneous symmetry breaking.
Introduce new scalar field and potential with nonzero vev
Add interaction terms between new scalar field and fields which should be massive
Introduce symmetry breaking and reparameterize the field (expand around vev)
In new parameterization, one obtains physical field H with a normal sign in front of the massive term and also massive terms for the gauge bosons, quarks and leptons.
1.2 Problems
#1 bottom mass
\mathcal{L}_Y^{\text{hadr}} = \frac{1}{\sqrt{2}} \left(v + H\right) \sum_{f=1}^n \left( h_D^f\bar{d}^f d^f + h_U^f\bar{u}^f u^f\right) \rightarrow m_b = \frac{v h_D^b}{\sqrt{2}} 
#2 the amplitude of pp \rightarrow ZH


2. Two Higgsdoublet model
2.1 Motivation
Although SM can give very precise predictions it is incomplete since it is not explaining, for example, baryon asymmetry, dark matter and hierarchy problem. That is why extensions are needed. One of the possible extensions is a Supersymmetric Standard Model which requires in its minimal case two Higgsdoublets. 2HDM is particularly interesting because:
 this is the minimum extension
 the model provides an additional source of CPviolation and therefore potentially can explain baryonic asymmetry in nature
2.2 Realization
The extension is built with use of the same logic as in Higgs sector of SM.
 We are introducing two Higgsdoublets and corresponding potential
 Introducing the spontaneous symmetry breaking which gives masses to gauge bosons and leptons
 In case of 2HDM 5 physical fields are left (2 charged + 2 scalars + 1 pseudoscalar)
The model contains two free parameters \tan \beta (\beta  relation angle to diagonalize mass matrixes) and M_A  the mass of pseudoscalar Higgs boson
3. Tools
3.1 In our work we are going to use the following computing tools:
 SusHi + PDFset to calculate cross section b\bar{b} \rightarrow A up to NNLO
 2HDMC to calculate decay cross section A \rightarrow Zh
All tools were successfully installed. Here you can find script to operate with SusHi (runSCRIPT. pl + necessary lib slharoutines.pm) out file sample (report.txt) and plot, showing the cross section dependence on the 2HDM free parametrs \tan \beta and M_A.
Cross section  Branching Ratio 

3.2 Corrected run files + output files
 run  script to run SusHi + 2HDMC for M_A = \{250, 400, 1000\} GeV in the following parameter diapasons: \tan \beta \in \left[0.5, 60\right] and \sin \left(\alpha  \beta\right) \in [1, 1]. The script produces SusHi and 2HDMC output files for all calls, grouped in the folders by the M_A value.
 getDATA.pl  script to extract data from the output files. The output of this particular script consists of the following files:
 CS_250.txt, CS_400.txt, CS_1000.txt  b\bar{b} \rightarrow A cross sections for M_A = \{250, 400, 1000\}
 width_250.txt, width_400.txt, width_1000.txt  total decay widths
 BR_250_5_5.txt, BR_400_5_5.txt, BR_1000_5_5.txt  branching ratios for A \rightarrow b\bar{b}
 BR_250_6_6.txt, BR_400_6_6.txt, BR_1000_6_6.txt  branching ratios for A \rightarrow t\bar{t}
 BR_250_23_25.txt, BR_400_23_25.txt, BR_1000_23_25.txt  branching ratios for A \rightarrow ZH_{125}
 reports.tar.gz  all output files
 slharoutines.pm  library to manipulate with SLHA  files, corrected (close file calls are added)
4. Study of Adecays
Three first columns correspond to the M_H, M_C \sim 3 \times 10^2\ \text{GeV}, the last  M_H, M_C \sim 3 \times 10^3\ \text{GeV}.
A \rightarrow ZH_{125}
At M_A = 250\ \text{GeV} topquark production (M_t = 172\ \text{GeV}) is prohibited by kinematic conditions while bottomquark production (M_b = 4.18\ \text{GeV}) as well as ZH_{125} production (M_Z = 91\ \text{GeV},\ M_h = 125\ \text{GeV}) are allowed. The picture we see can be understood by taking into account two facts:
 At M_A = 400\ \text{GeV} topquark production is allowed.
 At M_A = 1000\ \text{GeV} we see qualitatively the same picture as before, but with higher values for \Gamma_{\text{tot}}, which is also make sense since the higher mass we have the more decay channels occur.  We noticed that the total width of A \rightarrow Zh for the pseudoscalar mass M_A = 1000\ \text{GeV} is very large which is caused by the dominance of decays into the heavy Higgs and charged Higgs. For the masses \sim 10^3\ \text{Gev} picture is different. It has the same onionlike structure but smaller values for the decay width.  
From this graph one can clearly see that branching ration of A \rightarrow Zh is suppressed at \cos\ \left( \alpha  \beta\right) \approx 0 and enhanced for higher absolute values of \cos\ \left( \alpha  \beta\right). The decrease of branching ratio with \tan\ \beta rise is explained by the dominance of the A \rightarrow b\bar{b} process in this region.  The BR graph for M_A = 400\ \text{GeV} is differs from the previous one:
 On the last graph for M_A = 1000\ \text{GeV} one can see cosinesquaredlike dependence with minimum around the \cos\ \left( \alpha  \beta\right) \approx 0. Also the picture is almost independent of the \tan\ \beta value. This can be explained with the symmetry between A \rightarrow t\bar{t} and A \rightarrow b\bar{b} branching ratios which "cancel" each other (see graphs bellow). 
A \rightarrow b\bar{b}
A \rightarrow t\bar{t}
A \rightarrow W^+H^
for the light Higgs scenario (result is the same for A \rightarrow W^H^+)

5. Study of the BreitWigner approximation
For each mass I tried to find a couple of parameters [\tan\ \beta, \sin\ (\alpha \beta)], when \Gamma_\text{tot} not too large and BR not to small.
5.1 First variant of BreitWigner formula
Mathematica code  Breit_Wigner.nb
5.2 Comparison with the modified BreitWigner formula
Solid line corresponds to the results obtained by the first variant of the BreitWigner formula. Orange dots  to results of the modified formula. For M_A = 250\ \text{GeV} and M_A = 400\ \text{GeV} both formulas are in a good agreement, while for M_A = 1000\ \text{GeV} there is a shift of \sim 100\ \text{GeV}. 
 txt files with data for the modified formula  reports.tar.gz
 Mathematica code  Breit_Wigner_2.nb (the change of paths is needed)
5.2 Comparison with nlo and nnlo
Results for the first variant of the BreitWigner formula:
5.3 Study of the scale dependence
I have run Sushi for the following pares of scale parameters:
Renormalization scale  2  2  1  1  1  0.5  0.5 
Factorization scale  0.5  0.25  0.5  0.25  0.125  0.25  0.125 
and obtained the following diapasons for lo, nlo and nnlo:
leading order  next to the leading order 
next to the next leading order  all in the one graph 
I am still confused since nlo and nnlo domains do not intersect.
6. Kfactor
Basically Kfactor is a constant, as it was expected.
7. Analytical amplitude calculation with FeynArts+FormCalc
7.1 pp \rightarrow A \rightarrow Zh comparison with SusHi + 2HDMC
Mathematica notebook  ME_bbZh.nb
After squaring and summation over polarizations and color indexes I've got the following result for the amplitude:
(1)  M^2 = \frac{3\pi^2 \alpha_\text{fine str}^2 \times \cos (\beta\alpha)^2 \times M_b^2 \times Y_3 \times \left(2s^2\left(\frac{M_h^2 Y_3}{M_Z^2} + Y_3\right)  Y_3\left(s\left(\frac{M_h^4}{M_Z^2}  2M_h^2 +M_Z^2\right) + \frac{s^3}{M_Z^2}\right)\right)}{9\times2 \cos^2 \theta_W \times M_W^2 \sin^4 \theta_W \times \left((s  M_A^2)^2 + M_A^2 \Gamma_A^2\right)}. 
Herew the factor \frac{1}{9} appears because of the averaging over color indexes for two initial quarks. To derive the cross section one should apply the following expression (from Peskin's book)
d \sigma = \frac{1}{2s}\times \frac{\bar{k}_1}{16 \pi^2 \sqrt{s}}\times M(s)^2d \Omega, 
where s = (p_1 + p_2)^2 is the CM Energy squared and \bar{k}_1 is the magnitude of the momentum of one of the outgoing particles. In our case we do not have an angle dependence, so after angular integration
\sigma = \frac{1}{s}\times \frac{\bar{k}_1}{8 \pi \sqrt{s}}\times M(s)^2. 
To get expression for \bar{k} we use
s = \left[ \begin{pmatrix} E_Z\\ \bar{k} \end{pmatrix} + \begin{pmatrix} E_h\\ \bar{k} \end{pmatrix} \right]^2 = (E_Z + E_h)^2 = m_Z^2 + m_h^2 + 2\bar{k}^2 + 2\sqrt{m_Z^2 + \bar{k}^2}\sqrt{m_h^2 + \bar{k}^2} 
and obtain
\bar{k} \to \frac{\sqrt{2 s m_h^22 m_h^2 m_Z^2+m_h^42 s m_Z^2+m_Z^4+s^2}}{2 \sqrt{s}}. 
This leads to the answer for the cross section
\sigma(b\bar{b} \rightarrow Zh) = \frac{M^2}{8 \pi s \sqrt{s}} \times \sqrt{\frac{(M_Z^2  M_h^2 +s)^2}{4s}  M_Z^2}. 
In order to obtain the numerical answer we used \alpha_{\text{fine str}},\ M_W,\ M_Z,\ \theta_W on electroweak scale and M_b on the M_A. Also to calculate hadronic cross section \sigma(pp \rightarrow Zh) we performed the integration with PDFs:
\sigma (pp \rightarrow Zh) = 2\ \int_0^1 dx_1 dx_2\ \theta\left(\hat{s}x_1 x_2  (M_h + M_Z)^2\right)\ f_1(x_1, Q) f_2 (x_2, Q)\ \sigma (b\bar{b} \rightarrow Zh)(\hat{s}x_1 x_2), 
where Q is factorization scale chosen to be \frac{1}{4} M_A, \hat{s} is collier CM energy, theta function appears to satisfy kinematic conditions. In the table bellow you can find our results compared with answer of the BW formula.
M_A, \text{GeV}  \sigma_{\text{SusHi}}(M_A) \times \text{BR}, \text{pb}  2\int dx_1 dx_2 f_1(x_1) f_2(x_2) \sigma_{\text{anal}} (x_1 x_2 \hat{s}),  \text{Ratio} 

250  0.167587  0.16754  0.999718 
400  0.138187  0.135506  0.980601 
1000  0.00017508  0.00011373  0.64959 
As you can see the numerical results for 250 and 400 GeV are in a good correspondence with analytical calculation. Slight divination of the ratio from the unit most probably caused by the values of electroweak parameters.
7.2 s dependence
To get analytical result for \frac{d \sigma}{d \sqrt{s}} we performed variable transformation x_1 = \sqrt{\frac{s}{\hat{s}}} e^{y},\ x_2 = \sqrt{\frac{s}{\hat{s}}} e^{y} with Jacobian
J = \begin{vmatrix} \sqrt{\frac{s}{\hat{s}}}e^{y} && \sqrt{\frac{s}{\hat{s}}}e^{y} \\ \frac{1}{2\sqrt{\hat{s}s}}e^{y} && \frac{1}{2\sqrt{\hat{s}s}}e^{y} \end{vmatrix} = \frac{1}{\hat{s}} \hspace{10pt} \Rightarrow \hspace{10pt} dx_1 dx_2 = \frac{1}{\hat{s}}ds dy = \frac{2\sqrt{s}}{\hat{s}} d\sqrt{s} dy. 
This leads to the following expression for cross section derivative:
\frac{d\sigma(pp \rightarrow Zh)(s)}{d \sqrt{s}} = \frac{4\sqrt{s}}{\hat{s}}\sigma(b\bar{b} \rightarrow Zh)(s)\int_{\infty}^{\infty}dy f_1\left(\sqrt{\frac{s}{\hat{s}}} e^{y} \right) f_2\left(\sqrt{\frac{s}{\hat{s}}} e^{y}\right). 
Now we can compare the results obtained previously with analytical calculation.
Mathematica notebook  int.nb
b\bar{b} \rightarrow Zh
This process consists of the following contributions 
