Line data Source code
1 : //-----------------------------------------------------------------------
2 : // File and Version Information:
3 : // $Id: EvtPto3PAmpFactory.cpp,v 1.3 2009-03-16 15:44:04 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) 1998 Caltech, UCSB
12 : //
13 : // Module creator:
14 : // Alexei Dvoretskii, Caltech, 2001-2002.
15 : //-----------------------------------------------------------------------
16 : #include "EvtGenBase/EvtPatches.hh"
17 :
18 : // AmpFactory for building a P -> 3P decay
19 : // (pseudoscalar to three pseudoscalars)
20 :
21 : #include <assert.h>
22 : #include <math.h>
23 : #include <stdio.h>
24 : #include <stdlib.h>
25 :
26 : #include "EvtGenBase/EvtId.hh"
27 : #include "EvtGenBase/EvtPDL.hh"
28 : #include "EvtGenBase/EvtConst.hh"
29 : #include "EvtGenBase/EvtComplex.hh"
30 : #include "EvtGenBase/EvtCyclic3.hh"
31 : #include "EvtGenBase/EvtSpinType.hh"
32 : #include "EvtGenBase/EvtPto3PAmp.hh"
33 : #include "EvtGenBase/EvtNonresonantAmp.hh"
34 : #include "EvtGenBase/EvtFlatAmp.hh"
35 : #include "EvtGenBase/EvtLASSAmp.hh"
36 : #include "EvtGenBase/EvtPto3PAmpFactory.hh"
37 : #include "EvtGenBase/EvtPropBreitWigner.hh"
38 : #include "EvtGenBase/EvtPropFlatte.hh"
39 : #include "EvtGenBase/EvtPropBreitWignerRel.hh"
40 : #include "EvtGenBase/EvtDalitzResPdf.hh"
41 : #include "EvtGenBase/EvtDalitzFlatPdf.hh"
42 :
43 : using namespace EvtCyclic3;
44 : #include <iostream>
45 :
46 : void EvtPto3PAmpFactory::processAmp(EvtComplex c, std::vector<std::string> vv, bool conj)
47 : {
48 0 : if(_verbose) {
49 :
50 0 : printf("Make %samplitude\n",conj ? "CP conjugate" : "");
51 : unsigned i;
52 0 : for(i=0;i<vv.size();i++) printf("%s\n",vv[i].c_str());
53 0 : printf("\n");
54 0 : }
55 :
56 : EvtAmplitude<EvtDalitzPoint>* amp = 0;
57 : EvtPdf<EvtDalitzPoint>* pdf = 0;
58 0 : std::string name;
59 : Pair pairRes=AB;
60 :
61 : size_t i;
62 : /*
63 : Experimental amplitudes
64 : */
65 0 : if(vv[0] == "PHASESPACE") {
66 :
67 0 : pdf = new EvtDalitzFlatPdf(_dp);
68 0 : amp = new EvtFlatAmp<EvtDalitzPoint>();
69 0 : name = "NR";
70 : }
71 0 : else if (!vv[0].find("NONRES")) {
72 : double alpha=0;
73 : EvtPto3PAmp::NumType typeNRes=EvtPto3PAmp::NONRES;
74 0 : if (vv[0]=="NONRES_LIN") {
75 : typeNRes=EvtPto3PAmp::NONRES_LIN;
76 0 : pairRes=strToPair(vv[1].c_str());
77 0 : }
78 0 : else if (vv[0]=="NONRES_EXP") {
79 : typeNRes=EvtPto3PAmp::NONRES_EXP;
80 0 : pairRes = strToPair(vv[1].c_str());
81 0 : alpha = strtod(vv[2].c_str(),0);
82 : }
83 0 : else assert(0);
84 0 : pdf = new EvtDalitzFlatPdf(_dp);
85 0 : amp = new EvtNonresonantAmp( &_dp, typeNRes, pairRes, alpha);
86 0 : }
87 0 : else if (vv[0]=="LASS" || vv[0]=="LASS_ELASTIC" || vv[0]=="LASS_RESONANT") {
88 0 : pairRes = strToPair(vv[1].c_str());
89 0 : double m0 = strtod(vv[2].c_str(),0);
90 0 : double g0 = strtod(vv[3].c_str(),0);
91 0 : double a = strtod(vv[4].c_str(),0);
92 0 : double r = strtod(vv[5].c_str(),0);
93 0 : double cutoff = strtod(vv[6].c_str(),0);
94 0 : pdf = new EvtDalitzResPdf(_dp,m0,g0,pairRes);
95 0 : amp = new EvtLASSAmp( &_dp, pairRes, m0, g0, a, r, cutoff, vv[0]);
96 0 : }
97 :
98 : /*
99 : Resonant amplitudes
100 : */
101 0 : else if(vv[0] == "RESONANCE") {
102 : EvtPto3PAmp* partAmp = 0;
103 :
104 : // RESONANCE stanza
105 :
106 0 : pairRes = strToPair(vv[1].c_str());
107 : EvtSpinType::spintype spinR = EvtSpinType::SCALAR;
108 : double mR, gR;
109 0 : name = vv[2];
110 0 : EvtId resId = EvtPDL::getId(vv[2]);
111 0 : if(_verbose) printf("Particles %s form %sresonance %s\n",
112 0 : vv[1].c_str(),vv[2].c_str(), conj ? "(conj) " : "");
113 :
114 : // If no valid particle name is given, assume that
115 : // it is the spin, the mass and the width of the particle.
116 :
117 0 : if(resId.getId() == -1) {
118 :
119 0 : switch(atoi(vv[2].c_str())) {
120 :
121 0 : case 0: { spinR = EvtSpinType::SCALAR; break; }
122 0 : case 1: { spinR = EvtSpinType::VECTOR; break; }
123 0 : case 2: { spinR = EvtSpinType::TENSOR; break; }
124 0 : case 3: { spinR = EvtSpinType::SPIN3; break; }
125 0 : case 4: { spinR = EvtSpinType::SPIN4; break; }
126 0 : default: { assert(0); break; }
127 : }
128 :
129 0 : mR = strtod(vv[3].c_str(),0);
130 0 : gR = strtod(vv[4].c_str(),0);
131 : i = 4;
132 0 : }
133 : else {
134 :
135 : // For a valid particle get spin, mass and width
136 :
137 0 : spinR = EvtPDL::getSpinType(resId);
138 0 : mR = EvtPDL::getMeanMass(resId);
139 0 : gR = EvtPDL::getWidth(resId);
140 : i = 2;
141 :
142 : // It's possible to specify mass and width of a particle
143 : // explicitly
144 :
145 0 : if(vv[3] != "ANGULAR") {
146 :
147 0 : if(_verbose)
148 0 : printf("Setting m(%s)=%s g(%s)=%s\n",
149 0 : vv[2].c_str(),vv[3].c_str(),vv[2].c_str(),vv[4].c_str());
150 :
151 0 : mR = strtod(vv[3].c_str(),0);
152 0 : gR = strtod(vv[4].c_str(),0);
153 : i = 4;
154 0 : }
155 : }
156 :
157 : // ANGULAR stanza
158 :
159 0 : if(vv[++i] != "ANGULAR") {
160 :
161 0 : printf("%s instead of ANGULAR\n",vv[i].c_str());
162 0 : exit(0);
163 : }
164 0 : Pair pairAng = strToPair(vv[++i].c_str());
165 0 : if(_verbose) printf("Angle is measured between particles %s\n",vv[i].c_str());
166 :
167 : // TYPE stanza
168 :
169 0 : std::string typeName = vv[++i];
170 0 : assert(typeName == "TYPE");
171 0 : std::string type = vv[++i];
172 0 : if(_verbose) printf("Propagator type %s\n",vv[i].c_str());
173 :
174 0 : if(type == "NBW") {
175 :
176 0 : EvtPropBreitWigner prop(mR,gR);
177 0 : partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::NBW);
178 0 : }
179 0 : else if(type == "RBW_ZEMACH") {
180 :
181 0 : EvtPropBreitWignerRel prop(mR,gR);
182 0 : partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_ZEMACH);
183 0 : }
184 0 : else if(type == "RBW_KUEHN") {
185 :
186 0 : EvtPropBreitWignerRel prop(mR,gR);
187 0 : partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_KUEHN);
188 0 : }
189 0 : else if(type == "RBW_CLEO") {
190 :
191 0 : EvtPropBreitWignerRel prop(mR,gR);
192 0 : partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::RBW_CLEO);
193 0 : }
194 0 : else if(type == "FLATTE") {
195 :
196 0 : double m1a = _dp.m( first(pairRes) );
197 0 : double m1b = _dp.m( second(pairRes) );
198 : // 2nd channel
199 0 : double g2 = strtod(vv[++i].c_str(),0);
200 0 : double m2a = strtod(vv[++i].c_str(),0);
201 0 : double m2b = strtod(vv[++i].c_str(),0);
202 0 : EvtPropFlatte prop( mR, gR, m1a, m1b, g2, m2a, m2b );
203 0 : partAmp = new EvtPto3PAmp(_dp,pairAng,pairRes,spinR,prop,EvtPto3PAmp::FLATTE);
204 0 : }
205 0 : else assert(0);
206 :
207 : // Optional DVFF, BVFF stanzas
208 :
209 0 : if(i < vv.size() - 1) {
210 0 : if(vv[i+1] == "DVFF") {
211 : i++;
212 0 : if(vv[++i] == "BLATTWEISSKOPF") {
213 :
214 0 : double R = strtod(vv[++i].c_str(),0);
215 0 : partAmp->set_fd(R);
216 : }
217 0 : else assert(0);
218 0 : }
219 : }
220 :
221 0 : if(i < vv.size() - 1) {
222 0 : if(vv[i+1] == "BVFF") {
223 : i++;
224 0 : if(vv[++i] == "BLATTWEISSKOPF") {
225 :
226 0 : if(_verbose) printf("BVFF=%s\n",vv[i].c_str());
227 0 : double R = strtod(vv[++i].c_str(),0);
228 0 : partAmp->set_fb(R);
229 : }
230 0 : else assert(0);
231 0 : }
232 : }
233 :
234 : const int minwidths=5;
235 : //Optional resonance minimum and maximum
236 0 : if(i < vv.size() - 1) {
237 0 : if(vv[i+1] == "CUTOFF") {
238 : i++;
239 0 : if(vv[i+1] == "MIN") {
240 : i++;
241 0 : double min = strtod(vv[++i].c_str(),0);
242 0 : if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<" "<<minwidths<<std::endl;
243 : //ensure against cutting off too close to the resonance
244 0 : assert( min<(mR-minwidths*gR) );
245 0 : partAmp->setmin(min);
246 0 : }
247 0 : else if (vv[i+1] == "MAX") {
248 : i++;
249 0 : double max = strtod(vv[++i].c_str(),0);
250 0 : if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<" "<<minwidths<<std::endl;
251 : //ensure against cutting off too close to the resonance
252 0 : assert( max>(mR+minwidths*gR) );
253 0 : partAmp->setmax(max);
254 : }
255 0 : else assert(0);
256 : }
257 : }
258 :
259 : //2nd iteration in case min and max are both specified
260 0 : if(i < vv.size() - 1) {
261 0 : if(vv[i+1] == "CUTOFF") {
262 : i++;
263 0 : if(vv[i+1] == "MIN") {
264 : i++;
265 0 : double min = strtod(vv[++i].c_str(),0);
266 0 : if(_verbose) std::cout<<"CUTOFF MIN = "<<min<<std::endl;
267 : //ensure against cutting off too close to the resonance
268 0 : assert( min<(mR-minwidths*gR) );
269 0 : partAmp->setmin(min);
270 0 : }
271 0 : else if (vv[i+1] == "MAX") {
272 : i++;
273 0 : double max = strtod(vv[++i].c_str(),0);
274 0 : if(_verbose) std::cout<<"CUTOFF MAX = "<<max<<std::endl;
275 : //ensure against cutting off too close to the resonance
276 0 : assert( max>(mR+minwidths*gR) );
277 0 : partAmp->setmax(max);
278 : }
279 0 : else assert(0);
280 : }
281 : }
282 :
283 :
284 : i++;
285 :
286 0 : pdf = new EvtDalitzResPdf(_dp,mR,gR,pairRes);
287 0 : amp = partAmp;
288 0 : }
289 :
290 0 : assert(amp);
291 0 : assert(pdf);
292 :
293 0 : if(!conj) {
294 0 : _amp->addOwnedTerm(c,amp);
295 : }
296 : else {
297 0 : _ampConj->addOwnedTerm(c,amp);
298 : }
299 :
300 0 : double scale = matchIsobarCoef(_amp, pdf, pairRes);
301 0 : _pc->addOwnedTerm(abs2(c)*scale,pdf);
302 :
303 0 : _names.push_back(name);
304 0 : }
305 :
306 : double EvtPto3PAmpFactory::matchIsobarCoef(EvtAmplitude<EvtDalitzPoint>* amp,
307 : EvtPdf<EvtDalitzPoint>* pdf,
308 : EvtCyclic3::Pair ipair) {
309 :
310 : // account for differences in the definition of amplitudes by matching
311 : // Integral( c'*pdf ) = Integral( c*|A|^2 )
312 : // to improve generation efficiency ...
313 :
314 0 : double Ipdf = pdf->compute_integral(10000).value();
315 : double Iamp2 = 0;
316 :
317 :
318 0 : EvtCyclic3::Pair jpair = EvtCyclic3::next(ipair);
319 0 : EvtCyclic3::Pair kpair = EvtCyclic3::next(jpair);
320 :
321 : // Trapezoidal integral
322 : int N=10000;
323 :
324 0 : double di = (_dp.qAbsMax(ipair) - _dp.qAbsMin(ipair))/((double) N);
325 :
326 0 : double siMin = _dp.qAbsMin(ipair);
327 :
328 0 : double s[3]; // playing with fire
329 0 : for(int i=1; i<N; i++) {
330 :
331 0 : s[ipair] = siMin + di*i;
332 0 : s[jpair] = _dp.q(jpair, 0.9999, ipair, s[ipair]);
333 0 : s[kpair] = _dp.bigM()*_dp.bigM() - s[ipair] - s[jpair]
334 0 : + _dp.mA()*_dp.mA() + _dp.mB()*_dp.mB() + _dp.mC()*_dp.mC();
335 :
336 0 : EvtDalitzPoint point( _dp.mA(), _dp.mB(), _dp.mC(),
337 0 : s[EvtCyclic3::AB], s[EvtCyclic3::BC], s[EvtCyclic3::CA]);
338 :
339 0 : if (!point.isValid()) continue;
340 :
341 0 : double p = point.p(other(ipair), ipair);
342 0 : double q = point.p(first(ipair), ipair);
343 :
344 0 : double itg = abs2( amp->evaluate(point) )*di*4*q*p;
345 0 : Iamp2 += itg;
346 :
347 0 : }
348 0 : if (_verbose) std::cout << "integral = " << Iamp2 << " pdf="<<Ipdf << std::endl;
349 :
350 0 : assert(Ipdf>0 && Iamp2>0);
351 :
352 0 : return Iamp2/Ipdf;
353 0 : }
|
