diff --git a/Lectures_my/EMPP/Lecture1/images/dupa.png b/Lectures_my/EMPP/Lecture1/images/dupa.png
new file mode 100644
index 0000000..70db5f7
--- /dev/null
+++ b/Lectures_my/EMPP/Lecture1/images/dupa.png
Binary files differ
diff --git a/Lectures_my/EMPP/Lecture1/images/result_weight.png b/Lectures_my/EMPP/Lecture1/images/result_weight.png
new file mode 100644
index 0000000..afda7e1
--- /dev/null
+++ b/Lectures_my/EMPP/Lecture1/images/result_weight.png
Binary files differ
diff --git a/Lectures_my/EMPP/Lecture1/mchrzasz.tex b/Lectures_my/EMPP/Lecture1/mchrzasz.tex
index f0dacdd..d3fc7eb 100644
--- a/Lectures_my/EMPP/Lecture1/mchrzasz.tex
+++ b/Lectures_my/EMPP/Lecture1/mchrzasz.tex
@@ -393,7 +393,7 @@
 $\longmapsto$ Let's use the {\color{BrickRed}{statistical}} world and estimate the uncertainty of an integral in this case :)\\
 $\rightarrowtail$ A variance of a MC integral:
 \begin{align*}
-V(\hat{I}) = \dfrac{1}{n} = \dfrac{1}{n} \Big\lbrace E(f^2) - E^2(f) \Big\rbrace = \dfrac{1}{n} \Big\lbrace \dfrac{1}{b-a} \int_a^b f^2(x)dx - I^2 \Big\rbrace
+V(\hat{I}) = \dfrac{1}{n} \Big\lbrace E(f^2) - E^2(f) \Big\rbrace = \dfrac{1}{n} \Big\lbrace \dfrac{1}{b-a} \int_a^b f^2(x)dx - I^2 \Big\rbrace
 \end{align*}
     \begin{alertblock}{}                                                                                                                                                                                                                                                                                                
 $\looparrowright$ To calculate $V(\hat{I})$ one needs to know the value of $I$!
@@ -521,18 +521,264 @@
 \end{footnotesize}
 \end{frame}
 
+\begin{frame}\frametitle{Large Number Theorem}
+\begin{footnotesize}
+
+\end{footnotesize}
+\end{frame}
+
+
+\begin{frame}\frametitle{Crude Monte Carlo method of integration}
+\begin{footnotesize}
+$\Rrightarrow$ {\color{MidnightBlue}{Crude Monte Carlo method of integration is based on Large Number Theorem (LNT): }}\\
+\begin{align*}
+\dfrac{1}{N} \sum_{i=1}^N f(x_i) \xrightarrow{N\to \infty} \dfrac{1}{b-a}\int_a^b f(x)dx =E(f)
+\end{align*}
+$\Rrightarrow$ The standard deviation can be calculated:
+\begin{align*}
+\sigma = \dfrac{1}{\sqrt{N}} \sqrt{\Big[ E(f^2) -E^2(f)\Big] }
+\end{align*}
+
+$\Rrightarrow$ From LNT we have:
+\begin{align*}
+P= \int w(x) \rho(x) dx = \int_0^{\pi/2} (\frac{l}{L} \cos x ) \frac{2}{\pi} dx= \dfrac{2l}{\pi L}  \xrightarrow{N\to \infty} \frac{1}{N}\sum_{i=1}^N w(x_i)
+\end{align*}
+$\Rrightarrow$ Important comparison between ''Hit and mishit'' and Crude \mc~methods. One can analytically calculate:
+\begin{columns}
+
+\column{2.5in}
+\begin{align*}
+\hat{\sigma}^{{\rm{Crude}}} \simeq \dfrac{1.57}{\sqrt{N}}
+\end{align*}
+
+\column{2.5in}
+\begin{align*}
+\hat{\sigma}^{{\rm{Hit~and~mishit}}} \simeq \dfrac{2.37}{\sqrt{N}}
+\end{align*}
+
+
+\end{columns}
+
+$\Rrightarrow$ Crude \mc~is \textbf{always} better then ''Hit and mishit'' method.
+
+
+\end{footnotesize}
+\end{frame}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+
+
+\begin{frame}\frametitle{Classical methods of variance reduction}
+\begin{footnotesize}
+
+$\Rrightarrow$ In Monte Carlo methods the statistical uncertainty is defined as:
+\begin{align*}
+\sigma = \dfrac{1}{\sqrt{N}}\sqrt{V(f)}
+\end{align*}
+$\Rrightarrow$ Obvious conclusion:
+\begin{itemize}
+\item To reduce the uncertainty one needs to increase $N$.\\
+$\rightrightarrows$ Slow convergence. In order to reduce the error by factor of 10 one needs to simulate factor of 100 more points!
+\end{itemize}
+$\Rrightarrow$ How ever the other handle ($V(f)$) can be changed! $\longrightarrow$ Lot's of theoretical effort goes into reducing this factor.\\
+$\Rrightarrow$ We will discuss {\color{Mahogany}{four}} classical methods of variance reduction:
+\begin{enumerate}
+\item Stratified sampling.
+\item Importance sampling.
+\item Control variates.
+\item Antithetic variates.
+\end{enumerate}
 
 
 
 
 
