From d222eab9cbcc58c6bb83f43daf59262d9d57a25b Mon Sep 17 00:00:00 2001
From: Morfouace <morfouac@ipno.in2p3.fr>
Date: Tue, 12 Apr 2016 11:16:26 -0400
Subject: [PATCH] ! New EventGenerator to read an energy distribution and an !
 angular ditribution from an TH1F histogram to simulate p, d, t...

---
 Inputs/EventGenerator/3He.source             |   8 +-
 Inputs/EventGenerator/EnergyDist_proton.root | Bin 0 -> 5155 bytes
 Inputs/EventGenerator/ThetaDist_proton.root  | Bin 0 -> 6833 bytes
 Inputs/EventGenerator/alpha.source           |   6 +-
 Inputs/EventGenerator/deuton.source          |   4 +-
 Inputs/EventGenerator/proton.pbuu            |  14 +
 Inputs/EventGenerator/proton.source          |   2 +-
 Inputs/EventGenerator/triton.source          |   4 +-
 NPLib/Physics/NPReaction.h                   |  16 +-
 NPSimulation/Core/CMakeLists.txt             |   2 +-
 NPSimulation/Core/EventGeneratorIsotropic.cc |  18 +-
 NPSimulation/Core/EventGeneratorpBUU.cc      | 256 +++++++++++++++++++
 NPSimulation/Core/EventGeneratorpBUU.hh      |  66 +++++
 NPSimulation/Core/ParticleStack.cc           |   3 +-
 NPSimulation/Core/PrimaryGeneratorAction.cc  |  19 +-
 NPSimulation/PhysicsListOption.txt           |   6 +-
 Projects/Lassa/Analysis.cxx                  | 102 +++++---
 Projects/Lassa/Analysis.h                    |  15 +-
 Projects/Lassa/RunToTreat.txt                |   4 +-
 Projects/macros/GeometricalEfficiency.C      |   2 +-
 20 files changed, 473 insertions(+), 74 deletions(-)
 create mode 100644 Inputs/EventGenerator/EnergyDist_proton.root
 create mode 100644 Inputs/EventGenerator/ThetaDist_proton.root
 create mode 100644 Inputs/EventGenerator/proton.pbuu
 create mode 100644 NPSimulation/Core/EventGeneratorpBUU.cc
 create mode 100644 NPSimulation/Core/EventGeneratorpBUU.hh

diff --git a/Inputs/EventGenerator/3He.source b/Inputs/EventGenerator/3He.source
index 2c4bf8b61..d85a5bce8 100644
--- a/Inputs/EventGenerator/3He.source
+++ b/Inputs/EventGenerator/3He.source
@@ -1,13 +1,13 @@
-%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 %			   Energy are given in MeV , Position in mm				  %	
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
-	EnergyLow=  120	
-	EnergyHigh= 200
+	EnergyLow=  0	
+	EnergyHigh= 520
 	HalfOpenAngleMin= 0
-	HalfOpenAngleMax= 12
+	HalfOpenAngleMax= 90
 	x0= 0	
 	y0= 0	
 	z0= 0	
