00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023 #include "EvtGenBase/EvtPatches.hh"
00024
00025 #include <stdlib.h>
00026 #include "EvtGenBase/EvtRandom.hh"
00027 #include "EvtGenBase/EvtParticle.hh"
00028 #include "EvtGenBase/EvtGenKine.hh"
00029 #include "EvtGenBase/EvtPDL.hh"
00030 #include "EvtGenBase/EvtReport.hh"
00031 #include "EvtGenModels/EvtBtoXsEtap.hh"
00032 #include <string>
00033 #include "EvtGenBase/EvtConst.hh"
00034 using std::endl;
00035
00036 EvtBtoXsEtap::~EvtBtoXsEtap() {}
00037
00038 void EvtBtoXsEtap::getName(std::string& model_name){
00039
00040 model_name="BTOXSETAP";
00041
00042 }
00043
00044 EvtDecayBase* EvtBtoXsEtap::clone(){
00045
00046 return new EvtBtoXsEtap;
00047
00048 }
00049
00050
00051 void EvtBtoXsEtap::init(){
00052
00053
00054
00055 checkNArg(0);
00056 }
00057
00058 void EvtBtoXsEtap::initProbMax(){
00059
00060 noProbMax();
00061
00062 }
00063
00064 void EvtBtoXsEtap::decay( EvtParticle *p ){
00065
00066
00067
00068
00069
00070
00071
00072
00073 double m_b;
00074 int i;
00075 p->makeDaughters(getNDaug(),getDaugs());
00076 EvtParticle *pdaug[MAX_DAUG];
00077
00078 for(i=0;i<getNDaug();i++){
00079 pdaug[i]=p->getDaug(i);
00080 }
00081
00082 static EvtVector4R p4[MAX_DAUG];
00083 static double mass[MAX_DAUG];
00084
00085 m_b = p->mass();
00086
00087
00088
00089 mass[1] = EvtPDL::getMass(getDaug(1));
00090
00091 double xbox, ybox, min, max,hichfit;
00092 min=0.493;
00093 max=4.3;
00094 const double TwoPi = EvtConst::twoPi;
00095 int Xscode = EvtPDL::getStdHep(getDaug(0));
00096
00097
00098
00099
00100 double par[6];
00101 if ((Xscode == 30343) || (Xscode == -30343) ||
00102 (Xscode == 30353) || (Xscode == -30353)) {
00103 min=0.6373;
00104
00105
00106 par[0]=2.36816;
00107
00108 par[1]=0.62325725;
00109
00110 par[2]=2.2;
00111
00112 par[3]=-0.2109375;
00113
00114 par[4]=2.7;
00115
00116 par[5]=0.54;
00117
00118
00119
00120
00121
00122
00123
00124
00125
00126
00127
00128
00129 } else {
00130 report(DEBUG,"EvtGen") << "In EvtBtoXsEtap: Particle with id " << Xscode << " is not a Xsd/u particle"<<endl;
00131 return;
00132 }
00133
00134 double boxheight=par[5];
00135 double boxwidth=max-min;
00136
00137 mass[0]=0.0;
00138 while ((mass[0] > max) || (mass[0] < min)){
00139 xbox = EvtRandom::Flat(boxwidth)+min;
00140 ybox=EvtRandom::Flat(boxheight);
00141 if (xbox<par[2]) {
00142
00143 hichfit=(1/sqrt(TwoPi*par[1]))*exp(-0.5*pow((xbox-par[0])/par[1],2));
00144
00145
00146
00147
00148
00149 } else {
00150 hichfit=par[3]*pow((xbox-par[4]),2)+par[5];
00151
00152 }
00153 if (ybox>hichfit) {
00154 mass[0]=0.0;
00155 } else {
00156 mass[0]=xbox;
00157 }
00158 }
00159
00160
00161
00162 EvtGenKine::PhaseSpace( getNDaug(), mass, p4, m_b );
00163
00164 for(i=0;i<getNDaug();i++){
00165 pdaug[i]->init( getDaugs()[i], p4[i] );
00166 }
00167
00168 return ;
00169 }
00170
00171