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 #include "EvtGenBase/EvtPatches.hh"
00034
00035 #include <stdlib.h>
00036 #include "EvtGenBase/EvtRandom.hh"
00037 #include "EvtGenModels/EvtBtoXsgammaAliGreub.hh"
00038 #include <string>
00039 #include "EvtGenBase/EvtConst.hh"
00040 #include "EvtGenBase/EvtParticle.hh"
00041 #include "EvtGenBase/EvtGenKine.hh"
00042 #include "EvtGenBase/EvtPDL.hh"
00043 #include "EvtGenBase/EvtReport.hh"
00044 using std::endl;
00045
00046
00047 EvtBtoXsgammaAliGreub::~EvtBtoXsgammaAliGreub(){}
00048
00049 void EvtBtoXsgammaAliGreub::init(int nArg, double* args){
00050
00051 if ((nArg - 1) != 0) {
00052
00053 report(ERROR,"EvtGen") << "EvtBtoXsgamma generator model "
00054 << "EvtBtoXsgammaAliGreub expected "
00055 << "zero arguments but found: "<<nArg-1<<endl;
00056 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
00057 ::abort();
00058
00059 }
00060
00061 }
00062
00063 double EvtBtoXsgammaAliGreub::GetMass(int Xscode) {
00064
00065
00066
00067
00068
00069
00070
00071
00072 double min=0.64;
00073 double max=4.5;
00074 double xbox, ybox, alifit;
00075 double mass=0.0;
00076
00077 double par[18];
00078 if ((Xscode == 30343) || (Xscode == -30343) ||
00079 (Xscode == 30353) || (Xscode == -30353)) {
00080 min=0.6373;
00081
00082 par[0]=-2057.2380371094;
00083 par[1]=2502.2556152344;
00084 par[2]=1151.5632324219;
00085 par[3]=0.82431584596634;
00086 par[4]=-4110.5234375000;
00087 par[5]=8445.6757812500;
00088 par[6]=-3034.1894531250;
00089 par[7]=1.1557708978653;
00090 par[8]=1765.9311523438;
00091 par[9]=1.3730158805847;
00092 par[10]=0.51371538639069;
00093 par[11]=2.0056934356689;
00094 par[12]=37144.097656250;
00095 par[13]=-50296.781250000;
00096 par[14]=27319.095703125;
00097 par[15]=-7408.0678710938;
00098 par[16]=1000.8093261719;
00099 par[17]=-53.834449768066;
00100 } else if ((Xscode == 30363) || (Xscode == -30363)) {
00101 min = 0.9964;
00102 par[0]=-32263.908203125;
00103 par[1]=57186.589843750;
00104 par[2]=-24230.728515625;
00105 par[3]=1.1155973672867;
00106 par[4]=-12161.131835938;
00107 par[5]=20162.146484375;
00108 par[6]=-7198.8564453125;
00109 par[7]=1.3783323764801;
00110 par[8]=1995.1691894531;
00111 par[9]=1.4655895233154;
00112 par[10]=0.48869228363037;
00113 par[11]=2.1038570404053;
00114 par[12]=55100.058593750;
00115 par[13]=-75201.703125000;
00116 par[14]=41096.066406250;
00117 par[15]=-11205.986328125;
00118 par[16]=1522.4024658203;
00119 par[17]=-82.379623413086;
00120 } else {
00121 report(DEBUG,"EvtGen") << "In EvtBtoXsgammaAliGreub: Particle with id " << Xscode << " is not a Xss particle"<<endl;
00122 return 0;
00123 }
00124
00125 double boxheight=par[8];
00126 double boxwidth=max-min;
00127
00128 while ((mass > max) || (mass < min)){
00129 xbox = EvtRandom::Flat(boxwidth)+min;
00130 ybox=EvtRandom::Flat(boxheight);
00131 if (xbox<par[3]) {
00132 alifit=par[0]+par[1]*xbox+par[2]*pow(xbox,2);
00133 } else if (xbox<par[7]) {
00134 alifit=par[4]+par[5]*xbox+par[6]*pow(xbox,2);
00135 } else if (xbox<par[11]) {
00136 alifit=par[8]*exp(-0.5*pow((xbox-par[9])/par[10],2));
00137 } else {
00138 alifit=par[12]+par[13]*xbox+par[14]*pow(xbox,2)+par[15]*pow(xbox,3)+par[16]*pow(xbox,4)+par[17]*pow(xbox,5);
00139 }
00140 if (ybox>alifit) {
00141 mass=0.0;
00142 } else {
00143 mass=xbox;
00144 }
00145 }
00146 return mass;
00147 }