diff --git a/Inputs/EventGenerator/EnergyDist_proton.root b/Inputs/EventGenerator/EnergyDist_proton.root
new file mode 100644
index 0000000000000000000000000000000000000000..7c3d13b3d50b0dd2be68f12e171589a6ab4598e0
GIT binary patch
literal 5155
zcmbt&bx;)0*Y*n1-AGF#-67q$q>>BLy~M&Iy>v>W^a4@>QX-+C^s;obNOwvs2q=;-
z`kQ&@op0Wm@0<DV%$a-c+<DG(X6_%)xsSJ}rw;&d&<X$m>;V8uHUNNP?>?8gj|cam
z`0zgoPBs8QunKq}j)!>$;3wP8qA~BMW@|^^-~K<XV8Gu(+D0_A)c~P?mG9XA04zgw
z7dHoDu!Dz#w^M+Xn}GM<cK**Z0P8RA9~J=cx`!#=$3OgixBmaPrm(;K<;J!ChnM@u
zz4{mb&`?uY9S;g|u=DYL<t|L|-%YRtxj<a(<>Z86AZ=W4Y=AW`!2dy$VUSOf3H3?~
zYdm5@JxWabCs`%?+lQ5^-qMy!9qtD8>R-jqzr-|(ya?17^gr{OJSkYsd=9E1Isaf>
zHz^Jp)aa|wXmNJdL6v8V>^Tp!a}UCOTx)8}?uMBgI-16oSC+5fH&3q)pp?CYY$>dI
zd}_SKe6NKHivm@Ji3Tw^N(ys?yexFP(a(@s+Wy{y{pHmzW#N;;@>ZizPMH|lg7wx%
z>M5)#Ij;ts>(xb@3vvn7_{1`wyILco3Tt%Exl(d1K<dcFOje!QL*{%ltSUX^WiH_m
zU-eSMJCAUgCI3&_s;MEY;#;-8b%q^qx52yY;nM5)K;sqEQO1v~H##>R?4rGGvtU$b
z?Ft=3n{hxhs;h-JNK&pKv_vBb_0=|D`Hx-MLVi{Zs%5YC`ebJ(lJ#}@<=s-5^j;9R
ztxbMB&C0a1k5Uf}Yog98A@F<Tr~Rm;1}fiDX9ij8HA6TvP3)CdQyvGI7Q2DCrHfQw
z=0&8Ei@YS?fCi(MJH=ZGa720tKYLC--(KqeGp?ERhH1f<mJ{hu75J9O+tSfq?03B8
z_>-?og+%zGi%XAnn=YuG6438Vf<rnkeu37Xen3c0P_1YK@iiQP;dWYGB_Vs3wPis=
zA2${DG?yuXKMU&1%(*yPvLe#2bkclu2!L6kG+H}P(&1=~^p{&%2Yt$uE|#V;+pYoD
zS?R<@sa;N5R}U-)T=mjjY(?y@Cy(kt^G(g0w~9@rR$d&isLbXj&}8jwx&#ZK`!Af*
zb7>6_7gldB1u6I~f2AM}`&N|I^Gmu%ChQ%#2U*5fk&~|4njl>w%^|fAEv8KVfkEEi
zR14{bRH>!rN?Nb+v=W*gxJ+mRPnXmtziDWar|@bWr16i`ngAQf>(~urIJ$^5r{*5J
zw91*UzvMq*j%B}HbB#~wr$%>%iFxt66gP=!k!Wo^2GJ={WvVY!m92U7WIT46wbEKF
z`-o^F09IG(9cXIM!Rt%A=H4MJ2(^YWddNV<et{gVoxBF2)tMCqpyKtHl`Bn}te-*T
zbEZIhDDVdnew<SRM}bldn~q>25uVXybhEz~SHiJMMgnstcyX~BB5JmQ)K{pj^Yc1n
z>uxH?Yidr&t*3SsvAEv<*iDrZ`Qx#`hbeap0fP%C4FV%f+`H@|uM>(XXu13K4SbE7
zTA+~&{&aJFKUre5SzoGQ4rx+i0cYA?AEJnQK_J%d-`~6A#|O7?+Cb3KXXue8|9&AQ
zVYJ~KAIUal!LaPGHr40OK<6BIWfpytCnleR?rZ0cZsTGMzLvxZa^yLXcj}TxTkr4_
zFw40fb(`Ze&ydU0y4ji6BMaT5e_vhJzsg5a1MQ%K3d%$7Zn+ipm1#+z{kYxi3Frl+
zOz&#DTh*y&sv*bQ`bj7lZdE!iFBF&8m-m~2huxn?ICUsf(`_l!%+h+Xp%e~!zHirD
z3e-#9^$`BGGFs*1kZ7zN{8LMafz1!;-(QlX>G2vQorj<K2ye6^eKc$WKQIgWPIIT!
zpLF9!?_l$22gM0Vrwdi8mmH`><8crp3%L>#Sg(Dfkh50|Z-eWO^Sq-C<fGHd83;Qn
z^cT(ScAbcsq~8l~26V6D<EX(Xv6<M<3YAk;v6!CiUda@t1c(~KmB)zbu(B9^tkzis
zRW`aTqAeO$Mnm=ceZ`q@rFS92N;C2}kAfH*hftb&2AV-LO8Yo?mD_Hh8lUk`>8T7g
zX=gz=@qV}?F<|yM_26jbH|(drhwO-#5^l!`IHzigpMHl$$+8iaUbwPqOwZ?XvK6U!
zf98!wOt_rey;V>B>c~SEUSB-tJBd!q?YND&(B?VWGZ72Tss6J`zCQMgoLNpO*dFVo
z$N_;vOo2{7U1#;h{~l&pKL4=I*p&M@v?L(UYzIbiGP8_8Z@Gr}#ccasl69cs{b8Bg
z;^pxx*NUfuM+i?o))3gk6EIX<5T<dOE;5opASP&0C7oRw5c0^Y*J52-QhcnSJYd^q
zWa>L<hrcoOkHE;-Z>|EFc>hWN^-m~--c?-R7J1fta)mTa*;cHRI4?wRRy^~9_I1C}
zdYAGwgvcv2k}Lpm<4FEHM{jZeJ3^k~O-H8G7;q6aT{J|xexMisl!0>;t_>H?{o>By
zmFI0cw5rR)=@F18C7pXId3Aew2l(`(^*?fZuUl-c|8z_7uWtPnVSqiBp|%Uuhwz1u
zw}XwlgSVE4qbIql=PO&c02VhF53qy%{RZ{9m#!=VfGgX*bfG{?o(v{b3QR$VNU$4R
z$|di)$eE~DI?@+tKrx@p*yqZU1Z;DGmuks;eQKvC`gCJ>80l-i`wZq>*4`Kv_2d%c
z8Z(aONw_|qK9eGbW$GmciR>9xxc?r-^NV~aI@zfijGpJ$v;c394llQK1j(bNo!lnc
z1gI)v#|h<&h|H%?m-!nfmC2=jN8)vHJT(0;RpuSYl5D<w8&80S?(wM`?$s&W%$uBe
zkf^t@qDQl=;K{5(`71$bf#>wH3LwEkzT!Y)M`SW%NMJ^wm%Ii{N+C6ZY(f#QK_<hV
zlCACapM&e*0m3c35O`M_*mLRpcH$tiJ!dcag@2g%^K?w?e9MoB@v%4Ir%()-6?N;X
zNPY?a$k_+oRR7u?RaLR14Pu78<_C`K0KQ@xx@joJ4eKK*1zssf>hcx-fOs^a%eyE~
zXL?$e;+;zD368VOIIyRT90Ik;;<Ekm+Y#azUBn!1Dx?a;Y3q`tp4gXm#6tu3d0C$n
z5)wk9lkIp&$qv28u$IqTp@F}Mb)LOEI<3B9Q2(_y+en0poHfUXf7ru{S2F~s)+a7M
z#@wlj2og##hd6y>UpD_mPO*Mc<p`hEl9Mp$sd%Q!Ce?~#?^r8IDOyJ$`n$|je{JE}
zl$JM=fZ_6}vvMkXuJdAwc5D&uvBhag>Y<Dk6;+MD_(nHDI8tCNK;*2)dAu*fYz=tb
z`I9u~o#z7ni8rlqUMs>>a|Whig9wh-Qh+9&T8f()(^dt~$zen96ji+Le(8Ks*!d&B
zbP2g$kgZWW8;gMy%DZm09}{81W7cmSL?C{al)GKBzXNpFZsQVor(#GJ!*N?5l9Sno
zHS^$XPvjBxmVA%W@GC_qu3?+06gHQDA%?L>A*X0|iytjvlBOCR6N6>%liMp=3^{vj
z%liv1G-ZUzh+15Y*QPvP)nq?ZsgVKGPPTCxiGybTY+r{^=wyqRPBahj=<Cl&(F}CH
z)T$Vi{!{$kwJT_VN5Tf86D9X4a)bLzh%;9z4s4t>U|dj@Dd>4~cuVx;BSGE%tC6}8
zB-Hl!DJ`MF_(DxFsa_)UjL<6QCDo5!vw-=!QYd-z{$vD{(fhbwIM4<{bM!-3-|OI$
zV6BRRWVw^B?HSj6#|;IuJ+p+LmhJP~asAS3W6eLx_yo^E)|JvE%a><jR6zQNFt3ft
z#vw2AxeNp!o1qYIJ7eEmlVsi4-Ju^^qOQera+)AFuzf{;`zK@3-YB1rnd#)4uzk59
z|LDW5XyeD?hYom{9-zB7{k`A{+-^?Huq@Shk}_4&#%HV<um{~#qCO+`(p6X&FglIT
zYtJh99pecT`l&z0o?5%{I_lSD<j1WstW-Xhc)f-_%j~>`$=={Pl>eFLsqFq@PRwMB
zS=yNTSp;5gIHiG}H`lf(U-4J^!LVP=gE6Hxh#=Z<NL&aX#d3iFGj(&-*;_`z1fAPk
z^YYLEPwUlvBJr)NniuOtk5+9vbrN>*&)2IN25$DgO*ft{HQ2Nr)aA`0>K;d&%$+TW
z9TH7i<<5gTL4iju-0n+)2`l*#);zGXa?YLIo*Q&Wj0sV~?nVsgw*u(9J1HIZ1ucXE
zW;6pcKm<tTr?aqYw}k%Otu_(Os-K7Od;&AN9S-oKZ*V=|*R)(d07{s5edY5fh{%F|
z!Q>bq+g2~NyCD;tpQ#X}R=72J9g)DPPr}Xnvz_D-jRLtu(tELR9c_v|BU^)BD?9GQ
zmWxwkfqHe}EX&@oyj>Fql(^YB#a`muZP^WuY-)xynmiBMoaO;e6BJ&2Kq@*po4k{`
z8_4CgnLC<2>BJ*Py%<~)8BtcWgWLD|eLhPU_jihN{v98mOFVW=&F^B*)#BP%5XnyQ
zRP4m@$cxhAqM9`ysTpeG@aW@aK1c#a)~RJqoP!~{)4XKw<9q&j-gy{@y{T&WyRP@g
zu`>yFopR=0!3y%<);b3*NtxkrIUmGZtKrBO4Y9*LJINy}5LWJ;C(vxG!`&giVyP^t
zkt%JfF~PL(hG!6HNXg2}#_!+!eI_<Xt7W`6ej1a<-!gL*E3D5pQazjUA7K3;=Ddls
z`#G$IMTq8hL5K!xd%Zj&z>Gkc6H|-_8Z66;Q1U}2Z@-lEO}`@1)z#KcjP-_U4u5jf
zaCDeU{>ag#33KAH$D_Mc{N7$xIo%r5p&Nq*C1aFHi(ctHPTyKAALkhP<JP%`!*248
zb#>6fdv4JEs1ooVjVN)FY_>j+N`3bO^X1iy)fVGmeO#v4q>%`u{95t;MQm7WzZ>^G
zv$4I}R&GNs2ck|gGbtQ*CEnQ*z0!Hdu(r!)mQ2iDQHd#CeG$l5Nn83w#q^<W_mJD_
z+fyGpUw2Pp<20%EFLHGCQ~0Jv&fi$Ym_{c3^NPLPfBb58Ug!nWROoD;6vXh){Ajvq
zmGqnqSt6q_JLzFx!<cPG%o3H-aBHGMUHg+Q_*<QNDkNeFA15L&CV#bMx#=con=pB#
zASCB6(0gIP++!8&Q6THA;LanM2x0n1CPsZ>9r9cwLCcfFSBRdUnDB9;KBtlZ42@WQ
z3(g#@=_KsKrU@MT)o*;ADURV+FKM@aOv|}&V?iD4T)>hGN--4axehaUa;a)ekR6lk
zz7l$2J8qI}qzeP1vvySpx_Rxl=;NOlZv_?4BnLVJj5G^}1gyvNt8mh9kR*&S<qn&O
zQwDFK;oS9{>~68F^SlzOx=$AoKKDjyHpqSclEa1hb1J~4s$CKL$WS(Kp(Ei)d+Akv
z+?(RZFn8pKh;h9oEF$z~HVbNoki?SOIU{}xzZR3u<F0I&P9U$pxA;{9UsX2-C&W$b
zAZ(rtDeS=lB+{73ELC6sZpne^x6k@|(NQ+>68teNcfXQ5!Z=uqCt^N6joe$DTIA%(
ztO}ePpUpTed2O1zGlf|pI!74Xf2z!gZ{esoY~^omdz&OfQ2dk%e<D^v`h<IXC0Yr>
z$Z$@RdZ`E<r+vFXT3CHoQc)QBjZFSk++`t&XMIL_O^(%>nWjY(tAYmhQ309ow(I!K
zkr{DRtYk5MdzpZJuZ!&V0j(b8Vs9?Zftsqj42pxt=!@k;C>ANiRO_?hbkaB$>e0ji
zP5;L@L*u+59@;_61rsYXrKDoKnE8p(j|c*%dQf+%+V=D;;oJNo7LHqd<Q6IM<K_*&
zs@z&Pvqmg%vb03z&!m?H%Wh$)iR#9RK{_sT>XARcBZ9Dp=chJjNA+vLlqrM+Jg+M(
z1g?K%2}-`Wvc<Wwr*chg4E|~Z6o<IY$4sTDji6aPo%MaFCx952=Le(0QDg4nPR8!z
zZ)_lhl`3zd<Kj}=OI~xPk()HGBkI#uzgnQZS?jv7>*-^qD9_x@t)vNI2^B@3C&Wb{
zGy@>-0cV35kv;?9&(-RDDH2>6R|0<s!e=Hqk0OOJ4=(d-)^2-q`XrWFTE6Qj_$ht$
z!aG<fMxZZmhvZPogSghjRr_qPt#+m0=ddnwbNWqd5RkrawHu-muoH4JO20DF$;-p@
z5r^%tWHtTrHsS;_4h!{oHLeB_mAT+@EZACc^F749{jf|VsLoE#ybQ6}6CD`kL`KGr
z^l(lhPOR~1y4d>wOS$;IWM(Q7)uG&9W0bP?lzU|oe-&NFCsN36l2-@!H{>{yb4ZMP
z&QESmq5RDPwk#A427myK=beVmiIa@!C27qU0-U+Rv9!L3UajUcMnKc)jI@GKspfOM
zn;oMkg;NLtg<jJNsB7^jCb6Bl<I0Tlyq+c-;wj&_NY>m6wP?b3mef<D4nc178DJ;4
zHIEt7X&FUtr%FT27_^#d=?%(&V6JNI4x5soPr|hSviwQ<LiaxRC#C2<x5Xv<4dQy4
zkN|t*2#JXcxk}{Dc<c8|U-57YC5snGPzM+4GCyo06tHMh-*rF9d;Rmol>W{4|2(mK
vFYDiY|G#-+_c!-Hfd3CK>%M@@KSKQfLHO6|yZ=!9uY~t>lx%+DZvp=W4Y_{3

literal 0
HcmV?d00001

diff --git a/Inputs/EventGenerator/ThetaDist_proton.root b/Inputs/EventGenerator/ThetaDist_proton.root
new file mode 100644
index 0000000000000000000000000000000000000000..897d7d38dd876a4abe2a48ab1f0e8bc52d4eb587
GIT binary patch
literal 6833
zcmbVxWmFtp5ake<5ZpC`2Or!9cbDKgxD9SY7=pVy!2<*cp5PGN-GW1apn*Ve4ZeKY
zANy<NY@bu7Uw6N%dw=z<SJ&Oe1r7lGIRgLyRsaBb&2zXuuO*)Y={eL;{^#PACjfvK
z4nX={hOY!*$J<FKIqe~&>o|V?`u}uA1O6+dr9<+%2_W#V@-rF$0MwSVbF?PXhFQbS
zW$ZlQrmpTTa2IEG_y1V^Z(aj{|5f={2>^Je1yFmAf7Jj0vb+BqMbwD@m|$r9%NP2`
z<mkVA6m11AIW!MfYfHGhrxQ0W{69)6)IN5uc2=UITz-ZsuVhdG8W;cs(zuSFtM(kx
zQ;j0^8yGS-jrwm9rN8K;DY%3m5=m;0SAD`0iX<${wfU8#XsMi+=zXL0P#EZ>8iS->
z*`ZOz#YpG}na!qo55#x#+ZU;P!WYK)O&Z#x8C4AVB<g~Oz9~Mjn3?<0WAZrnZuwex
z_hWWt{35gK@5p;}QMehT<urorMdZ<*8;RJS@i&;iVNTPx??#*kMkbr=nIDsqZZWkA
zH(W*RXUM!FLOR`ZH8ANWtXZl(E7n!e-k-gnWW$l0O%A4~r*sv<OHyg26jOPRSQyTT
zr5XnbcC^>5_y@+_yon+cgJ+sRJqUgzIkf*eXzfq)f6I0ENdU|PCs-xBJSH}}635$0
zUruZAP`Z`<x)s0ps6K6*ovWIYaVL2|RUR#G*qKZA^~$EL5Ke)o5!o_DpCITcA~X2w
zsd{{fjfpMii0$OZmRmh-vK|ID#bV$q6Z*B-Y?8iQ*0mQu1MY}VoW9$dvbYD-UwRqr
zUCZ6WgG;pXUYl+=*ETRXKz-DDgVGIomVdlvMY3w-6KGmI<;Gc#e`>3itCK{0GyhVx
zi+9=AQ#NtHlwh2T-M7^Ab@V+~3I}>(B3JW`F}o1|%J?>ZLi+2y`J2P=W)6glI(l1}
zaz=<01VL22d|=yz(s+V#yMq&y(4x(=xwdFe@J(w)Qywz&q5a1BfQ63REBE4}quQog
z)$|kYDYZNYsBI^VC}OoZ!wh`F4rac9k!2NsW=Npm;b-%MZj`3FjdeCSEe;48Z5kB^
z@Z`J(I*Pu_8{IkS?)my6H#o#Ar5A|V_9d{jgpe5g-SLJC=S@j5*B{hGw!>s&+_Z75
zQTwLhtIVU1Uj#9fn;p{Bl$F&wpHLI!7HO%XK-oEBjA#BPkkiXAwK%Ohixo>CEDcXJ
zMM8b<)e<YWa&`e3yXhxHoNpcF8JyYbZj<^FK{eZCQqtXSoGJniqcq&3t)W$zk?)_X
z#Bk%i{z?5I?M?ys8g>}zYbfv~c7#K3x_BH*`(_?e?ypmr^>!1c=M!lE#9%j>p6@ql
zuidmN8)R4V+pJzT!o~rPh?c$dd(iIu5*LptDl0En^gZcLsPxhB#at7i)<ZC)e1ZJ;
z$t*5Ki@IV&IO^~9M++g+Wr->g;`}auDvIA1GeL*SrC${bAE;1Z0sG89(n0b*yT4{e
zQQ<$%YqGB}_oxFj15L=O%lK)7-gOe>ui`b=pIe9>{`Eth+itP*@sITP8J3-UB>bU^
zSMWuip!DpN_y;#wqO+_HW@(A4v(0?Qmymz92s1opEA)J7{2d1RIKqHxcgWLUe=ABQ
zy5@l<D!`eg%WqKcn<*F5&bNJ-Cv`;oy$waSU#D3+!%RamKEdC6QMMB{#q1&q`!P*C
zqKGkFBrjnRpPkdYS~QtZcsdFuzfOu7m=B(w=2LP;*O=D%esT>>%^;BviTGys-A(`t
zf?wp%3Y~3IDES%1-om<Fkdcy4$tD)vRy;H%@J3kPu(mKD?tY;wzwNii*~|gC_YQ_?
zBf4rXbUS$v16q4SVa)%~pX~Ar<=51&EZ;9`!?!9!5_QT)J|pxGe$-wMf=o^6cJ9lR
z^i5XH+XHg+qo(uQEi+1KPV>L^Y3cK_iBd=P_T8s&=}9hXw{sMfi7|tH!KUdv69(KR
zO6jwz!w2mxIp0xk*0JelDpFM*m4=i2m9dTO72}$bccmFr0^4j0lD#m7G7$8w$2!wB
zJ@^YoG#x8NeY=IJ5E9MJF)n(w6nsYR4+yVM_^Riru=_&%s(?}%zY_c^P~hU9kW$gA
zzWAmOrH>$snBdY#b&CfQr8Z7N3s1(ra(~=WMo9#}#25wObXzdepa$1lJbY;5IC~D)
zK#!XH5u58uCsI3@)EJmF<b6mGiIgYR=h1P^5^3^>+QLQ>GI9R&GoY)jk!*)ZeC(Xg
zEcYYUd2fO_W!}s4hbhjrT=cpc;Le{CGs?1XgLLefC#ErDJ#*(Lv@s{*u{cFeR(tb~
zp}X1!&b%kgV)9pYyuh8R0`rM%q?tK%k768~+}As^y_~~nGi1Sn8z^HDg6U%rZ;>(c
zY%_-1qjmla<lCOus?enp%pGI@j43|ORtHWxWX`W%UEk2JEJ#V{W9x`x=c!`}$YZx|
zIJ4e%|0UR6!|X;B;>$_1$KJw7J)Jn;op#IvEm1uxGVA(rEbD-UIg=(A-Au}XHN}_%
z_B&rTAEK#aONe5<`8cBu&BiQ_CI!WfLr&`evCOep0+_}TNyA*CBy#8;@;qvoH6F82
z&fBK0NZaEQvpy~X{#SHLA$Vw{u%J5Wkn=X+F8vh`pR2q0L?;U*flF-U+Z1ruAr^-f
za=q_*gxxdj2VwQ@%)GV+kswc$koiEwy15c8DF@L&RMqRCUIf#YWUS?qC1I_w_he_w
zl02C}^|1|yttpz-TZ;r#hFt3^%U#1d-@l!;PcLri7NGxniDwxWbowX56938Ye;o$E
z3aG8}+y-E3!QHLRovhszoo!t3Wn4Th9KER>?VKH~t)6cl@Mkd=#RE7AK8x|3k+CCa
z!BY(6u9&2pQVbsEH#Q<*Ok<8tQxH<M7kn2^*l+k1JhCIiU_Cm2aD*v1?;f6&Nu)R!
z)uiwpMFE4?DTu$6>h^qkFke;!neB&+-An`_>@^|qemehmNKk>o_kwXP^k*6VnfcLx
zZP?iG-G=hm@z3Tp(cN5+Bf(u0!;4UP;64pu$XljC@V@u1U+T0wNa!k@z~p2#@Z<^y
z_MKKJc>he2W8&ena;e5=s)wbzyxu0!(#vAsR&9SjcFaSH$0N`EIaAj=rgq=j$bG>K
zZ(;YphYqtb=P|Pp=TU^QkuQ~Cb7g&&rF%MKBxz|%Hq<UBTO#q$n@iM=bNwsXwfmJb
z?=}f4<_V+_wMkRU-+Dt`HyAG-%E0h$YoS5HtUVnH@hj*gS`-X>&=2yC4S7rfpf9sP
zK70P=E34J|eeimQL5_dZc0Zr3K63xPgt1v*WUei_49AxJuq%aX_Mj_0=*>8Yvg%Ib
zp~;P+w>m$Xk^+k5Ume+zW|?-<QD4U_twV=Q1pTz7f-WNyaun0yRd@u||MoUTGksU&
z&ym!)7^pREqlL=bp+k#LRJ<{gpWLc)Yd<mQ+#y&vsh{oJxe17Z&8H@6qfBR<sd4MQ
zl{3)%L!_J=9D8VN#g575C$0^pat(D)gp!Aw@6<WP?U$f)lJB)aw<PL$&4LRfs=m}K
zp|4({8Q>uxPQ;$!Q0H@~I!4bY@}n{IwM|q?95Mpva+*{<W$0KjHdor#$E-q5+78)5
zFVGf!6?qkFEL+DRD5N={?bEk_WA+u(3wyrzPpiXg>o#l!dAy-_7bdLNVJ@ZXyoA?P
z$w4LQGJ%xi9(l`zM+q(8c$^lykTE4AB)Cz?^4SN2yf5OyH>dPM<*z^YNvwzKOXbz|
z{^8VSq`DsnNOUc;B{Hp!6<ZtmEF)*S0A)!ruhFnl7Wu%hBGdcTUgs}df0dO%Qlnym
zzZ}l4s1ujQ*CE)NmycE#&S#jO^8#oc=Mjb}@*%W5CyB$pJgr?I$KpkZ5A9ge9cg^9
z5yk=OmgJht5aPS(2x29r7Hf<xmHoYj%Fi|5<H=@K?MGr3OQzr1U!-Z$NW~osd)S}?
zWhA!NcH-HMiZXPct#@Azd<u`6!nj90(e}EJwoG@h5W4X3aDD{Yuj%^h`kdV@r4NG@
zRYukX2$=&|i|(Q57@CG=%yy`zey&;3pBW@^_*_(5LH(bMqsF<Rz&g#kLZ;>G?Z?jn
zciqdwqNW5{?<KV+kNHzu?u);}1Ja7Y;61sn%Un<2o?G`Xo{9K=2GU_h^YWX`2&`X{
zi&lUWYrM@AW;~QK{I|a{aJ>34p*$v&N=W{By;Eha_^>zwx?b8^;K8W~>7VGm1%VN@
zAy>8Wr9c<LZp>=Y?>(prHgmZ=E*v(9wB{HrdJhb4G1X%r)|+@jjkmtyKX||esIQfR
z;v3_!6eZ)l)RB*H{fx-l1P;Z>L10l!`9B3I0oT7U)_!rQpE_X8Z5w@2z50#L&azS)
zX+p0M!T%<q=uop@H-aSieRHFw1%4+sPVqqbS~1Ka@iAG!0m*OA``f1?PpD7;1Sq+C
zr!+YRo^IV*yFnM@0~h4+Yj*~${BQ|TSN{}_-sJj4sk3UxHCHM74JLjTD9_umdzo`#
zYT1mfUf^n~VY2gff@*WL{!KeowN))qTSaga7$F!Ute#A2?xxKS7PpXHX1V;Z`SyuN
zwp;bZEb>US$gyQ?%rRQ_fM!+~Lva+GW2$V^w2&pW`7g;As6~_odwbY3>}s-!i$YbN
zno$w&P9|8g!7Q?E2p;miURAfb99v-w%|DB~KS(?(VPa#5h=S|$vF3Yl#?64(i?o*9
z!UNW<cd>v%Xu`+FoSC0{P6{s*zd95}<F>N&so6vky?GZEf1>Cu?bB=LRmlbdVQe@V
ztm~q%+qQF4j%MTrUCYfls~Xh{G~7fFS4a9vhNx}e01G6Mk3L6hRmhfPy{tUAnd%$s
z8D`e<=VG)WQOhbZ5vhaEpwZ+QjhR2fO^s~uWShI}6lKd1PuygX3V?vE!0jrL&vMLw
zu1ZIS4r6)fPaE*hCoP!^oaaAZ*fdD2Hj)L4{J{@;9cVz5F^BVdC%OxIBf)7p8$z0d
zl$9zZ_JN|0#_K0`I#j@YQ=Y-5xTJ&;*&D)aZ_p@RSL+;+k(^(|XUn#Z$9oZ^e3u-q
zDNO$Ol8ks(*6Q~#_@25~IrW@5jNN?VBT@uKqW-B4lHPuDj4pQbu3u>fL7cI@TK7H{
zjl>)(`=N5c2V15=TEiKIs-z%BYgG3jK70ZVKg&M~%#nv#&CoTo=i?W`;SU#2Jw~a5
zv@r>jJmhu{ERRY!dvFTUg%`%FV-_ML^dV|Jqm^5!b#U@`cQ9<-BuMi#MO=m*u#$zh
zgN0Wt?el43ft%0O&qmnlW`oCU@i}5rl0;zt_13NZ&w2+>d2OpY;yy0Mt`_nEHEK~)
z$?`v_;VbDCve%lp&R|-pS3PvT4hwApj@U6OdM5AWHy4%{j|KCzvZ|$DPB?H^i|5&1
zG<PM0yhlK9c-_o&+hKJ!$6ZuuRA!;hLVUVeC2j{>(aLlffh|k<h5TQ|$RJf1C$_a+
z5A3P#4p$qqAA(zIYwPE-ofFMN-`HmINio(AAs<aMoFgjtcwU|b+vW7u;KbRzyEx>(
z?*LJ=OWeMY35deXX$rmQcYJ2r-@-BF0>7;(xg=qLOtLp2J`7LpgvMOt+uiBZQ3}!s
z+vE*)&pRwkb4cD5g#i>CE?LMBnvagokL)eyo(rB**wNva?PO(`QsSOFyuPf1D8r35
zh!fc$-6VbX4Wayv0n>?#jq)gc@Cw)HG%MVUd@#2UyOhiv|MF$+Y&D6KBJ9n(LpJV{
zdOGQ(>{$84fXnnG4^gV{#OSw7lIlEMQc@EG7AkG%JGhLPA&Ih9re5}#>$a1*@}n@m
z^>;MbKsBDOT@GI@-^Vl&j8>_OBP#J5vU>4IB~FTiOA_2$33qkU)KSj7%7mna$<(+r
zJF2Ru59@yOXZ3GRW%I(ibTH0HhPq~~$p)a<Eo<$4@h1v173=c*W|fdw`Ruz}0-%?5
z?|DK`A^*qytRfaDM`+87R}gTK54!z~r{p_1t21&k;-qMM4=$ke((>SbcJ`3Q8?zw_
z3QXlnSK^9(RYP+%ip+>=9zwY&y>19a54()OyQcEJkz!e+YT`^wdCiWR;h7Yy>MI>q
zQEltE6xEoW1E22pS2n9Dt_NmJ0=&`k;bFM4GVHtlm!l)&2Q@?8A_ifkeQXW_m1QsE
zjYkz)^iwNCe2qgkE%4fD@;g&M3WnFLG;V-P(cUwhOj5%Y^j{^s+7PaVZ8E#bu|`w8
zM(P;7VQ0Kxsl<(Ni<IncaN1)<l^mzF|7{1)#uN#gPKR<Drgng97x2Yf%mydWB;}ov
zIEqzl?<xf!`Z6<@c;1SbPmJaI-Ot#TEaiCdbhls$T_s3Ns{P6vJF^5x8;5&v_Xu(o
z*E1vpcI)}DZmVZ0XUqKRQIIy+^{<&|s>+TEn&(GnC_6S)r+V3PoM(gu$Los4evHbp
zy=h836wc`^_nL`WkCG|L3TO3WoLagldiv%)B9j5xo7^18GpEYe(&+}%`VF;Qmel!?
zx@WRmRfe39PUl5Jx~0{+N2GGcQ}Cm3Uw=C?oL--qLoA(+WCj+Hb0C3Kk!B<<QzsHZ
zarJ{O843fmzKa?pjG5CofY}Md-Gi=bT1DW6SB1J5V-y8@zy6SieL*cTukBLx?8w>p
z(T`WHcBp-l_^N#cQZLo3qV0DJM0F2GE@R?5iL0o}hGkKZImWTTlQr6KwfPRsAziGy
zmxLGzyaq#OCEUrz8;8PH(&!0G2#EO>Qq7gIx0P-%f6`wc&8W#wp6r*0?k}1{zZJ)R
ze$3V0$9W+~O!l;ln*zH}$mrXVVk=JgG;MuX-0SZlBmOvZ(#a*p89jh=0*oy<ZPzVa
zOhAVvcl=PlnSJg_-r8#lrt>L;OuG&*#g~u<xu|7(Fsx-8nY1~}y^@ngY&qt?X~y>;
zve)|gX#Q)Vwu&HPx>_P_KG1rzWLE|`-u8}$YqNj^XD<-YG~2w|UY;GGk!6~^Bc9~+
zh#QeV)=}V`ciaj<NWY!<1eirxcO(YT6q8gm-)xL7j`iN)^9N;Q_2bi4@(MVC(;=^P
z4L-JJrb)*=w&d#p{rF={;F*-bfB@mY3U5k;wM~VQzvi0|%;wCIyK>NN?F?Q6jXDfh
zWAHYz$ux+_R+iZyMc0lsT4whS(so^~1xVwcAutIobw!M<Ae@dk9DSrmm;exf2Or%G
zbJGPhXwk^aNVW)1yqF9&#cpg^sal>hQ?_m|ZMYvUBnbH6mmVo(`7KUu-bQ~)?xy~2
zDnB3WzM3K#sYswf-zuF$ES)IEpU+Qq=#7nYm<dQ&4A8$788Y5?vxst9IS<&E122#B
zuqpP(1SOC%QtD|L1?y?K^cMr_w2sS-we4;8DApgy9CA$w2P`)GM1LIn?$kF7*%Cn5
z+&HRADEc1e3}_X+C0z_zN>!D|`%sFYigb22!H~-bJ?WeFr=KhvQEPYT&n*%nGxPJY
zfyNAKEp(4{9CbMa@@@@x!`JJkE|uB9JN#wTc`_ez^B6wyXB&3593M05ft775&|ob>
zPKdx~UxX?MeOVZL{${a(7x@YsR}biRNhbr%BqFyJmm+_DP@@@m52`cb4xnLWRM<7|
z(sn<CfT_mvGJ0#<h&r+8TKAkPw&5j)%g#aIbbRnP5?_wo9EAOk3tY^wMyj%?l$o<3
z{c7-yq~j5!#sJH=Bs<Vii=bUCJPgaE0My)2WbS2v+5R_X%U!LK(5PRYjwgAIrtu^G
zJ$;-3$(RJE#~t#K`y`OYW{l*ZTrWYajB-pp*UA*r#Y75PYD^?t324W7@K&9~;7~-(
zO~+n2atvyf-rS*o{i}uRfm0~?oCB?oJ_zRgB;_^7UMPvQ{t+87=XO_})ZY06_ig@l
z#U>i>b;Cv1skM*X53BB-OA~2h2vR%i7Y5CPIQn5Zt|Y4Nfzt__fLs|d8)2Sm7FLlq
z_L*2iv!R>(S{_SEX`WV2u~xvS#<lswLSN`XyKD$$^9j#Ls?}8xsD4~#*Q%MMZgL2{
zVic2GlN;$%8GHTbk9t$V)j~YAw;QyTV=Ux0NlW*MunZku5K0iLMHvmcV&@_|<YSB_
z-&txbeFp;9RD#WJ;iK08{akSklM}y3pt}2#Zl5#A4rj4?2%*pb13BoxuTH+DRHJNG
z9MpK_F6W#tR96}2$Yi^may5OAuv{p(IyO9{=iUv&Zdu`G6n+T|wpHU@7L?M_s*yKT
zDlOd*r;;Ni@*TZn39|f+oPPEB{D7wRb-oNke=AV_msZPK_aQ9$-8dynt{b^5k|I0M
z0DeFY#G^<?LF0_)*GBP9kJKLRq%{uZXi*i<e4nfLfEn}tr?gWfY%;+Rs3g<f^R$uw
zd=TP@nKMfD{z{gHFMKeTAttu5*LZfp@kDSus-$f3gVV{Al<8+w*`^-vU8l-R*QT&&
z+q$o7tAEpzv9>v^a8Q(w2<KOEr(^Ka?`5?)D%x$PGE3t<xgQe2CSLom&C`RKVDHkU
zoSWWH8KaAaI;})S=u@K5xD6)`yz;b>wQASmnM^+Q!J;SYq)VQ{_Z<Ga*WI&7DzfZm
zf>cNGH?#-E-JDI%fhtcmUu6kZmG@qJ{K1nAywuLy%*Gbp>TUr%x85EmrZiDPYS@t(
zciSpMANxe0Cg}F-!Va5E=mQa^Gip}_qwX;;=K~ec>f!2n-s1VsmJ$0mh52X8JR2zg
trZE2pTjqJg^SI{!*GGBImiz~b{=Y!~IZn@`pZ^o+GbDM97hfvie*mLRkl+9S

literal 0
HcmV?d00001

diff --git a/Inputs/EventGenerator/alpha.source b/Inputs/EventGenerator/alpha.source
index 02b1c4dfe..efbad0a97 100644
--- a/Inputs/EventGenerator/alpha.source
+++ b/Inputs/EventGenerator/alpha.source
@@ -4,10 +4,10 @@
 %			   Energy are given in MeV , Position in mm				  %	
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
-  EnergyLow= 5
-  EnergyHigh= 5 
+  EnergyLow= 0
+  EnergyHigh= 585 
   HalfOpenAngleMin= 0
-  HalfOpenAngleMax= 5
+  HalfOpenAngleMax= 90
   x0= 0
   y0= 0 
   z0= 0
diff --git a/Inputs/EventGenerator/deuton.source b/Inputs/EventGenerator/deuton.source
index fd3680952..2730c05f1 100644
--- a/Inputs/EventGenerator/deuton.source
+++ b/Inputs/EventGenerator/deuton.source
@@ -5,9 +5,9 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
 	EnergyLow=  0	
-	EnergyHigh= 9
+	EnergyHigh= 250
 	HalfOpenAngleMin= 0
-	HalfOpenAngleMax= 45
+	HalfOpenAngleMax= 90
 	x0= 0	
 	y0= 0	
 	z0= 0	
diff --git a/Inputs/EventGenerator/proton.pbuu b/Inputs/EventGenerator/proton.pbuu
new file mode 100644
index 000000000..d90331838
--- /dev/null
+++ b/Inputs/EventGenerator/proton.pbuu
@@ -0,0 +1,14 @@
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%%%%%%%%%% An Isotropic Source to be used as EventGenerator %%%%%%%%%%%
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+%			   Energy are given in MeV , Position in mm   %	
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+pBUU
+	AngleHistPath= ThetaDist_proton.root spectrum2
+	EnergyHistPath= EnergyDist_proton.root spectrum1
+	x0= 0	
+	y0= 0	
+	z0= 0	
+	Particle= proton
+%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
+% Supported particle type: proton, neutron, deuton, triton, He3 , alpha 
diff --git a/Inputs/EventGenerator/proton.source b/Inputs/EventGenerator/proton.source
index fe74556d6..9bb498759 100644
--- a/Inputs/EventGenerator/proton.source
+++ b/Inputs/EventGenerator/proton.source
@@ -5,7 +5,7 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
 	EnergyLow=  0	
-	EnergyHigh= 150
+	EnergyHigh= 190
 	HalfOpenAngleMin= 0
 	HalfOpenAngleMax= 90
 	x0= 0	
diff --git a/Inputs/EventGenerator/triton.source b/Inputs/EventGenerator/triton.source
index 5a9aa2a5a..a3430584a 100644
--- a/Inputs/EventGenerator/triton.source
+++ b/Inputs/EventGenerator/triton.source
@@ -5,9 +5,9 @@
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 Isotropic
 	EnergyLow=  0	
-	EnergyHigh= 12
+	EnergyHigh= 280
 	HalfOpenAngleMin= 0
-	HalfOpenAngleMax= 45
+	HalfOpenAngleMax= 90
 	x0= 0	
 	y0= 0	
 	z0= 0	
diff --git a/NPLib/Physics/NPReaction.h b/NPLib/Physics/NPReaction.h
index 0869dae2d..61dd7a2fc 100644
--- a/NPLib/Physics/NPReaction.h
+++ b/NPLib/Physics/NPReaction.h
@@ -9,23 +9,17 @@
 
 /*****************************************************************************
  *                                                                           *
- * Original Author :  Adrien MATTA    contact address: a.matta@surrey.ac.uk  *
+ * Original Author : Pierre MORFOUACE contact: morfouac@nscl.msu.edu         *
  *                                                                           *
- * Creation Date   : March 2009                                              *
- * Last update     : January 2011                                            *
+ * Creation Date   : April 2016                                              *
+ * Last update     : April 2016                                              *
  *---------------------------------------------------------------------------*
  * Decription:                                                               *
- *  This class deal with Two Body transfert Reaction                         *
+ *  This class simulate particule with a give energy                         *
+ *  and angular distirubtion in the lab.                                     *
  *  Physical parameter (Nuclei mass) are loaded from the nubtab03.asc file   *
  *  (2003 nuclear table of isotopes mass).                                   *
  *                                                                           *
- *  KineRelativistic: Used in NPSimulation                                   *
- *  A relativistic calculation is made to compute Light and Heavy nuclei     *
- *  angle given the Theta CM angle.                                          *
- *                                                                           *
- *  ReconstructRelativistic: Used in NPAnalysis                              *
- *  A relativistic calculation is made to compute Excitation energy given the*
- *  light angle and energy in Lab frame.                                     *
  *                                                                           *
  *---------------------------------------------------------------------------*
  * Comment:                                                                  *
diff --git a/NPSimulation/Core/CMakeLists.txt b/NPSimulation/Core/CMakeLists.txt
index 0549e0593..ade632670 100644
--- a/NPSimulation/Core/CMakeLists.txt
+++ b/NPSimulation/Core/CMakeLists.txt
@@ -1,2 +1,2 @@
-add_library(NPSCore SHARED CalorimeterScorers.cc EventAction.cc EventGeneratorParticleDecay.cc ObsoleteGeneralScorers.cc PrimaryGeneratorAction.cc Target.cc Chamber.cc EventGeneratorBeam.cc EventGeneratorTwoBodyReaction.cc Particle.cc PrimaryGeneratorActionMessenger.cc NPSVDetector.cc DetectorConstruction.cc EventGeneratorGammaDecay.cc MaterialManager.cc ParticleStack.cc SiliconScorers.cc PhotoDiodeScorers.cc VEventGenerator.cc DetectorMessenger.cc EventGeneratorIsotropic.cc MyMagneticField.cc PhysicsList.cc SteppingVerbose.cc NPSDetectorFactory.cc RunAction.cc)
+add_library(NPSCore SHARED CalorimeterScorers.cc EventAction.cc EventGeneratorParticleDecay.cc ObsoleteGeneralScorers.cc PrimaryGeneratorAction.cc Target.cc Chamber.cc EventGeneratorBeam.cc EventGeneratorTwoBodyReaction.cc EventGeneratorpBUU.cc Particle.cc PrimaryGeneratorActionMessenger.cc NPSVDetector.cc DetectorConstruction.cc EventGeneratorGammaDecay.cc MaterialManager.cc ParticleStack.cc SiliconScorers.cc PhotoDiodeScorers.cc VEventGenerator.cc DetectorMessenger.cc EventGeneratorIsotropic.cc MyMagneticField.cc PhysicsList.cc SteppingVerbose.cc NPSDetectorFactory.cc RunAction.cc)
 target_link_libraries(NPSCore ${ROOT_LIBRARIES} ${Geant4_LIBRARIES} ${NPLib_LIBRARIES} -lNPInitialConditions -lNPInteractionCoordinates)
diff --git a/NPSimulation/Core/EventGeneratorIsotropic.cc b/NPSimulation/Core/EventGeneratorIsotropic.cc
index e1cdcd423..9a4f8714e 100644
--- a/NPSimulation/Core/EventGeneratorIsotropic.cc
+++ b/NPSimulation/Core/EventGeneratorIsotropic.cc
@@ -40,14 +40,16 @@ using namespace CLHEP;
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 EventGeneratorIsotropic::EventGeneratorIsotropic(){
-  m_EnergyLow    =  0  ;
-  m_EnergyHigh   =  0  ;
-  m_x0           =  0  ;
-  m_y0           =  0  ;
-  m_z0           =  0  ;
-  m_particle     =  NULL;
-  m_ParticleStack = ParticleStack::getInstance();
-  m_ExcitationEnergy = 0 ;
+    m_EnergyLow    =  0  ;
+    m_EnergyHigh   =  0  ;
+    m_x0           =  0  ;
+    m_y0           =  0  ;
+    m_z0           =  0  ;
+    m_SigmaX        = 0;
+    m_SigmaY        = 0;
+    m_particle     =  NULL;
+    m_ParticleStack = ParticleStack::getInstance();
+    m_ExcitationEnergy = 0 ;
 }
 
 
diff --git a/NPSimulation/Core/EventGeneratorpBUU.cc b/NPSimulation/Core/EventGeneratorpBUU.cc
new file mode 100644
index 000000000..eeda16409
--- /dev/null
+++ b/NPSimulation/Core/EventGeneratorpBUU.cc
@@ -0,0 +1,256 @@
+/*****************************************************************************
+ * Copyright (C) 2009-2013   this file is part of the NPTool Project         *
+ *                                                                           *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
+ * For the list of contributors see $NPTOOL/Licence/Contributors             *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre MORFOUACE  contact address: morfouac@nscl.msu.edu       *
+ *                                                                           *
+ * Creation Date  : April 2016                                             *
+ * Last update    :                                                          *
+ *---------------------------------------------------------------------------*
+ * Decription:                                                               *
+ *  This event Generator is used to simulated pBUU ion Source           *
+ *  Very usefull to figure out Geometric Efficacity of experimental Set-Up   *
+ *---------------------------------------------------------------------------*
+ * Comment:                                                                  *
+ *                                                                           *
+ *                                                                           *
+ *****************************************************************************/
+// C++
+#include<limits>
+
+// G4 headers
+#include "G4ParticleTable.hh"
+#include "G4IonTable.hh"
+// G4 headers including CLHEP headers
+// for generating random numbers
+#include "Randomize.hh"
+
+// NPS headers
+#include "EventGeneratorpBUU.hh"
+
+// NPl headers
+#include "RootOutput.h"
+#include "NPNucleus.h"
+
+using namespace CLHEP;
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorpBUU::EventGeneratorpBUU(){
+    m_x0            =  0  ;
+    m_y0            =  0  ;
+    m_z0            =  0  ;
+    m_SigmaX        = 0;
+    m_SigmaY        = 0;
+    m_particle      =  NULL;
+    m_ParticleStack = ParticleStack::getInstance();
+    m_ExcitationEnergy = 0 ;
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+EventGeneratorpBUU::~EventGeneratorpBUU(){
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorpBUU::ReadConfiguration(string Path,int){
+  ////////General Reading needs////////
+  string LineBuffer;
+  string DataBuffer;
+  
+  bool ReadingStatus = false ;
+    bool check_AngleHistPath = false;
+    bool check_EnergyHistPath = false;
+    bool check_x0 = false ;
+    bool check_y0 = false ;
+    bool check_z0 = false ;
+    bool check_particle = false ;
+    bool check_ExcitationEnergy = false ;
+  
+  ////////Reaction Setting needs///////
+  string particle   ;
+  //////////////////////////////////////////////////////////////////////////////////////////
+  ifstream ReactionFile;
+  ReactionFile.open(Path.c_str());
+  
+  if (ReactionFile.is_open()) {}
+  else {
+    return;
+  }
+  
+  while (!ReactionFile.eof()) {
+    //Pick-up next line
+    getline(ReactionFile, LineBuffer);
+    
+      if (LineBuffer.compare(0, 4, "pBUU") == 0) {
+          G4cout << "///////////////////////////////////////////////////" << G4endl ;
+          G4cout << "pBUU Source Found" << G4endl ;
+          ReadingStatus = true;
+      }
+      
+    
+    
+    while (ReadingStatus)
+      {
+      ReactionFile >> DataBuffer;
+      //Search for comment Symbol %
+      if (DataBuffer.compare(0, 1, "%") == 0) {   ReactionFile.ignore ( std::numeric_limits<std::streamsize>::max(), '\n' );}
+      
+ 
+      else if  (DataBuffer== "AngleHistPath=") {
+          check_AngleHistPath = true ;
+          TString FileName,HistName;
+          ReactionFile >> FileName >> HistName;
+          G4cout << "Reading Angle Distribution file: " << FileName << G4endl;
+
+          string GlobalPath = getenv("NPTOOL");
+          TString StandardPath = GlobalPath + "/Inputs/EventGenerator/" + FileName;
+          
+          TFile *f1 = new TFile(StandardPath);
+          fAngleHist = (TH1F*) f1->FindObjectAny(HistName);
+          if(!fAngleHist){
+              G4cout << "Error: Histogramm " << HistName << " not found in file " << FileName << G4endl;
+              exit(1);
+          }
+      }
+          
+      else if  (DataBuffer== "EnergyHistPath=") {
+          check_EnergyHistPath = true ;
+          TString FileName,HistName;
+          ReactionFile >> FileName >> HistName;
+          G4cout << "Reading Energy Distribution file: " << FileName << G4endl;
+          
+          string GlobalPath = getenv("NPTOOL");
+          TString StandardPath = GlobalPath + "/Inputs/EventGenerator/" + FileName;
+
+          TFile *f2 = new TFile(StandardPath);
+          fEnergyHist = (TH1F*) f2->FindObjectAny(HistName);
+          if(!fEnergyHist){
+              G4cout << "Error: Histogramm " << HistName << " not found in file " << FileName << G4endl;
+              exit(1);
+          }
+      }
+
+      
+      else if (DataBuffer == "x0=") {
+        check_x0 = true ;
+        ReactionFile >> DataBuffer;
+        m_x0 = atof(DataBuffer.c_str()) * mm;
+        G4cout << "x0 " << m_x0 << " mm" << G4endl;
+      }
+      
+      else if (DataBuffer == "y0=") {
+        check_y0 = true ;
+        ReactionFile >> DataBuffer;
+        m_y0 = atof(DataBuffer.c_str()) * mm;
+        G4cout << "y0 " << m_y0 << " mm" << G4endl;
+      }
+      
+      else if (DataBuffer == "z0=" ) {
+        check_z0 = true ;
+        ReactionFile >> DataBuffer;
+        m_z0 = atof(DataBuffer.c_str()) * mm;
+        G4cout << "z0 " << m_z0 << " mm" << G4endl;
+      }
+      
+      else if (DataBuffer == "SigmaX=" ) {
+        ReactionFile >> DataBuffer;
+        m_SigmaX = atof(DataBuffer.c_str()) * mm;
+        G4cout << "SigmaX " << m_SigmaX << " mm" << G4endl;
+      }
+  
+      else if (DataBuffer == "SigmaY=" ) {
+        ReactionFile >> DataBuffer;
+        m_SigmaY = atof(DataBuffer.c_str()) * mm;
+        G4cout << "SigmaY " << m_SigmaY << " mm" << G4endl;
+      }
+     
+      else if (DataBuffer=="Particle=" || DataBuffer=="particle=") {
+        check_particle = true ;
+        ReactionFile >> m_particleName;
+        G4cout << "Particle : " << m_particleName << G4endl ;
+        
+        // Case of light particle
+             if(m_particleName=="proton"){ m_particleName="1H"  ; check_ExcitationEnergy = true ;}
+        else if(m_particleName=="deuton"){ m_particleName="2H"  ; check_ExcitationEnergy = true ;}
+        else if(m_particleName=="triton"){ m_particleName="3H"  ; check_ExcitationEnergy = true ;}
+        else if(m_particleName=="alpha") { m_particleName="4He" ; check_ExcitationEnergy = true ;}
+        else if(m_particleName=="gamma") { check_ExcitationEnergy = true ;}
+        else if(m_particleName=="neutron") { check_ExcitationEnergy = true ;}
+          else { check_ExcitationEnergy = true ;}
+
+      }
+      
+      else if (DataBuffer=="ExcitationEnergy=") {
+        check_ExcitationEnergy = true ;
+        ReactionFile >> DataBuffer;
+        m_ExcitationEnergy = atof(DataBuffer.c_str()) * MeV;
+        G4cout << "ExcitationEnergy : " << m_ExcitationEnergy << G4endl ;
+      }
+      
+      //   If no pBUU Token and no comment, toggle out
+      else
+        {ReadingStatus = false; G4cout << "WARNING : Wrong Token Sequence: Getting out " << G4endl ;}
+      
+      ///////////////////////////////////////////////////
+      //   If all Token found toggle out
+      if(    check_EnergyHistPath && check_AngleHistPath && check_x0 && check_y0 && check_z0 && check_particle && check_ExcitationEnergy)
+        ReadingStatus = false ;
+      
+      }
+    
+  }
+  
+  if(    !check_AngleHistPath || !check_EnergyHistPath || !check_x0 || !check_y0 || !check_z0 || !check_particle )
+    {G4cout << "ERROR : Token Sequence Incomplete, pBUU definition can not be Fonctionnal" << G4endl ; exit(1);}
+  
+
+  
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorpBUU::GenerateEvent(G4Event*){
+  
+    if(m_particle==NULL){
+        if(m_particleName=="gamma" || m_particleName=="neutron" ||  m_particleName=="opticalphoton"){
+            m_particle =  G4ParticleTable::GetParticleTable()->FindParticle(m_particleName.c_str());
+        }
+        else{
+            NPL::Nucleus* N = new NPL::Nucleus(m_particleName);
+            m_particle = G4ParticleTable::GetParticleTable()->GetIonTable()->GetIon(N->GetZ(), N->GetA(),m_ExcitationEnergy);
+            delete N;
+        }
+    }
+    
+    G4double theta = fAngleHist->GetRandom()*deg;
+    G4double particle_energy = fEnergyHist->GetRandom() / MeV;
+    G4double phi             = RandFlat::shoot() * 2 * pi;
+    
+    // Direction of particle, energy and laboratory angle
+    G4double momentum_x = sin(theta) * cos(phi)  ;
+    G4double momentum_y = sin(theta) * sin(phi)  ;
+    G4double momentum_z = cos(theta)             ;
+  
+    G4double x0 = RandGauss::shoot(m_x0,m_SigmaX);
+    G4double y0 = RandGauss::shoot(m_y0,m_SigmaY);
+
+    Particle particle(m_particle, theta,particle_energy,G4ThreeVector(momentum_x, momentum_y, momentum_z),G4ThreeVector(x0, y0, m_z0));
+  
+  
+    m_ParticleStack->AddParticleToStack(particle);
+}
+
+
+
+//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
+void EventGeneratorpBUU::InitializeRootOutput(){
+  
+}
diff --git a/NPSimulation/Core/EventGeneratorpBUU.hh b/NPSimulation/Core/EventGeneratorpBUU.hh
new file mode 100644
index 000000000..ff4b6096d
--- /dev/null
+++ b/NPSimulation/Core/EventGeneratorpBUU.hh
@@ -0,0 +1,66 @@
+#ifndef EventGeneratorpBUU_h
+#define EventGeneratorpBUU_h
+/*****************************************************************************
+ * Copyright (C) 2009-2013   this file is part of the NPTool Project         *
+ *                                                                           *
+ * For the licensing terms see $NPTOOL/Licence/NPTool_Licence                *
+ * For the list of contributors see $NPTOOL/Licence/Contributors             *
+ *****************************************************************************/
+
+/*****************************************************************************
+ * Original Author: Pierre MORFOUACE  contact address: morfouac@nscl.msu.edu       *
+ *                                                                           *
+ * Creation Date  : April 2016                                             *
+ * Last update    :                                                          *
+ *---------------------------------------------------------------------------*
+ * Decription:                                                               *
+ *  This event Generator is used to simulated ion Source from pBUU simulation         *
+ *  Very usefull to figure out Geometric Efficacity of experimental Set-Up   *
+ *---------------------------------------------------------------------------*
+ * Comment:                                                                  *
+ *                                                                           *
+ *                                                                           *
+ *****************************************************************************/
+// C++ header
+#include <string>
+using namespace std;
+
+using namespace CLHEP;
+
+// G4 headers
+#include "G4Event.hh"
+
+// NPS headers
+#include "VEventGenerator.hh"
+#include "ParticleStack.hh"
+
+// ROOT Headers
+#include "TH1F.h"
+#include "TFile.h"
+
+class EventGeneratorpBUU : public NPS::VEventGenerator{
+public:     // Constructor and destructor
+    EventGeneratorpBUU() ;
+    virtual ~EventGeneratorpBUU();
+  
+public:     // Inherit from VEventGenerator Class
+    void ReadConfiguration(string,int)               ;
+    void GenerateEvent(G4Event*) ;
+    void InitializeRootOutput()                  ;
+  
+private:    // Source parameter from input file
+    G4double               m_x0               ;  // Vertex Position X
+    G4double               m_y0               ;  // Vertex Position Y
+    G4double               m_z0               ;  // Vertex Position Z
+    G4double               m_SigmaX           ;
+    G4double               m_SigmaY           ;
+    G4ParticleDefinition*  m_particle         ;  // Kind of particle to shoot 
+    G4double               m_ExcitationEnergy ;  // Excitation energy of the emitted particle
+    string                 m_particleName     ;
+    ParticleStack*         m_ParticleStack    ;
+    
+    TH1F* fAngleHist;
+    TH1F* fEnergyHist;
+  
+};
+#endif
diff --git a/NPSimulation/Core/ParticleStack.cc b/NPSimulation/Core/ParticleStack.cc
index b93028c47..09c012db2 100644
--- a/NPSimulation/Core/ParticleStack.cc
+++ b/NPSimulation/Core/ParticleStack.cc
@@ -96,7 +96,8 @@ void ParticleStack::AddBeamParticleToStack(Particle& particle){
   m_ParticleStack.push_back(particle);
   // Incident beam parameter
   m_InitialConditions-> SetIncidentParticleName   (particle.GetParticleDefinition()->GetParticleName());
-  m_InitialConditions-> SetIncidentInitialKineticEnergy  (particle. GetParticleThetaCM());
+  //m_InitialConditions-> SetIncidentInitialKineticEnergy  (particle. GetParticleThetaCM());
+    m_InitialConditions-> SetIncidentInitialKineticEnergy  (particle. GetParticleKineticEnergy());
   
   G4ThreeVector U(1,0,0);
   G4ThreeVector V(0,1,0);
diff --git a/NPSimulation/Core/PrimaryGeneratorAction.cc b/NPSimulation/Core/PrimaryGeneratorAction.cc
index 321070127..a857be5a1 100644
--- a/NPSimulation/Core/PrimaryGeneratorAction.cc
+++ b/NPSimulation/Core/PrimaryGeneratorAction.cc
@@ -42,6 +42,7 @@
 // Event Generator Class
 #include "EventGeneratorTwoBodyReaction.hh"
 #include "EventGeneratorIsotropic.hh"
+#include "EventGeneratorpBUU.hh"
 #include "EventGeneratorBeam.hh"
 #include "EventGeneratorGammaDecay.hh"
 #include "EventGeneratorParticleDecay.hh"
@@ -79,9 +80,10 @@ void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent){
 
 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
-  bool check_Isotropic            = false;
-  bool check_TwoBodyReaction      = false;
-  bool check_Beam                 = false;
+    bool check_Isotropic            = false;
+    bool check_pBUU                 = false;
+    bool check_TwoBodyReaction      = false;
+    bool check_Beam                 = false;
 
   // You can have more than one of those
   int   alreadyiInstantiate_GammaDecay = 0;
@@ -122,6 +124,17 @@ void PrimaryGeneratorAction::ReadEventGeneratorFile(string Path){
       m_EventGenerator.push_back(myEventGenerator);
     }
 
+      //Search for pBUU source
+    else if (LineBuffer.compare(0, 4, "pBUU") == 0  && !check_pBUU) {
+        check_pBUU = true;
+        NPS::VEventGenerator* myEventGenerator = new EventGeneratorpBUU();
+        EventGeneratorFile.close();
+        myEventGenerator->ReadConfiguration(Path);
+        EventGeneratorFile.open(Path.c_str());
+        myEventGenerator->InitializeRootOutput();
+        m_EventGenerator.push_back(myEventGenerator);
+    }
+
     //Search for Beam
     else if (LineBuffer.compare(0, 4, "Beam") == 0  && !check_Beam) {
       check_Beam = true;
diff --git a/NPSimulation/PhysicsListOption.txt b/NPSimulation/PhysicsListOption.txt
index 6e6fc7c04..42fbf0b5f 100644
--- a/NPSimulation/PhysicsListOption.txt
+++ b/NPSimulation/PhysicsListOption.txt
@@ -1,10 +1,10 @@
 EmPhysicsList Option4
 DefaultCutOff 1
-IonBinaryCascadePhysics 0
+IonBinaryCascadePhysics 1
 EmExtraPhysics 0 
 HadronElasticPhysics 0
 StoppingPhysics 0
 OpticalPhysics 0
-HadronPhysicsQGSP_BIC_HP 0
-Decay 0
+HadronPhysicsQGSP_BIC_HP 1
+Decay 1
 
diff --git a/Projects/Lassa/Analysis.cxx b/Projects/Lassa/Analysis.cxx
index d8bb65d16..0c1783034 100644
--- a/Projects/Lassa/Analysis.cxx
+++ b/Projects/Lassa/Analysis.cxx
@@ -40,10 +40,18 @@ void Analysis::Init(){
     InitialConditions = new TInitialConditions();
     InitOutputBranch();
     InitInputBranch();
-    EDelta = 0.;
     totalEvents = 0;
     detectedEvents = 0;
     peakEvents = 0;
+    
+    f_proton = new TF1("f_proton","1 - TMath::Exp(-(x-0.0918309)/27.7746)",0,10);
+    f_deuton = new TF1("f_deuton","1 - TMath::Exp(-(x-0.0434552)/21.134)",0,10);
+    f_triton = new TF1("f_triton","1 - TMath::Exp(-(x-0.0353106)/17.9059)",0,10);
+    
+    //	Energy loss table: the G4Table are generated by the simulation
+    Proton_CsI = EnergyLoss("proton_CsI.G4table","G4Table",100 );
+    Deuton_CsI = EnergyLoss("deuteron_CsI.G4table","G4Table",100 );
+    Triton_CsI = EnergyLoss("triton_CsI.G4table","G4Table",100 );
 }
 
 ////////////////////////////////////////////////////////////////////////////////
@@ -53,18 +61,17 @@ void Analysis::TreatEvent(){
     
     totalEvents++;
 
-    double BeamEnergy = 120.0;
-    //double BeamEnergy = InitialConditions->GetIncidentInitialKineticEnergy();
+    
+    double ParticleEnergy = InitialConditions->GetKineticEnergy(0);
 
-    //cout << BeamEnergy << endl;
     
-    EDelta = 2.0;
+    double EDelta = 2.0;
 
-    //cout << Lassa->ThickSi_E.size() << endl;
+    if(Lassa->ThickSi_E.size()>0) InitialEnergy = ParticleEnergy;
     ///////////////////////////LOOP on Lassa Hit//////////////////////////////////
     if(Lassa->ThickSi_E.size() == 1){
-
         detectedEvents++;
+        
 
         //for(unsigned int countLassa = 0 ; countLassa < Lassa->ThickSi_E.size(); countLassa++){
         TelescopeNumber = Lassa->TelescopeNumber[0];
@@ -80,33 +87,58 @@ void Analysis::TreatEvent(){
         double Y_target = InitialConditions->GetIncidentPositionY();
         double Z_target = InitialConditions->GetIncidentPositionZ();
 
-        TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
+        TVector3 PositionOnTarget = TVector3(0,0,0);
+        //TVector3 PositionOnTarget = TVector3(X_target,Y_target,Z_target);
         TVector3 HitDirection = PositionOnLassa-PositionOnTarget;
         TVector3 HitDirectionUnit = HitDirection.Unit();
 
+        TVector3 BeamDirection = TVector3(0,0,1);
         //TVector3 BeamDirection = InitialConditions->GetBeamDirection();
         //double XBeam = BeamDirection.X();
         //double YBeam = BeamDirection.Y();
         //double ZBeam = BeamDirection.Z();
 
-        //ThetaLab = BeamDirection.Angle(HitDirection);
-        ThetaLab = 0;//ThetaLab/deg;
-
-        PhiLab = HitDirection.Phi()/deg;
+        ThetaLab = BeamDirection.Angle(HitDirection);
+        PhiLab = HitDirection.Phi();
 
         E_ThickSi = Lassa->ThickSi_E[0];
       
-
-        if(Lassa->CsI_E.size()>=1){  
-
-            E_CsI = Lassa->CsI_E[0];
+        ELab = E_ThickSi;
+        ELab_nucl = E_ThickSi;
+        
+        //R_alpha = (4.59+1.192*pow(E_ThickSi,1.724));//*cos(ThetaLab);//Lise
+        R_alpha = (4.90+1.17883*pow(E_ThickSi,1.73497));//*cos(ThetaLab);//Geant4
         
-            ELab = E_ThickSi + E_CsI;
+        if(Lassa->CsI_E.size()==1){
+            E_CsI = Lassa->CsI_E[0];
+            ELab += E_CsI;
+            //Try to simulate the nuclear reaction loss
+            //ThicknessCsI = Proton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+            //ThicknessCsI = Deuton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+            //ThicknessCsI = Triton_CsI.EvaluateMaterialThickness(0*MeV, Lassa->CsI_E[0]*MeV, 200*millimeter, 0.1*millimeter);
+            
+            //double eval = f_proton->Eval(ThicknessCsI/10);
+            //double eval = f_deuton->Eval(ThicknessCsI/10);
+            /*double eval = f_triton->Eval(ThicknessCsI/10);
+            double Random_value = Rand.Uniform(0,1);
+            
+            if(Random_value>eval) ELab_nucl += E_CsI;
+            else ELab_nucl += Rand.Uniform(0,E_CsI);*/
+            
+        }
     
-            if(ELab > (BeamEnergy-EDelta) && ELab < (BeamEnergy+EDelta) && ELab > 0.){
-                peakEvents++;
-            }
+        if(fabs(ParticleEnergy-ELab)<EDelta){
+            peakEvents++;
         }
+        else{
+            ELab = -100;
+            //E_CsI = -100;
+        }
+        
+        if(fabs(ParticleEnergy-ELab_nucl)>EDelta) ELab_nucl = -100;
+        
+        ThetaLab = ThetaLab/deg;
+        PhiLab = PhiLab/deg;
     }
 }
 ////////////////////////////////////////////////////////////////////////////////
@@ -128,13 +160,15 @@ void Analysis::End(){
 
 ////////////////////////////////////////////////////////////////////////////////
 void Analysis::InitOutputBranch() {
-    //RootOutput::getInstance()->GetTree()->Branch("Ex",&Ex,"Ex/D");
+    RootOutput::getInstance()->GetTree()->Branch("ThicknessCsI",&ThicknessCsI,"ThicknessCsI/D");
     RootOutput::getInstance()->GetTree()->Branch("ELab",&ELab,"ELab/D");
+    RootOutput::getInstance()->GetTree()->Branch("ELab_nucl",&ELab_nucl,"ELab_nucl/D");
     RootOutput::getInstance()->GetTree()->Branch("ThetaLab",&ThetaLab,"ThetaLab/D");
     RootOutput::getInstance()->GetTree()->Branch("PhiLab",&PhiLab,"PhiLab/D");
-    //  RootOutput::getInstance()->GetTree()->Branch("ThetaCM",&ThetaCM,"ThetaCM/D");
-    //  RootOutput::getInstance()->GetTree()->Branch("totalEvents",&totalEvents,"totalEvents/I");
-    //  RootOutput::getInstance()->GetTree()->Branch("detectedEvents",&detectedEvents,"detectedEvents/I");
+    RootOutput::getInstance()->GetTree()->Branch("InitialEnergy",&InitialEnergy,"InitialEnergy/D");
+    RootOutput::getInstance()->GetTree()->Branch("E_ThickSi",&E_ThickSi,"E_ThickSi/D");
+    RootOutput::getInstance()->GetTree()->Branch("E_CsI",&E_CsI,"E_CsI/D");
+    RootOutput::getInstance()->GetTree()->Branch("R_alpha",&R_alpha,"R_alpha/D");
     //  RootOutput::getInstance()->GetTree()->Branch("peakEvents",&peakEvents,"peakEvents/I");
 }
 
@@ -147,15 +181,19 @@ void Analysis::InitInputBranch(){
 
 ////////////////////////////////////////////////////////////////////////////////     
 void Analysis::ReInitValue(){
-    E_ThickSi = 0.;
-    E_CsI = 0.;
-    ELab = 0.;
-    ThetaLab = -1000;
-    PhiLab = -1000;
-    X = -1000;
-    Y = -1000;
-    Z = -1000;
+    E_ThickSi = -100;
+    E_CsI =-100;
+    ELab = -100;
+    ELab_nucl = -100;
+    ThetaLab = -100;
+    PhiLab = -100;
+    X = -100;
+    Y = -100;
+    Z = -100;
     TelescopeNumber = -1;
+    InitialEnergy = -100;
+    ThicknessCsI = -1;
+    R_alpha = -100;
 }
 
 ////////////////////////////////////////////////////////////////////////////////
diff --git a/Projects/Lassa/Analysis.h b/Projects/Lassa/Analysis.h
index a4de214c7..c60f8cb4d 100644
--- a/Projects/Lassa/Analysis.h
+++ b/Projects/Lassa/Analysis.h
@@ -27,6 +27,7 @@
 #include "NPEnergyLoss.h"
 #include "NPReaction.h"
 #include "TRandom3.h"
+#include "TF1.h"
 class Analysis: public NPL::VAnalysis{
   public:
     Analysis();
@@ -43,14 +44,17 @@ class Analysis: public NPL::VAnalysis{
 
   private:
     double ELab;
+    double ELab_nucl;
     double E_ThickSi;
     double E_CsI;
-    double EDelta;
     double PhiLab;
     double ThetaLab;
     double X,Y,Z;
     double TelescopeNumber;
     double thresholdEnergy;
+    double InitialEnergy;
+    double ThicknessCsI;
+    double R_alpha;
 
     int totalEvents;
     int detectedEvents;
@@ -62,6 +66,15 @@ class Analysis: public NPL::VAnalysis{
     // intermediate variable
     TRandom3 Rand;
     
+    TF1* f_proton;
+    TF1* f_deuton;
+    TF1* f_triton;
+    
+    //Energy loss table
+    NPL::EnergyLoss Proton_CsI;
+    NPL::EnergyLoss Deuton_CsI;
+    NPL::EnergyLoss Triton_CsI;
+    
     TLassaPhysics* Lassa;
 	TInitialConditions* InitialConditions;
 };
diff --git a/Projects/Lassa/RunToTreat.txt b/Projects/Lassa/RunToTreat.txt
index a50a6ea67..599a23270 100755
--- a/Projects/Lassa/RunToTreat.txt
+++ b/Projects/Lassa/RunToTreat.txt
@@ -1,5 +1,7 @@
 TTreeName 
 	SimulatedTree
 RootFileName 
-	../../Outputs/Simulation/dumb.root
+  ../../Outputs/Simulation/lassa_pid_p.root
+  ../../Outputs/Simulation/lassa_pid_d.root
+	../../Outputs/Simulation/lassa_pid_t.root
 
diff --git a/Projects/macros/GeometricalEfficiency.C b/Projects/macros/GeometricalEfficiency.C
index 843843e63..b39ffb055 100644
--- a/Projects/macros/GeometricalEfficiency.C
+++ b/Projects/macros/GeometricalEfficiency.C
@@ -48,7 +48,7 @@
 using namespace std;
 
 
-void GeometricalEfficiency(const char * fname = "myResult"){
+void GeometricalEfficiency(const char * fname = "lassa_proton"){
   // Open output ROOT file from NPTool simulation run
   TString path = gSystem->Getenv("NPTOOL");
   path += "/Outputs/Simulation/";
-- 
GitLab