Line data Source code
1 : //-----------------------------------------------------------------------
2 : // File and Version Information:
3 : // $Id: EvtBtoKD3P.cpp,v 1.2 2009-04-02 15:22:28 robbep Exp $
4 : //
5 : // Environment:
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.
9 : //
10 : // Copyright Information:
11 : // Copyright (C) 2003, Colorado State University
12 : //
13 : // Module creator:
14 : // Abi soffer, CSU, 2003
15 : //-----------------------------------------------------------------------
16 : #include "EvtGenBase/EvtPatches.hh"
17 :
18 : // Decay model that does the decay B+->D0K, D0->3 psudoscalars
19 :
20 : #include <assert.h>
21 :
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"
31 : using std::endl;
32 :
33 : //------------------------------------------------------------------
34 0 : EvtBtoKD3P::EvtBtoKD3P() :
35 0 : _model1(0),
36 0 : _model2(0),
37 0 : _decayedOnce(false)
38 0 : {
39 0 : }
40 :
41 : //------------------------------------------------------------------
42 0 : EvtBtoKD3P::EvtBtoKD3P(const EvtBtoKD3P & other) :
43 0 : EvtDecayAmp( other ){
44 0 : }
45 :
46 : //------------------------------------------------------------------
47 0 : EvtBtoKD3P::~EvtBtoKD3P(){
48 0 : }
49 :
50 : //------------------------------------------------------------------
51 : EvtDecayBase * EvtBtoKD3P::clone(){
52 0 : return new EvtBtoKD3P();
53 0 : }
54 :
55 : //------------------------------------------------------------------
56 : std::string EvtBtoKD3P::getName(){
57 0 : return "BTOKD3P";
58 : }
59 :
60 : //------------------------------------------------------------------
61 : void EvtBtoKD3P::init(){
62 0 : checkNArg(2); // r, phase
63 0 : checkNDaug(3); // K, D0(allowed), D0(suppressed).
64 : // The last two daughters are really one particle
65 :
66 : // check that the mother and all daughters are scalars:
67 0 : checkSpinParent ( EvtSpinType::SCALAR);
68 0 : checkSpinDaughter(0,EvtSpinType::SCALAR);
69 0 : checkSpinDaughter(1,EvtSpinType::SCALAR);
70 0 : checkSpinDaughter(2,EvtSpinType::SCALAR);
71 :
72 : // Check that the B dtr types are K D D:
73 :
74 : // get the parameters:
75 0 : _r = getArg(0);
76 0 : double phase = getArg(1);
77 0 : _exp = EvtComplex(cos(phase), sin(phase));
78 0 : }
79 :
80 : //------------------------------------------------------------------
81 : void EvtBtoKD3P::initProbMax(){
82 0 : setProbMax(1); // this is later changed in decay()
83 0 : }
84 :
85 : //------------------------------------------------------------------
86 : void EvtBtoKD3P::decay(EvtParticle *p){
87 : // tell the subclass that we decay the daughter:
88 0 : _daugsDecayedByParentModel = true;
89 :
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
93 : const int KIND = 0;
94 : const int D1IND = 1;
95 : const int D2IND = 2;
96 :
97 : // generate kinematics of daughters (K and D):
98 0 : EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)};
99 0 : p->initializePhaseSpace(2, tempDaug);
100 :
101 : // Get the D daughter particle and the decay models of the allowed
102 : // and suppressed D modes:
103 0 : EvtParticle * theD = p->getDaug(D1IND);
104 0 : EvtPto3P * model1 = (EvtPto3P*)(EvtDecayTable::getInstance()->getDecayFunc(theD));
105 :
106 : // for the suppressed mode, re-initialize theD as the suppressed D alias:
107 0 : theD->init(getDaug(D2IND), theD->getP4());
108 0 : EvtPto3P * model2 = (EvtPto3P*)(EvtDecayTable::getInstance()->getDecayFunc(theD));
109 :
110 : // on the first call:
111 0 : if (false == _decayedOnce) {
112 0 : _decayedOnce = true;
113 :
114 : // store the D decay model pointers:
115 0 : _model1 = model1;
116 0 : _model2 = model2;
117 :
118 : // check the decay models of the first 2 daughters and that they
119 : // have the same final states:
120 0 : std::string name1=model1->getName();
121 0 : std::string name2=model2->getName();
122 :
123 0 : if (name1 != "PTO3P") {
124 0 : report(Severity::Error,"EvtGen")
125 0 : << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model"
126 0 : << endl
127 0 : << " but found to decay via " << name1.c_str()
128 0 : << " or " << name2.c_str()
129 0 : << ". Will terminate execution!" << endl;
130 0 : assert(0);
131 : }
132 :
133 0 : EvtId * daugs1 = model1->getDaugs();
134 0 : EvtId * daugs2 = model2->getDaugs();
135 :
136 : bool idMatch = true;
137 : int d;
138 0 : for (d = 0; d < 2; ++d) {
139 0 : if (daugs1[d] != daugs2[d]) {
140 : idMatch = false;
141 0 : }
142 : }
143 0 : if (false == idMatch) {
144 0 : report(Severity::Error,"EvtGen")
145 0 : << "D daughters of EvtBtoKD3P decay must decay to the same final state"
146 0 : << endl
147 0 : << " particles in the same order (not CP-conjugate order)," << endl
148 0 : << " but they were found to decay to" << endl;
149 0 : for (d = 0; d < model1->getNDaug(); ++d) {
150 0 : report(Severity::Error,"") << " " << EvtPDL::name(daugs1[d]).c_str() << " ";
151 : }
152 0 : report(Severity::Error,"") << endl;
153 0 : for (d = 0; d < model1->getNDaug(); ++d) {
154 0 : report(Severity::Error,"") << " " << EvtPDL::name(daugs2[d]).c_str() << " ";
155 : }
156 0 : report(Severity::Error,"") << endl << ". Will terminate execution!" << endl;
157 0 : assert(0);
158 : }
159 :
160 : // estimate the probmax. Need to know the probmax's of the 2
161 : // models for this:
162 0 : setProbMax(model1->getProbMax(0)
163 0 : + _r * _r * model2->getProbMax(0)
164 0 : + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0)));
165 :
166 0 : } // end of things to do on the first call
167 :
168 : // make sure the models haven't changed since the first call:
169 0 : if (_model1 != model1 || _model2 != model2) {
170 0 : report(Severity::Error,"EvtGen")
171 0 : << "D daughters of EvtBtoKD3P decay should have only 1 decay modes, "
172 0 : << endl
173 0 : << " but a new decay mode was found after the first call" << endl
174 0 : << " Will terminate execution!" << endl;
175 0 : assert(0);
176 : }
177 :
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 0 : EvtPdfSum<EvtDalitzPoint> * pc1 = model1->getPC();
184 0 : EvtPdfSum<EvtDalitzPoint> * pc2 = model2->getPC();
185 0 : EvtPdfSum<EvtDalitzPoint> pc;
186 0 : pc.addTerm(1.0, *pc1);
187 0 : pc.addTerm(1.0, *pc2);
188 :
189 : // from this combined cover function, generate the Dalitz point:
190 0 : EvtDalitzPoint x = pc.randomPoint();
191 :
192 : // get the aptitude for each of the models on this point and add them up:
193 0 : EvtComplex amp1 = model1->amplNonCP(x);
194 0 : EvtComplex amp2 = model2->amplNonCP(x);
195 0 : EvtComplex amp = amp1 + amp2 * _r * _exp;
196 :
197 : // get the value of the cover function for this point and set the
198 : // relative amplitude for this decay:
199 :
200 0 : double comp = sqrt(pc.evaluate (x));
201 0 : vertex (amp/comp);
202 :
203 : // Make the daughters of theD:
204 0 : bool massTreeOK = theD->generateMassTree();
205 0 : if (massTreeOK == false) {return;}
206 :
207 : // Now generate the p4's of the daughters of theD:
208 0 : std::vector<EvtVector4R> v = model2->initDaughters(x);
209 :
210 0 : if(v.size() != theD->getNDaug()) {
211 0 : report(Severity::Error,"EvtGen")
212 0 : << "Number of daughters " << theD->getNDaug()
213 0 : << " != " << "Momentum vector size " << v.size()
214 0 : << endl
215 0 : << " Terminating execution." << endl;
216 0 : assert(0);
217 : }
218 :
219 : // Apply the new p4's to the daughters:
220 0 : for(unsigned int i=0; i<theD->getNDaug(); ++i){
221 0 : theD->getDaug(i)->init(model2->getDaugs()[i], v[i]);
222 : }
223 0 : }
224 :
|