+\end{footnotesize}
+\end{frame}
 
 
 
 
+\begin{frame}\frametitle{Stratified sampling}
+\begin{footnotesize}
+$\Rrightarrow$ The most intuitive method of variance reduction. The idea behind it is to divide the function in different ranges and to use the Riemann integral property:
+\begin{align*}
+I = \int_0^1 f(u) du = \int_0^a f(u)du + \int_a^1 f(u) du,~ 0<a<1.
+\end{align*}
+
+$\Rrightarrow$ The reason for this method is that in smaller ranges the integration function is more flat. And it's trivial to see that the more flatter you get the smaller uncertainty. 
+$\rightrightarrows$ A constant function would have zero uncertainty!
+
+    \begin{exampleblock}{General schematic:}                                                                                                                                                                                                                                                                                                
+Let's take our integration domain and divide it in smaller domains. In the $j^{th}$ domain with the volume $w_j$ we simulate $n_j$ points from uniform distribution. We sum the function values in each of the simulated points for each of the domain. Finally we sum them with weights proportional to $w_i$ and anti-proportional to $n_i$. 
+        \end{exampleblock} 
+
+\end{footnotesize}
+\end{frame}
+
+\begin{frame}\frametitle{Stratified sampling - mathematical details}
+\begin{footnotesize}
+Let's define our integrals and domains:
+\begin{align*}
+I=\int_{\Omega} f(x) dx,{~}{~} \Omega=\bigcup_{i=1}^k w_i
+\end{align*}
+The integral over $j^{th}$ domain:
+\begin{align*}
+I_j=\int_{w_j} f(x) dx,{~}{~} \Rightarrow I = \sum_{j=1}^k I_i
+\end{align*}
+$\rightrightarrows$ $p_j$ uniform distribution in the $w_j$ domain: $dp_j=\frac{dx}{w_j}$.\\
+$\rightrightarrows$ The integral is calculated based on crude \mc~method. The estimator is equal:
+\begin{align*}
+\hat{I}_j = \frac{w_j}{n_j} \sum_{i=1}^{n_j}f(x_j^i)
+\end{align*}
+Now the total integral is just a sum:
+\begin{align*}
+\hat{I} = \sum_{j=1}^k \hat{I}_j = \sum_{j=1}^k \frac{w_j}{n_j} \sum_{i=1}^{n_j} f(x_j^{(i)}), 
+\end{align*}
+\begin{columns}
+\column{0.2in}
+\column{2in}
+Variance:
+$V(\hat{I})=\sum_{j=1}^k \dfrac{w_j^2}{n_j}V_j(f)$,
+\column{3in}
+and it's estimator:
+$\hat{V} (\hat{I}) = \sum_{j=1}^k \dfrac{w_j^2}{n_j} \hat{V}_j(f)$
+\end{columns}
+
+\end{footnotesize}
+\end{frame}
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{frame}\frametitle{Importance sampling}
+\begin{footnotesize}
+$\Rrightarrow$ If the function is changing rapidly in its domain one needs to use a more elegant method: make the function more stable.\\
+$\rightrightarrows$ The solution is from first course of mathematical analysis: change the integration variable :)
+\begin{align*}
+f(x)dx \longrightarrow \frac{f(x)}{g(x)}dG(x),~{\rm{where}}~g(x)=\frac{dG(x)}{dx}
+\end{align*}
+    \begin{exampleblock}{Schematic:}                                                                                                                                                                                                                                                                                                
+\begin{itemize}
+\item Generate the distribution from $G(x)$ instead of $\mathcal{U}$.
+\item For each generate point calculate the weight: $w(x)=\frac{f(x)}{g(x)}$.
+\item We calculate the expected value $\hat{E}(w)$ and its variance $\hat{V}_G(w)$ for the whole sample.
+\end{itemize}
+        \end{exampleblock} 
 
 
 
