]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //----------------------------------------------------------------------- |
2 | // File and Version Information: | |
3 | // $Id: EvtBtoKD3P.cc,v 1.7 2009/02/19 03:23:29 ryd 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 | EvtBtoKD3P::EvtBtoKD3P() : | |
35 | _model1(0), | |
36 | _model2(0), | |
37 | _decayedOnce(false) | |
38 | { | |
39 | } | |
40 | ||
41 | //------------------------------------------------------------------ | |
42 | EvtBtoKD3P::EvtBtoKD3P(const EvtBtoKD3P & other){ | |
43 | } | |
44 | ||
45 | //------------------------------------------------------------------ | |
46 | EvtBtoKD3P::~EvtBtoKD3P(){ | |
47 | } | |
48 | ||
49 | //------------------------------------------------------------------ | |
50 | EvtDecayBase * EvtBtoKD3P::clone(){ | |
51 | return new EvtBtoKD3P(); | |
52 | } | |
53 | ||
54 | //------------------------------------------------------------------ | |
55 | std::string EvtBtoKD3P::getName(){ | |
56 | return "BTOKD3P"; | |
57 | } | |
58 | ||
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 | |
64 | ||
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); | |
70 | ||
71 | // Check that the B dtr types are K D D: | |
72 | ||
73 | // get the parameters: | |
74 | _r = getArg(0); | |
75 | double phase = getArg(1); | |
76 | _exp = EvtComplex(cos(phase), sin(phase)); | |
77 | } | |
78 | ||
79 | //------------------------------------------------------------------ | |
80 | void EvtBtoKD3P::initProbMax(){ | |
81 | setProbMax(1); // this is later changed in decay() | |
82 | } | |
83 | ||
84 | //------------------------------------------------------------------ | |
85 | void EvtBtoKD3P::decay(EvtParticle *p){ | |
86 | // tell the subclass that we decay the daughter: | |
87 | _daugsDecayedByParentModel = true; | |
88 | ||
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 | |
92 | const int KIND = 0; | |
93 | const int D1IND = 1; | |
94 | const int D2IND = 2; | |
95 | ||
96 | // generate kinematics of daughters (K and D): | |
97 | EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)}; | |
98 | p->initializePhaseSpace(2, tempDaug); | |
99 | ||
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)); | |
104 | ||
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)); | |
108 | ||
109 | // on the first call: | |
110 | if (false == _decayedOnce) { | |
111 | _decayedOnce = true; | |
112 | ||
113 | // store the D decay model pointers: | |
114 | _model1 = model1; | |
115 | _model2 = model2; | |
116 | ||
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(); | |
121 | ||
122 | if (name1 != "PTO3P") { | |
123 | report(ERROR,"EvtGen") | |
124 | << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model" | |
125 | << endl | |
126 | << " but found to decay via " << name1.c_str() | |
127 | << " or " << name2.c_str() | |
128 | << ". Will terminate execution!" << endl; | |
129 | assert(0); | |
130 | } | |
131 | ||
132 | EvtId * daugs1 = model1->getDaugs(); | |
133 | EvtId * daugs2 = model2->getDaugs(); | |
134 | ||
135 | bool idMatch = true; | |
136 | int d; | |
137 | for (d = 0; d < 2; ++d) { | |
138 | if (daugs1[d] != daugs2[d]) { | |
139 | idMatch = false; | |
140 | } | |
141 | } | |
142 | if (false == idMatch) { | |
143 | report(ERROR,"EvtGen") | |
144 | << "D daughters of EvtBtoKD3P decay must decay to the same final state" | |
145 | << endl | |
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() << " "; | |
150 | } | |
151 | report(ERROR,"") << endl; | |
152 | for (d = 0; d < model1->getNDaug(); ++d) { | |
153 | report(ERROR,"") << " " << EvtPDL::name(daugs2[d]).c_str() << " "; | |
154 | } | |
155 | report(ERROR,"") << endl << ". Will terminate execution!" << endl; | |
156 | assert(0); | |
157 | } | |
158 | ||
159 | // estimate the probmax. Need to know the probmax's of the 2 | |
160 | // models for this: | |
161 | setProbMax(model1->getProbMax(0) | |
162 | + _r * _r * model2->getProbMax(0) | |
163 | + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0))); | |
164 | ||
165 | } // end of things to do on the first call | |
166 | ||
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, " | |
171 | << endl | |
172 | << " but a new decay mode was found after the first call" << endl | |
173 | << " Will terminate execution!" << endl; | |
174 | assert(0); | |
175 | } | |
176 | ||
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); | |
187 | ||
188 | // from this combined cover function, generate the Dalitz point: | |
189 | EvtDalitzPoint x = pc.randomPoint(); | |
190 | ||
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; | |
195 | ||
196 | // get the value of the cover function for this point and set the | |
197 | // relative amplitude for this decay: | |
198 | ||
199 | double comp = sqrt(pc.evaluate (x)); | |
200 | vertex (amp/comp); | |
201 | ||
202 | // Make the daughters of theD: | |
203 | theD->generateMassTree(); | |
204 | ||
205 | // Now generate the p4's of the daughters of theD: | |
206 | std::vector<EvtVector4R> v = model2->initDaughters(x); | |
207 | ||
208 | if(v.size() != theD->getNDaug()) { | |
209 | report(ERROR,"EvtGen") | |
210 | << "Number of daughters " << theD->getNDaug() | |
211 | << " != " << "Momentum vector size " << v.size() | |
212 | << endl | |
213 | << " Terminating execution." << endl; | |
214 | assert(0); | |
215 | } | |
216 | ||
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]); | |
220 | } | |
221 | } | |
222 |