1 //-----------------------------------------------------------------------
2 // File and Version Information:
3 // $Id: EvtBtoKD3P.cpp,v 1.2 2009-04-02 15:22:28 robbep 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) :
46 //------------------------------------------------------------------
47 EvtBtoKD3P::~EvtBtoKD3P(){
50 //------------------------------------------------------------------
51 EvtDecayBase * EvtBtoKD3P::clone(){
52 return new EvtBtoKD3P();
55 //------------------------------------------------------------------
56 std::string EvtBtoKD3P::getName(){
60 //------------------------------------------------------------------
61 void EvtBtoKD3P::init(){
62 checkNArg(2); // r, phase
63 checkNDaug(3); // K, D0(allowed), D0(suppressed).
64 // The last two daughters are really one particle
66 // check that the mother and all daughters are scalars:
67 checkSpinParent ( EvtSpinType::SCALAR);
68 checkSpinDaughter(0,EvtSpinType::SCALAR);
69 checkSpinDaughter(1,EvtSpinType::SCALAR);
70 checkSpinDaughter(2,EvtSpinType::SCALAR);
72 // Check that the B dtr types are K D D:
74 // get the parameters:
76 double phase = getArg(1);
77 _exp = EvtComplex(cos(phase), sin(phase));
80 //------------------------------------------------------------------
81 void EvtBtoKD3P::initProbMax(){
82 setProbMax(1); // this is later changed in decay()
85 //------------------------------------------------------------------
86 void EvtBtoKD3P::decay(EvtParticle *p){
87 // tell the subclass that we decay the daughter:
88 _daugsDecayedByParentModel = true;
90 // the K is the 1st daughter of the B EvtParticle.
91 // The decay mode of the allowed D (the one produced in b->c decay) is 2nd
92 // The decay mode of the suppressed D (the one produced in b->u decay) is 3rd
97 // generate kinematics of daughters (K and D):
98 EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)};
99 p->initializePhaseSpace(2, tempDaug);
101 // Get the D daughter particle and the decay models of the allowed
102 // and suppressed D modes:
103 EvtParticle * theD = p->getDaug(D1IND);
104 EvtPto3P * model1 = (EvtPto3P*)(EvtDecayTable::getInstance()->getDecayFunc(theD));
106 // for the suppressed mode, re-initialize theD as the suppressed D alias:
107 theD->init(getDaug(D2IND), theD->getP4());
108 EvtPto3P * model2 = (EvtPto3P*)(EvtDecayTable::getInstance()->getDecayFunc(theD));
110 // on the first call:
111 if (false == _decayedOnce) {
114 // store the D decay model pointers:
118 // check the decay models of the first 2 daughters and that they
119 // have the same final states:
120 std::string name1=model1->getName();
121 std::string name2=model2->getName();
123 if (name1 != "PTO3P") {
124 report(ERROR,"EvtGen")
125 << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model"
127 << " but found to decay via " << name1.c_str()
128 << " or " << name2.c_str()
129 << ". Will terminate execution!" << endl;
133 EvtId * daugs1 = model1->getDaugs();
134 EvtId * daugs2 = model2->getDaugs();
138 for (d = 0; d < 2; ++d) {
139 if (daugs1[d] != daugs2[d]) {
143 if (false == idMatch) {
144 report(ERROR,"EvtGen")
145 << "D daughters of EvtBtoKD3P decay must decay to the same final state"
147 << " particles in the same order (not CP-conjugate order)," << endl
148 << " but they were found to decay to" << endl;
149 for (d = 0; d < model1->getNDaug(); ++d) {
150 report(ERROR,"") << " " << EvtPDL::name(daugs1[d]).c_str() << " ";
152 report(ERROR,"") << endl;
153 for (d = 0; d < model1->getNDaug(); ++d) {
154 report(ERROR,"") << " " << EvtPDL::name(daugs2[d]).c_str() << " ";
156 report(ERROR,"") << endl << ". Will terminate execution!" << endl;
160 // estimate the probmax. Need to know the probmax's of the 2
162 setProbMax(model1->getProbMax(0)
163 + _r * _r * model2->getProbMax(0)
164 + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0)));
166 } // end of things to do on the first call
168 // make sure the models haven't changed since the first call:
169 if (_model1 != model1 || _model2 != model2) {
170 report(ERROR,"EvtGen")
171 << "D daughters of EvtBtoKD3P decay should have only 1 decay modes, "
173 << " but a new decay mode was found after the first call" << endl
174 << " Will terminate execution!" << endl;
178 // get the cover function for each of the models and add them up.
179 // They are summed with coefficients 1 because we are willing to
180 // take a small inefficiency (~50%) in order to ensure that the
181 // cover function is large enough without getting into complications
182 // associated with the smallness of _r:
183 EvtPdfSum<EvtDalitzPoint> * pc1 = model1->getPC();
184 EvtPdfSum<EvtDalitzPoint> * pc2 = model2->getPC();
185 EvtPdfSum<EvtDalitzPoint> pc;
186 pc.addTerm(1.0, *pc1);
187 pc.addTerm(1.0, *pc2);
189 // from this combined cover function, generate the Dalitz point:
190 EvtDalitzPoint x = pc.randomPoint();
192 // get the aptitude for each of the models on this point and add them up:
193 EvtComplex amp1 = model1->amplNonCP(x);
194 EvtComplex amp2 = model2->amplNonCP(x);
195 EvtComplex amp = amp1 + amp2 * _r * _exp;
197 // get the value of the cover function for this point and set the
198 // relative amplitude for this decay:
200 double comp = sqrt(pc.evaluate (x));
203 // Make the daughters of theD:
204 bool massTreeOK = theD->generateMassTree();
205 if (massTreeOK == false) {return;}
207 // Now generate the p4's of the daughters of theD:
208 std::vector<EvtVector4R> v = model2->initDaughters(x);
210 if(v.size() != theD->getNDaug()) {
211 report(ERROR,"EvtGen")
212 << "Number of daughters " << theD->getNDaug()
213 << " != " << "Momentum vector size " << v.size()
215 << " Terminating execution." << endl;
219 // Apply the new p4's to the daughters:
220 for(unsigned int i=0; i<theD->getNDaug(); ++i){
221 theD->getDaug(i)->init(model2->getDaugs()[i], v[i]);