00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016 #include "EvtGenBase/EvtPatches.hh"
00017
00018
00019
00020 #include <assert.h>
00021
00022 #include "EvtGenModels/EvtBtoKD3P.hh"
00023 #include "EvtGenBase/EvtDecayTable.hh"
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtId.hh"
00026 #include "EvtGenBase/EvtRandom.hh"
00027 #include "EvtGenModels/EvtPto3P.hh"
00028
00029 #include "EvtGenBase/EvtDalitzPoint.hh"
00030 #include "EvtGenBase/EvtCyclic3.hh"
00031 using std::endl;
00032
00033
00034 EvtBtoKD3P::EvtBtoKD3P() :
00035 _model1(0),
00036 _model2(0),
00037 _decayedOnce(false)
00038 {
00039 }
00040
00041
00042 EvtBtoKD3P::EvtBtoKD3P(const EvtBtoKD3P & other){
00043 }
00044
00045
00046 EvtBtoKD3P::~EvtBtoKD3P(){
00047 }
00048
00049
00050 EvtDecayBase * EvtBtoKD3P::clone(){
00051 return new EvtBtoKD3P();
00052 }
00053
00054
00055 void EvtBtoKD3P::getName(std::string& model_name){
00056 model_name="BTOKD3P";
00057 }
00058
00059
00060 void EvtBtoKD3P::init(){
00061 checkNArg(2);
00062 checkNDaug(3);
00063
00064
00065
00066 checkSpinParent ( EvtSpinType::SCALAR);
00067 checkSpinDaughter(0,EvtSpinType::SCALAR);
00068 checkSpinDaughter(1,EvtSpinType::SCALAR);
00069 checkSpinDaughter(2,EvtSpinType::SCALAR);
00070
00071
00072
00073
00074 _r = getArg(0);
00075 double phase = getArg(1);
00076 _exp = EvtComplex(cos(phase), sin(phase));
00077 }
00078
00079
00080 void EvtBtoKD3P::initProbMax(){
00081 setProbMax(1);
00082 }
00083
00084
00085 void EvtBtoKD3P::decay(EvtParticle *p){
00086
00087 _daugsDecayedByParentModel = true;
00088
00089
00090
00091
00092 const int KIND = 0;
00093 const int D1IND = 1;
00094 const int D2IND = 2;
00095
00096
00097 EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)};
00098 p->initializePhaseSpace(2, tempDaug);
00099
00100
00101
00102 EvtParticle * theD = p->getDaug(D1IND);
00103 EvtPto3P * model1 = (EvtPto3P*)(EvtDecayTable::getDecayFunc(theD));
00104
00105
00106 theD->init(getDaug(D2IND), theD->getP4());
00107 EvtPto3P * model2 = (EvtPto3P*)(EvtDecayTable::getDecayFunc(theD));
00108
00109
00110 if (false == _decayedOnce) {
00111 _decayedOnce = true;
00112
00113
00114 _model1 = model1;
00115 _model2 = model2;
00116
00117
00118
00119 std::string name1;
00120 std::string name2;
00121 model1->getName(name1);
00122 model2->getName(name2);
00123
00124 if (name1 != "PTO3P") {
00125 report(ERROR,"EvtGen")
00126 << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model"
00127 << endl
00128 << " but found to decay via " << name1.c_str()
00129 << " or " << name2.c_str()
00130 << ". Will terminate execution!" << endl;
00131 assert(0);
00132 }
00133
00134 EvtId * daugs1 = model1->getDaugs();
00135 EvtId * daugs2 = model2->getDaugs();
00136
00137 bool idMatch = true;
00138 int d;
00139 for (d = 0; d < 2; ++d) {
00140 if (daugs1[d] != daugs2[d]) {
00141 idMatch = false;
00142 }
00143 }
00144 if (false == idMatch) {
00145 report(ERROR,"EvtGen")
00146 << "D daughters of EvtBtoKD3P decay must decay to the same final state"
00147 << endl
00148 << " particles in the same order (not CP-conjugate order)," << endl
00149 << " but they were found to decay to" << endl;
00150 for (d = 0; d < model1->getNDaug(); ++d) {
00151 report(ERROR,"") << " " << EvtPDL::name(daugs1[d]).c_str() << " ";
00152 }
00153 report(ERROR,"") << endl;
00154 for (d = 0; d < model1->getNDaug(); ++d) {
00155 report(ERROR,"") << " " << EvtPDL::name(daugs2[d]).c_str() << " ";
00156 }
00157 report(ERROR,"") << endl << ". Will terminate execution!" << endl;
00158 assert(0);
00159 }
00160
00161
00162
00163 setProbMax(model1->getProbMax(0)
00164 + _r * _r * model2->getProbMax(0)
00165 + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0)));
00166
00167 }
00168
00169
00170 if (_model1 != model1 || _model2 != model2) {
00171 report(ERROR,"EvtGen")
00172 << "D daughters of EvtBtoKD3P decay should have only 1 decay modes, "
00173 << endl
00174 << " but a new decay mode was found after the first call" << endl
00175 << " Will terminate execution!" << endl;
00176 assert(0);
00177 }
00178
00179
00180
00181
00182
00183
00184 EvtPdfSum<EvtDalitzPoint> * pc1 = model1->getPC();
00185 EvtPdfSum<EvtDalitzPoint> * pc2 = model2->getPC();
00186 EvtPdfSum<EvtDalitzPoint> pc;
00187 pc.addTerm(1.0, *pc1);
00188 pc.addTerm(1.0, *pc2);
00189
00190
00191 EvtDalitzPoint x = pc.randomPoint();
00192
00193
00194 EvtComplex amp1 = model1->amplNonCP(x);
00195 EvtComplex amp2 = model2->amplNonCP(x);
00196 EvtComplex amp = amp1 + amp2 * _r * _exp;
00197
00198
00199
00200
00201 double comp = sqrt(pc.evaluate (x));
00202 vertex (amp/comp);
00203
00204
00205 theD->generateMassTree();
00206
00207
00208 std::vector<EvtVector4R> v = model2->initDaughters(x);
00209
00210 if(v.size() != theD->getNDaug()) {
00211 report(ERROR,"EvtGen")
00212 << "Number of daughters " << theD->getNDaug()
00213 << " != " << "Momentum vector size " << v.size()
00214 << endl
00215 << " Terminating execution." << endl;
00216 assert(0);
00217 }
00218
00219
00220 int i;
00221 for(i=0; i<theD->getNDaug(); ++i){
00222 theD->getDaug(i)->init(model2->getDaugs()[i], v[i]);
00223 }
00224 }
00225