+\begin{itemize}
+
+\item If $g(x)$ is choose correctly the resulting variance can be much smaller.
+\item There are some mathematical requirements:
+\begin{itemize}
+\item $g(x)$ needs to be non-negative and analytically integrable on its domain.
+\item $G(x)$ invertible or there should be a direct generator of $g$ distribution.
+\end{itemize}
+\end{itemize}
+
+
+\end{footnotesize}
+\end{frame}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{frame}\frametitle{Importance sampling - Example}
+\begin{footnotesize}
+
+$\Rrightarrow$ Let's take our good old $\pi$ determination example.
+\vspace{0.3cm}
+\begin{columns}
+\column{0.1in}
+\column{3in}
+$\Rrightarrow$ Let's take here for simplicity: $L=l$. 
+\begin{itemize}
+\item Let's take a trivial linear weight function: $g(x)=\frac{4}{\pi}(1-\frac{2}{\pi}x)$
+\item It's invertible analytically: $G(x)=\frac{4}{\pi}x(1-\frac{x}{\pi})$
+\item The weight function:
+\begin{align*}
+w(x)=\frac{p(x)}{g(x)}=\frac{\pi}{4}\frac{\cos x}{1-2x/ \pi}
+\end{align*}
+\end{itemize}
+
+\column{2in}
+\includegraphics[width=0.95\textwidth]{images/result_weight.png}
+\end{columns}
+\begin{itemize}
+\item Now the new standard deviation is smaller:
+\end{itemize}
+\begin{align*}
+\sigma_{\pi}^{{\rm{IS}}} \simeq \frac{0.41}{\sqrt{N}} < \sigma_{\pi}\simeq \frac{1.52}{\sqrt{N}}
+\end{align*}
+\begin{itemize}
+\item Importance sampling has advantages:
+\begin{itemize}
+\item Big improvements of variance reduction.
+\item The only method that can cope with singularities.
+\end{itemize}
+\end{itemize}
+
+
+\end{footnotesize}
+\end{frame}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{frame}\frametitle{Control variates}
+\begin{footnotesize}
+$\Rrightarrow$ Control variates uses an other nice property of Riemann integral:
+\begin{align*}
+\int f(x) dx = \int [f(x)-g(x)]dx+ \int g(x)dx
+\end{align*}
+\begin{itemize}
+\item $g(x)$ needs to be analytically integrable.
+\item The uncertainty comes only from the integral: $\int [f(x)-g(x)]dx$. 
+\item Obviously: $V(f\to g) \xrightarrow{f\to g} 0$
+\end{itemize}
+$\Rrightarrow$ Advantages:\\
+\begin{itemize}
+\item Quite stable, immune to the singularities.
+\item $g(x)$ doesn't need to be invertible analytically.
+\end{itemize}
+$\Rrightarrow$ Disadvantage:\\
+\begin{itemize}
+\item Useful only if you know $\int g(x)dx$
+\end{itemize}
+
+\end{footnotesize}
+\end{frame}
+
+
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+
+\begin{frame}\frametitle{Antithetic variates}
+\begin{footnotesize}
+$\Rrightarrow$ In \mc~methods usually one uses the independent random variables. The Antithetic variates method on purpose uses a set of correlated variables (negative correlation is the important property):
+\begin{itemize}
+\item Let $f$ and $f\prime$ will be functions of x on the same domain.
+\item The variance: $V(f+f\prime)=V(f)+V(f\prime)+2 Cov(f,f\prime)$.
+\item If $Cov(f,f\prime)<0$ then you can reduce the variance.
+\end{itemize}
+$\Rrightarrow$ Advantages:
+\begin{itemize}
+\item If you can pick up $f$ and $f\prime$ so that they have negative correlation one can significantly reduce the variance!
+\end{itemize}
+
+$\Rrightarrow$ Disadvantages:
+\begin{itemize}
+\item There are no general methods to produce such a negative correlations.
+\item Hard to generalize this for multidimensional case.
+\item You can't generate events from $f(x)$ with this method.
+\end{itemize}
+
+
+
+\end{footnotesize}
+\end{frame}
+
 
 \backupbegin   
 
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails/ht.py b/Lectures_my/EMPP/MC_examples/Heads_tails/ht.py
new file mode 100644
index 0000000..fcb8078
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails/ht.py
@@ -0,0 +1,23 @@
+from ROOT import  * 
+import sys
+
+def main(argv=sys.argv):  
+
+
+    pi=3.14159265359
+    c=TCanvas("c1", "c1", 800,600)
+    fun=TF1("f1", "0.9*cos(x)", 0, pi/2)
+    fun.SetLineColor(kBlack)
+    fun.SetFillColor(1)
+    fun.SetTitle("")
+    fun.GetXaxis().SetTitle('x')  
+    fun.GetYaxis().SetTitle('y')
+
+    fun.Draw("FC")
+
+    c1.SaveAs("result.png")
+        
+
+
+if __name__ == "__main__":    
+    sys.exit(main())      
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails/ht.py~ b/Lectures_my/EMPP/MC_examples/Heads_tails/ht.py~
new file mode 100644
index 0000000..b9852a4
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails/ht.py~
@@ -0,0 +1,22 @@
+from ROOT import  * 
+import sys
+
+def main(argv=sys.argv):  
+
+
+    pi=3.14159265359
+    c=TCanvas("c1", "c1", 800,600)
+    fun=TF1("f1", "0.9*cos(x)", 0, pi/2)
+    fun.SetLineColor(kBlack)
+    fun.SetFillColor(1)
+    fun.SetTitle("")
+    fun.GetXaxis().SetTitle('x')  
+
+    fun.Draw("FC")
+
+    c1.SaveAs("result.png")
+        
+
+
+if __name__ == "__main__":    
+    sys.exit(main())      
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails/result.pdf b/Lectures_my/EMPP/MC_examples/Heads_tails/result.pdf
new file mode 100644
index 0000000..2a3cb7d
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails/result.pdf
Binary files differ
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails/result.png b/Lectures_my/EMPP/MC_examples/Heads_tails/result.png
new file mode 100644
index 0000000..3617873
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails/result.png
Binary files differ
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/ht.py b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/ht.py
new file mode 100644
index 0000000..d2d6407
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/ht.py
@@ -0,0 +1,54 @@
+from ROOT import  * 
+import sys
+
+def main(argv=sys.argv):  
+
+
+    pi=3.14159265359
+    c=TCanvas("c1", "c1", 800,600)
+ 
+
+    fun1=TF1("f2", "1.2732*(1-0.6366*x)", 0, pi/2) 
+    fun1.SetLineColor(kBlue)             
+    fun1.SetFillColor(1)                  
+    fun1.SetTitle("")                     
+    fun1.GetXaxis().SetTitle('x')         
+    fun1.GetYaxis().SetTitle('y')         
+                                     
+    fun1.Draw()
+
+
+
+    fun=TF1("f1", "0.9*cos(x)", 0, pi/2) 
+    fun.SetLineColor(kBlack)             
+    fun.SetFillColor(1)                  
+    fun.SetTitle("")                     
+    fun.GetXaxis().SetTitle('x')         
+    fun.GetYaxis().SetTitle('y')         
+    fun.Draw("SAME")  
+
+
+
+
+    fun2=TF1("f3", "0.785*cos(x)/(1-2*x/3.1415)", 0, pi/2)
+    fun2.SetLineColor(kRed)                        
+    fun2.SetFillColor(1)                            
+    fun2.SetTitle("")                               
+    fun2.GetXaxis().SetTitle('x')                   
+    fun2.GetYaxis().SetTitle('y')                   
+
+    fun2.Draw("SAME")
+
+    leg=TLegend(0.7, 0.6, 0.9, 0.8)
+    leg.AddEntry(fun1, "g(x)", "L")
+    leg.AddEntry(fun2, "w(x)", "L")
+    leg.AddEntry(fun, "p(x)", "L")
+
+
+    leg.Draw()
+    c1.SaveAs("result.png")
+        
+
+
+if __name__ == "__main__":    
+    sys.exit(main())      
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/ht.py~ b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/ht.py~
new file mode 100644
index 0000000..8cba83b
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/ht.py~
@@ -0,0 +1,52 @@
+from ROOT import  * 
+import sys
+
+def main(argv=sys.argv):  
+
+
+    pi=3.14159265359
+    c=TCanvas("c1", "c1", 800,600)
+ 
+
+    fun1=TF1("f2", "1.2732*(1-0.6366*x)", 0, pi/2) 
+    fun1.SetLineColor(kBlue)             
+    fun1.SetFillColor(1)                  
+    fun1.SetTitle("")                     
+    fun1.GetXaxis().SetTitle('x')         
+    fun1.GetYaxis().SetTitle('y')         
+                                     
+    fun1.Draw()
+
+
+
+    fun=TF1("f1", "0.9*cos(x)", 0, pi/2) 
+    fun.SetLineColor(kBlack)             
+    fun.SetFillColor(1)                  
+    fun.SetTitle("")                     
+    fun.GetXaxis().SetTitle('x')         
+    fun.GetYaxis().SetTitle('y')         
+    fun.Draw("SAME")  
+
+
+
+
+    fun2=TF1("f3", "0.785*cos(x)/(1-2*x/3.1415)", 0, pi/2)
+    fun2.SetLineColor(kRed)                        
+    fun2.SetFillColor(1)                            
+    fun2.SetTitle("")                               
+    fun2.GetXaxis().SetTitle('x')                   
+    fun2.GetYaxis().SetTitle('y')                   
+
+    fun2.Draw("SAME")
+
+    leg=TLegend(0.7, 0.6, 0.9, 0.8)
+    leg.AddEntry(fun1, "g(x)", "LEP")
+    
+
+    leg.Draw()
+    c1.SaveAs("result.png")
+        
+
+
+if __name__ == "__main__":    
+    sys.exit(main())      
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/result.pdf b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/result.pdf
new file mode 100644
index 0000000..2a3cb7d
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/result.pdf
Binary files differ
diff --git a/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/result.png b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/result.png
new file mode 100644
index 0000000..afda7e1
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Heads_tails_importance_sampling/result.png
Binary files differ
diff --git a/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/dupa.png b/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/dupa.png
new file mode 100644
index 0000000..70db5f7
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/dupa.png
Binary files differ
diff --git a/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/ht.py b/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/ht.py
new file mode 100644
index 0000000..bfa1be7
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/ht.py
@@ -0,0 +1,31 @@
+from ROOT import  * 
+import sys
+
+def main(argv=sys.argv):  
+
+
+    pi=3.14159265359
+    c=TCanvas("c1", "c1", 800,600)
+ 
+    NTOYS=1000000
+    rand=TRandom(123)    
+
+    hists=[]
+
+    for i in range(1,13, 2):
+        hist=TH1D("hist", "hist", 400, 0, 10)
+        for j in range(0, NTOYS):
+            tmp=0.
+            for k in range(0,i):
+                tmp+=rand.Rndm()
+            hist.Fill(tmp)
+        hists.append(hist)
+
+    for i in range(0,len(hists)): 
+        if(i==0): hists[i].DrawNormalized()
+        else: hists[i].DrawNormalized("SAME")
+
+    c.SaveAs("dupa.png")
+
+if __name__ == "__main__":    
+    sys.exit(main())      
diff --git a/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/ht.py~ b/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/ht.py~
new file mode 100644
index 0000000..d9f4fdf
--- /dev/null
+++ b/Lectures_my/EMPP/MC_examples/Large_Number_Theorem/ht.py~
@@ -0,0 +1,31 @@
+from ROOT import  * 
+import sys
+
+def main(argv=sys.argv):  
+
+
+    pi=3.14159265359
+    c=TCanvas("c1", "c1", 800,600)
+ 
+    NTOYS=10000
+    rand=TRandom(123)    
+
+    hists=[]
+
+    for i in range(1,13, 2):
+        hist=TH1D("hist", "hist", 400, 0, 10)
+        for j in range(0, NTOYS):
+            tmp=0.
+            for k in range(0,i):
+                tmp+=rand.Rndm()
+            hist.Fill(tmp)
+        hists.append(hist)
+
+    for i in range(0,len(hists)): 
+        if(i==0): hists[i].DrawNormalized()
+        else: hists[i].DrawNormalized("SAME")
+
+    c.SaveAs("dupa.png")
+
+if __name__ == "__main__":    
+    sys.exit(main())