]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //----------------------------------------------------------------------- |
2 | // File and Version Information: | |
0ca57c2f | 3 | // $Id: EvtBtoKD3P.cpp,v 1.2 2009-04-02 15:22:28 robbep Exp $ |
da0e9ce3 | 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 | //------------------------------------------------------------------ | |
0ca57c2f | 42 | EvtBtoKD3P::EvtBtoKD3P(const EvtBtoKD3P & other) : |
43 | EvtDecayAmp( other ){ | |
da0e9ce3 | 44 | } |
45 | ||
46 | //------------------------------------------------------------------ | |
47 | EvtBtoKD3P::~EvtBtoKD3P(){ | |
48 | } | |
49 | ||
50 | //------------------------------------------------------------------ | |
51 | EvtDecayBase * EvtBtoKD3P::clone(){ | |
52 | return new EvtBtoKD3P(); | |
53 | } | |
54 | ||
55 | //------------------------------------------------------------------ | |
56 | std::string EvtBtoKD3P::getName(){ | |
57 | return "BTOKD3P"; | |
58 | } | |
59 | ||
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 | |
65 | ||
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); | |
71 | ||
72 | // Check that the B dtr types are K D D: | |
73 | ||
74 | // get the parameters: | |
75 | _r = getArg(0); | |
76 | double phase = getArg(1); | |
77 | _exp = EvtComplex(cos(phase), sin(phase)); | |
78 | } | |
79 | ||
80 | //------------------------------------------------------------------ | |
81 | void EvtBtoKD3P::initProbMax(){ | |
82 | setProbMax(1); // this is later changed in decay() | |
83 | } | |
84 | ||
85 | //------------------------------------------------------------------ | |
86 | void EvtBtoKD3P::decay(EvtParticle *p){ | |
87 | // tell the subclass that we decay the daughter: | |
88 | _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 | EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)}; | |
99 | p->initializePhaseSpace(2, tempDaug); | |
100 | ||
101 | // Get the D daughter particle and the decay models of the allowed | |
102 | // and suppressed D modes: | |
103 | EvtParticle * theD = p->getDaug(D1IND); | |
0ca57c2f | 104 | EvtPto3P * model1 = (EvtPto3P*)(EvtDecayTable::getInstance()->getDecayFunc(theD)); |
da0e9ce3 | 105 | |
106 | // for the suppressed mode, re-initialize theD as the suppressed D alias: | |
107 | theD->init(getDaug(D2IND), theD->getP4()); | |
0ca57c2f | 108 | EvtPto3P * model2 = (EvtPto3P*)(EvtDecayTable::getInstance()->getDecayFunc(theD)); |
da0e9ce3 | 109 | |
110 | // on the first call: | |
111 | if (false == _decayedOnce) { | |
112 | _decayedOnce = true; | |
113 | ||
114 | // store the D decay model pointers: | |
115 | _model1 = model1; | |
116 | _model2 = model2; | |
117 | ||
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(); | |
122 | ||
123 | if (name1 != "PTO3P") { | |
124 | report(ERROR,"EvtGen") | |
125 | << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model" | |
126 | << endl | |
127 | << " but found to decay via " << name1.c_str() | |
128 | << " or " << name2.c_str() | |
129 | << ". Will terminate execution!" << endl; | |
130 | assert(0); | |
131 | } | |
132 | ||
133 | EvtId * daugs1 = model1->getDaugs(); | |
134 | EvtId * daugs2 = model2->getDaugs(); | |
135 | ||
136 | bool idMatch = true; | |
137 | int d; | |
138 | for (d = 0; d < 2; ++d) { | |
139 | if (daugs1[d] != daugs2[d]) { | |
140 | idMatch = false; | |
141 | } | |
142 | } | |
143 | if (false == idMatch) { | |
144 | report(ERROR,"EvtGen") | |
145 | << "D daughters of EvtBtoKD3P decay must decay to the same final state" | |
146 | << endl | |
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() << " "; | |
151 | } | |
152 | report(ERROR,"") << endl; | |
153 | for (d = 0; d < model1->getNDaug(); ++d) { | |
154 | report(ERROR,"") << " " << EvtPDL::name(daugs2[d]).c_str() << " "; | |
155 | } | |
156 | report(ERROR,"") << endl << ". Will terminate execution!" << endl; | |
157 | assert(0); | |
158 | } | |
159 | ||
160 | // estimate the probmax. Need to know the probmax's of the 2 | |
161 | // models for this: | |
162 | setProbMax(model1->getProbMax(0) | |
163 | + _r * _r * model2->getProbMax(0) | |
164 | + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0))); | |
165 | ||
166 | } // end of things to do on the first call | |
167 | ||
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, " | |
172 | << endl | |
173 | << " but a new decay mode was found after the first call" << endl | |
174 | << " Will terminate execution!" << endl; | |
175 | 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 | 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); | |
188 | ||
189 | // from this combined cover function, generate the Dalitz point: | |
190 | EvtDalitzPoint x = pc.randomPoint(); | |
191 | ||
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; | |
196 | ||
197 | // get the value of the cover function for this point and set the | |
198 | // relative amplitude for this decay: | |
199 | ||
200 | double comp = sqrt(pc.evaluate (x)); | |
201 | vertex (amp/comp); | |
202 | ||
203 | // Make the daughters of theD: | |
0ca57c2f | 204 | bool massTreeOK = theD->generateMassTree(); |
205 | if (massTreeOK == false) {return;} | |
da0e9ce3 | 206 | |
207 | // Now generate the p4's of the daughters of theD: | |
208 | std::vector<EvtVector4R> v = model2->initDaughters(x); | |
209 | ||
210 | if(v.size() != theD->getNDaug()) { | |
211 | report(ERROR,"EvtGen") | |
212 | << "Number of daughters " << theD->getNDaug() | |
213 | << " != " << "Momentum vector size " << v.size() | |
214 | << endl | |
215 | << " Terminating execution." << endl; | |
216 | assert(0); | |
217 | } | |
218 | ||
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]); | |
222 | } | |
223 | } | |
224 |