1 //-----------------------------------------------------------------------
2 // File and Version Information:
3 // $Id: EvtBtoKD3P.cc,v 1.7 2009/02/19 03:23:29 ryd Exp $
6 // This software is part of the EvtGen package developed jointly
7 // for the BaBar and CLEO collaborations. If you use all or part
8 // of it, please give an appropriate acknowledgement.
10 // Copyright Information:
11 // Copyright (C) 2003, Colorado State University
14 // Abi soffer, CSU, 2003
15 //-----------------------------------------------------------------------
16 #include "EvtGenBase/EvtPatches.hh"
18 // Decay model that does the decay B+->D0K, D0->3 psudoscalars
22 #include "EvtGenModels/EvtBtoKD3P.hh"
23 #include "EvtGenBase/EvtDecayTable.hh"
24 #include "EvtGenBase/EvtParticle.hh"
25 #include "EvtGenBase/EvtId.hh"
26 #include "EvtGenBase/EvtRandom.hh"
27 #include "EvtGenModels/EvtPto3P.hh"
28 #include "EvtGenBase/EvtReport.hh"
29 #include "EvtGenBase/EvtDalitzPoint.hh"
30 #include "EvtGenBase/EvtCyclic3.hh"
33 //------------------------------------------------------------------
34 EvtBtoKD3P::EvtBtoKD3P() :
41 //------------------------------------------------------------------
42 EvtBtoKD3P::EvtBtoKD3P(const EvtBtoKD3P & other){
45 //------------------------------------------------------------------
46 EvtBtoKD3P::~EvtBtoKD3P(){
49 //------------------------------------------------------------------
50 EvtDecayBase * EvtBtoKD3P::clone(){
51 return new EvtBtoKD3P();
54 //------------------------------------------------------------------
55 std::string EvtBtoKD3P::getName(){
59 //------------------------------------------------------------------
60 void EvtBtoKD3P::init(){
61 checkNArg(2); // r, phase
62 checkNDaug(3); // K, D0(allowed), D0(suppressed).
63 // The last two daughters are really one particle
65 // check that the mother and all daughters are scalars:
66 checkSpinParent ( EvtSpinType::SCALAR);
67 checkSpinDaughter(0,EvtSpinType::SCALAR);
68 checkSpinDaughter(1,EvtSpinType::SCALAR);
69 checkSpinDaughter(2,EvtSpinType::SCALAR);
71 // Check that the B dtr types are K D D:
73 // get the parameters:
75 double phase = getArg(1);
76 _exp = EvtComplex(cos(phase), sin(phase));
79 //------------------------------------------------------------------
80 void EvtBtoKD3P::initProbMax(){
81 setProbMax(1); // this is later changed in decay()
84 //------------------------------------------------------------------
85 void EvtBtoKD3P::decay(EvtParticle *p){
86 // tell the subclass that we decay the daughter:
87 _daugsDecayedByParentModel = true;
89 // the K is the 1st daughter of the B EvtParticle.
90 // The decay mode of the allowed D (the one produced in b->c decay) is 2nd
91 // The decay mode of the suppressed D (the one produced in b->u decay) is 3rd
96 // generate kinematics of daughters (K and D):
97 EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)};
98 p->initializePhaseSpace(2, tempDaug);
100 // Get the D daughter particle and the decay models of the allowed
101 // and suppressed D modes:
102 EvtParticle * theD = p->getDaug(D1IND);
103 EvtPto3P * model1 = (EvtPto3P*)(EvtDecayTable::getDecayFunc(theD));
105 // for the suppressed mode, re-initialize theD as the suppressed D alias:
106 theD->init(getDaug(D2IND), theD->getP4());
107 EvtPto3P * model2 = (EvtPto3P*)(EvtDecayTable::getDecayFunc(theD));
109 // on the first call:
110 if (false == _decayedOnce) {
113 // store the D decay model pointers:
117 // check the decay models of the first 2 daughters and that they
118 // have the same final states:
119 std::string name1=model1->getName();
120 std::string name2=model2->getName();
122 if (name1 != "PTO3P") {
123 report(ERROR,"EvtGen")
124 << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model"
126 << " but found to decay via " << name1.c_str()
127 << " or " << name2.c_str()
128 << ". Will terminate execution!" << endl;
132 EvtId * daugs1 = model1->getDaugs();
133 EvtId * daugs2 = model2->getDaugs();
137 for (d = 0; d < 2; ++d) {
138 if (daugs1[d] != daugs2[d]) {
142 if (false == idMatch) {
143 report(ERROR,"EvtGen")
144 << "D daughters of EvtBtoKD3P decay must decay to the same final state"
146 << " particles in the same order (not CP-conjugate order)," << endl
147 << " but they were found to decay to" << endl;
148 for (d = 0; d < model1->getNDaug(); ++d) {
149 report(ERROR,"") << " " << EvtPDL::name(daugs1[d]).c_str() << " ";
151 report(ERROR,"") << endl;
152 for (d = 0; d < model1->getNDaug(); ++d) {
153 report(ERROR,"") << " " << EvtPDL::name(daugs2[d]).c_str() << " ";
155 report(ERROR,"") << endl << ". Will terminate execution!" << endl;
159 // estimate the probmax. Need to know the probmax's of the 2
161 setProbMax(model1->getProbMax(0)
162 + _r * _r * model2->getProbMax(0)
163 + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0)));
165 } // end of things to do on the first call
167 // make sure the models haven't changed since the first call:
168 if (_model1 != model1 || _model2 != model2) {
169 report(ERROR,"EvtGen")
170 << "D daughters of EvtBtoKD3P decay should have only 1 decay modes, "
172 << " but a new decay mode was found after the first call" << endl
173 << " Will terminate execution!" << endl;
177 // get the cover function for each of the models and add them up.
178 // They are summed with coefficients 1 because we are willing to
179 // take a small inefficiency (~50%) in order to ensure that the
180 // cover function is large enough without getting into complications
181 // associated with the smallness of _r:
182 EvtPdfSum<EvtDalitzPoint> * pc1 = model1->getPC();
183 EvtPdfSum<EvtDalitzPoint> * pc2 = model2->getPC();
184 EvtPdfSum<EvtDalitzPoint> pc;
185 pc.addTerm(1.0, *pc1);
186 pc.addTerm(1.0, *pc2);
188 // from this combined cover function, generate the Dalitz point:
189 EvtDalitzPoint x = pc.randomPoint();
191 // get the aptitude for each of the models on this point and add them up:
192 EvtComplex amp1 = model1->amplNonCP(x);
193 EvtComplex amp2 = model2->amplNonCP(x);
194 EvtComplex amp = amp1 + amp2 * _r * _exp;
196 // get the value of the cover function for this point and set the
197 // relative amplitude for this decay:
199 double comp = sqrt(pc.evaluate (x));
202 // Make the daughters of theD:
203 theD->generateMassTree();
205 // Now generate the p4's of the daughters of theD:
206 std::vector<EvtVector4R> v = model2->initDaughters(x);
208 if(v.size() != theD->getNDaug()) {
209 report(ERROR,"EvtGen")
210 << "Number of daughters " << theD->getNDaug()
211 << " != " << "Momentum vector size " << v.size()
213 << " Terminating execution." << endl;
217 // Apply the new p4's to the daughters:
218 for(unsigned int i=0; i<theD->getNDaug(); ++i){
219 theD->getDaug(i)->init(model2->getDaugs()[i], v[i]);