/***************************************************************************** * Project: RooFit * * * * This code was autogenerated by RooClassFactory * *****************************************************************************/ // Your description goes here... #include "Riostream.h" #include "PndLmdcs_th.h" #include "RooAbsReal.h" #include "RooAbsCategory.h" #include #include "TMath.h" ClassImp(cs_th) cs_th::cs_th(const char *name, const char *title, RooAbsReal& _th, // RooAbsReal& _L, RooAbsReal& _Sigtot, RooAbsReal& _B, RooAbsReal& _Rho, double _Plab) : RooAbsPdf(name,title), th("th","th",this,_th), // L("L","L",this,_L), Sigtot("Sigtot","Sigtot",this,_Sigtot), B("B","B",this,_B), Rho("Rho","Rho",this,_Rho), Plab(_Plab) { M=0.938; //Plab=1.5; Elab=sqrt(M*M+Plab*Plab); S=2*M*M+2*M*Elab; pcm2=S/4.-M*M; gamma=(Elab+M)/sqrt(S); beta=Plab/Elab; lnP=log(Plab); pi=3.14159; hbarc2=0.389379; alpha=1./137.036; pkoef=10./(5.0677*5.0677); dt=1.0802e-5; } cs_th::cs_th(const cs_th& other, const char* name) : RooAbsPdf(other,name), th("th",this,other.th), // L("L",this,other.L), Sigtot("Sigtot",this,other.Sigtot), B("B",this,other.B), Rho("Rho",this,other.Rho), Plab(other.Plab) { M=0.938; //Plab=1.5; Elab=sqrt(M*M+Plab*Plab); S=2*M*M+2*M*Elab; pcm2=S/4.-M*M; gamma=(Elab+M)/sqrt(S); beta=Plab/Elab; lnP=log(Plab); pi=3.14159; hbarc2=0.389379; alpha=1./137.036; pkoef=10./(5.0677*5.0677); dt=1.0802e-5; } Double_t cs_th::evaluate() const { /*Double_t t = -1*x[0] ; Double_t sigtot = par[0]; Double_t B = par[1]; Double_t Rho = par[2]; Double_t L = par[3];*/ if(th<1.7) return 0; //************************************************************************************* //transform polar-angle th to momentum-transfer t //************************************************************************************* Double_t thetacm = 2*TMath::ATan(gamma*TMath::Tan(th/1000.)); Double_t t = 2*pcm2*(1-cos(thetacm)); Double_t jaco = 8.*gamma*gamma*pcm2*TMath::Tan(th/1000.)*4./TMath::Power(TMath::Cos(2.*th/1000.)+1.,2)/TMath::Power(1.+gamma*gamma*TMath::Tan(th/1000.)*TMath::Tan(th/1000.),2); //************************************************************************************* Double_t p1,p2,p3,dndt; Double_t G = 1/(1+TMath::Abs(t)/0.71)/(1+TMath::Abs(t)/0.71); Double_t G2 = G*G; Double_t G4 = G2*G2; Double_t del = TMath::Abs(t)/0.71; Double_t dd1 = 0.5*B*TMath::Abs(t)+4*del, dd2 = 4*del; Double_t logdd1 = TMath::Log(dd1), logdd2 = TMath::Log(dd2); Double_t delta = alpha*(0.577+logdd1+4*del*logdd2+2*del); //************************************************************************************* //p1 and p2 are Coulomb parts and interference parts //************************************************************************************* p1 = 4*pi*alpha*alpha*G4*hbarc2/(beta*beta)/(t*t);//Coulomb part p2 = alpha*Sigtot*G2*exp(-0.5*B*TMath::Abs(t))*(Rho*cos(delta)+sin(delta))/beta/TMath::Abs(t);//iExact // p2 = alpha*Sigtot*G2*exp(-0.5*B*TMath::Abs(t))*TMath::Sqrt(1+Rho*Rho)/beta/TMath::Abs(t); //Inter called by DPM //************************************************************************** //Parameterized formula for Hadron parts in DPM //*************************************************************************** p3 = Sigtot*Sigtot*(1+Rho*Rho)*exp(-1*B*TMath::Abs(t))/16./pi/hbarc2;//E760, typical dndt = (p1+p2+p3)*jaco*dt*2; return dndt; //return 1.0 ; } Double_t cs_th::evaluate(Double_t x) const { if(th<1.7) return 0; //************************************************************************************* //transform polar-angle th to momentum-transfer t //************************************************************************************* Double_t thetacm = 2*TMath::ATan(gamma*TMath::Tan(x/1000.)); Double_t t = 2*pcm2*(1-cos(thetacm)); Double_t jaco = 8.*gamma*gamma*pcm2*TMath::Tan(th/1000.)*4./TMath::Power(TMath::Cos(2.*th/1000.)+1.,2)/TMath::Power(1.+gamma*gamma*TMath::Tan(th/1000.)*TMath::Tan(th/1000.),2); //************************************************************************************* Double_t p1,p2,p3,dndt; Double_t G = 1/(1+TMath::Abs(t)/0.71)/(1+TMath::Abs(t)/0.71); Double_t G2 = G*G; Double_t G4 = G2*G2; Double_t del = TMath::Abs(t)/0.71; Double_t dd1 = 0.5*B*TMath::Abs(t)+4*del, dd2 = 4*del; Double_t logdd1 = TMath::Log(dd1), logdd2 = TMath::Log(dd2); Double_t delta = alpha*(0.577+logdd1+4*del*logdd2+2*del); //************************************************************************************* //p1 and p2 are Coulomb parts and interference parts //************************************************************************************* p1 = 4*pi*alpha*alpha*G4*hbarc2/(beta*beta)/(t*t);//Coulomb part p2 = alpha*Sigtot*G2*exp(-0.5*B*TMath::Abs(t))*(Rho*cos(delta)+sin(delta))/beta/TMath::Abs(t);//iExact // p2 = alpha*Sigtot*G2*exp(-0.5*B*TMath::Abs(t))*TMath::Sqrt(1+Rho*Rho)/beta/TMath::Abs(t); //Inter called by DPM //************************************************************************** //Parameterized formula for Hadron parts in DPM //*************************************************************************** p3 = Sigtot*Sigtot*(1+Rho*Rho)*exp(-1*B*TMath::Abs(t))/16./pi/hbarc2;//E760, typical dndt = (p1+p2+p3)*jaco*dt*2; return dndt; } /* Int_t cs_th::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char*) const { // LIST HERE OVER WHICH VARIABLES ANALYTICAL INTEGRATION IS SUPPORTED, // ASSIGN A NUMERIC CODE FOR EACH SUPPORTED (SET OF) PARAMETERS // THE EXAMPLE BELOW ASSIGNS CODE 1 TO INTEGRATION OVER VARIABLE X // YOU CAN ALSO IMPLEMENT MORE THAN ONE ANALYTICAL INTEGRAL BY REPEATING THE matchArgs // EXPRESSION MULTIPLE TIMES if (matchArgs(allVars,analVars,th)) return 1 ; return 0 ; } Double_t cs_th::analyticalIntegral(Int_t code, const char* rangeName) const { // RETURN ANALYTICAL INTEGRAL DEFINED BY RETURN CODE ASSIGNED BY getAnalyticalIntegral // THE MEMBER FUNCTION x.min(rangeName) AND x.max(rangeName) WILL RETURN THE INTEGRATION // BOUNDARIES FOR EACH OBSERVABLE x // assert(code==1) ; // return (x.max(rangeName)-x.min(rangeName)) ; return 1. ; } */