00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037
00038
00039
00040 #include "RooFit.h"
00041
00042 #include <math.h>
00043
00044
00045 #include "RooBukinPdf.h"
00046 #include "RooRealVar.h"
00047 #include "TMath.h"
00048
00049 ClassImp(RooBukinPdf)
00050
00051
00052
00053
00054 RooBukinPdf::RooBukinPdf(const char *name, const char *title,
00055 RooAbsReal& _x, RooAbsReal& _Xp,
00056 RooAbsReal& _sigp, RooAbsReal& _xi,
00057 RooAbsReal& _rho1, RooAbsReal& _rho2) :
00058
00059
00060 RooAbsPdf(name, title),
00061 x("x","x",this,_x),
00062 Xp("Xp","Xp",this,_Xp),
00063 sigp("sigp","sigp",this,_sigp),
00064 xi("xi","xi",this,_xi),
00065 rho1("rho1","rho1",this,_rho1),
00066 rho2("rho2","rho2",this,_rho2)
00067 {
00068
00069 consts = 2*sqrt(2*log(2.));
00070 }
00071
00072
00073
00074
00075 RooBukinPdf::RooBukinPdf(const RooBukinPdf& other, const char *name):
00076 RooAbsPdf(other,name),
00077 x("x",this,other.x),
00078 Xp("Xp",this,other.Xp),
00079 sigp("sigp",this,other.sigp),
00080 xi("xi",this,other.xi),
00081 rho1("rho1",this,other.rho1),
00082 rho2("rho2",this,other.rho2)
00083
00084 {
00085
00086 consts = 2*sqrt(2*log(2.));
00087 }
00088
00089
00090
00091
00092 Double_t RooBukinPdf::evaluate() const
00093 {
00094
00095
00096 double r1=0,r2=0,r3=0,r4=0,r5=0,hp=0;
00097 double x1 = 0,x2 = 0;
00098 double fit_result = 0;
00099
00100 hp=sigp*consts;
00101 r3=log(2.);
00102 r4=sqrt(TMath::Power(xi,2)+1);
00103 r1=xi/r4;
00104
00105 if(TMath::Abs(xi) > exp(-6.)){
00106 r5=xi/log(r4+xi);
00107 }
00108 else
00109 r5=1;
00110
00111 x1 = Xp + (hp / 2) * (r1-1);
00112 x2 = Xp + (hp / 2) * (r1+1);
00113
00114
00115 if(x < x1){
00116 r2=rho1*TMath::Power((x-x1)/(Xp-x1),2)-r3 + 4 * r3 * (x-x1)/hp * r5 * r4/TMath::Power((r4-xi),2);
00117 }
00118
00119
00120
00121 else if(x < x2) {
00122 if(TMath::Abs(xi) > exp(-6.)) {
00123 r2=log(1 + 4 * xi * r4 * (x-Xp)/hp)/log(1+2*xi*(xi-r4));
00124 r2=-r3*(TMath::Power(r2,2));
00125 }
00126 else{
00127 r2=-4*r3*TMath::Power(((x-Xp)/hp),2);
00128 }
00129 }
00130
00131
00132
00133 else {
00134 r2=rho2*TMath::Power((x-x2)/(Xp-x2),2)-r3 - 4 * r3 * (x-x2)/hp * r5 * r4/TMath::Power((r4+xi),2);
00135 }
00136
00137 if(TMath::Abs(r2) > 100){
00138 fit_result = 0;
00139 }
00140 else{
00141
00142 fit_result = exp(r2);
00143 }
00144
00145 return fit_result;
00146
00147 }