]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //-------------------------------------------------------------------------- |
2 | // | |
3 | // Environment: | |
4 | // This software is part of the EvtGen package developed jointly | |
5 | // for the BaBar and CLEO collaborations. If you use all or part | |
6 | // of it, please give an appropriate acknowledgement. | |
7 | // | |
8 | // Copyright Information: See EvtGen/COPYRIGHT | |
9 | // Copyright (C) 1998 Caltech, UCSB | |
10 | // | |
11 | // Module: EvtVectorIsr.cc | |
12 | // | |
13 | // Description: | |
14 | // This is a special decay model to generate e+e- -> phi gamma + soft gammas | |
15 | // using soft collinear ISR calculation from AfkQed | |
16 | // This is implemented as a decay of the VPHO. | |
17 | // | |
18 | // Modification history: | |
19 | // | |
20 | // Joe Izen Oct, 2005 Soft Colinear Photons (secondary ISR) ported from AfkQed | |
21 | // Joe Izen Dec 16, 2002 Fix cos_theta distribution - prevents boom at cos_theta=+/-1 | |
22 | // RYD/Adriano June 16, 1998 Module created | |
23 | // | |
24 | //------------------------------------------------------------------------ | |
25 | // | |
26 | #include "EvtGenBase/EvtPatches.hh" | |
27 | #include <stdlib.h> | |
28 | ||
29 | #include <math.h> | |
30 | #include <iostream> | |
31 | #include <iomanip> | |
32 | #include <sstream> | |
33 | ||
34 | ||
35 | #include "EvtGenBase/EvtParticle.hh" | |
36 | #include "EvtGenBase/EvtPhotonParticle.hh" | |
37 | #include "EvtGenBase/EvtRandom.hh" | |
38 | #include "EvtGenBase/EvtPDL.hh" | |
39 | #include "EvtGenBase/EvtAbsLineShape.hh" | |
40 | #include "EvtGenModels/EvtVectorIsr.hh" | |
41 | #include "EvtGenBase/EvtReport.hh" | |
42 | #include "EvtGenBase/EvtConst.hh" | |
43 | #include "EvtGenBase/EvtAbsLineShape.hh" | |
44 | #include <string> | |
45 | #include "EvtGenBase/EvtVector4C.hh" | |
46 | ||
47 | EvtVectorIsr::~EvtVectorIsr() {} | |
48 | ||
49 | std::string EvtVectorIsr::getName(){ | |
50 | ||
51 | return "VECTORISR"; | |
52 | } | |
53 | ||
54 | EvtDecayBase* EvtVectorIsr::clone(){ | |
55 | ||
56 | return new EvtVectorIsr; | |
57 | } | |
58 | ||
59 | void EvtVectorIsr::init(){ | |
60 | ||
61 | // check that there are 2 arguments | |
62 | ||
63 | checkNDaug(2); | |
64 | ||
65 | checkSpinParent(EvtSpinType::VECTOR); | |
66 | checkSpinDaughter(0,EvtSpinType::VECTOR); | |
67 | checkSpinDaughter(1,EvtSpinType::PHOTON); | |
68 | ||
69 | int narg = getNArg(); | |
70 | if ( narg > 4 ) checkNArg(4); | |
71 | ||
72 | csfrmn=1.; | |
73 | csbkmn=1.; | |
74 | fmax=1.2; | |
75 | firstorder=false; | |
76 | ||
77 | if ( narg > 0 ) csfrmn=getArg(0); | |
78 | if ( narg > 1 ) csbkmn=getArg(1); | |
79 | if ( narg > 2 ) fmax=getArg(2); | |
80 | if ( narg > 3 ) firstorder=true; | |
81 | } | |
82 | ||
83 | ||
84 | void EvtVectorIsr::initProbMax(){ | |
85 | ||
86 | noProbMax(); | |
87 | } | |
88 | ||
89 | void EvtVectorIsr::decay( EvtParticle *p ){ | |
90 | ||
91 | //the elctron mass | |
92 | double electMass=EvtPDL::getMeanMass(EvtPDL::getId("e-")); | |
93 | ||
94 | static EvtId gammaId=EvtPDL::getId("gamma"); | |
95 | ||
96 | EvtParticle *phi; | |
97 | EvtParticle *gamma; | |
98 | ||
99 | //4-mom of the two colinear photons to the decay of the vphoton | |
100 | EvtVector4R p4softg1(0.,0.,0.,0.); | |
101 | EvtVector4R p4softg2(0.,0.,0.,0.); | |
102 | ||
103 | ||
104 | //get pointers to the daughters set | |
105 | //get masses/initial phase space - will overwrite the | |
106 | //p4s below to get the kinematic distributions correct | |
107 | p->initializePhaseSpace(getNDaug(),getDaugs()); | |
108 | phi=p->getDaug(0); | |
109 | gamma=p->getDaug(1); | |
110 | ||
111 | //Generate soft colinear photons and the electron and positron energies after emission. | |
112 | //based on method of AfkQed and notes of Vladimir Druzhinin. | |
113 | // | |
114 | //function ckhrad(eb,q2m,r1,r2,e01,e02,f_col) | |
115 | //eb: energy of incoming electrons in CM frame | |
116 | //q2m: minimum invariant mass of the virtual photon after soft colinear photon emission | |
117 | //returned arguments | |
118 | //e01,e02: energies of e+ and e- after soft colinear photon emission | |
119 | //fcol: weighting factor for Born cross section for use in an accept/reject test. | |
120 | ||
121 | ||
122 | double wcm=p->mass(); | |
123 | double eb=0.5*wcm; | |
124 | ||
125 | //TO guarantee the collinear photons are softer than the ISR photon, require q2m > m*wcm | |
126 | double q2m=phi->mass()*wcm; | |
127 | double f_col(0.); | |
128 | double e01(0.); | |
129 | double e02(0.); | |
130 | double ebeam=eb; | |
131 | double wcm_new = wcm; | |
132 | double s_new = wcm*wcm; | |
133 | ||
134 | double fran = 1.; | |
135 | double f = 0; | |
136 | int m = 0; | |
137 | double largest_f=0;//only used when determining max weight for this vector particle mass | |
138 | ||
139 | if (!firstorder){ | |
140 | while (fran > f){ | |
141 | m++; | |
142 | ||
143 | int n=0; | |
144 | while (f_col == 0.){ | |
145 | n++; | |
146 | ckhrad(eb,q2m,e01,e02,f_col); | |
147 | if (n > 10000){ | |
148 | report(INFO,"EvtGen") << "EvtVectorIsr is having problems. Called ckhrad 10000 times.\n"; | |
149 | assert(0); | |
150 | } | |
151 | } | |
152 | ||
153 | //Effective beam energy after soft photon emission (neglecting electron mass) | |
154 | ebeam = sqrt(e01*e02); | |
155 | wcm_new = 2*ebeam; | |
156 | s_new = wcm_new*wcm_new; | |
157 | ||
158 | //The Vector mass should never be greater than wcm_new | |
159 | if (phi->mass() > wcm_new){ | |
160 | report(INFO,"EvtGen") << "EvtVectorIsr finds Vector mass="<<phi->mass()<<" > Weff=" << wcm_new<<". Should not happen\n"; | |
161 | assert(0); | |
162 | } | |
163 | ||
164 | //Determine Born cross section @ wcm_new for e+e- -> gamma V. We aren't interested in the absolute normalization | |
165 | //Just the functional dependence. Assuming a narrow resonance when determining cs_Born | |
166 | double cs_Born = 1.; | |
167 | if (EvtPDL::getMaxRange(phi->getId()) > 0.) { | |
168 | double x0 = 1 - EvtPDL::getMeanMass(phi->getId())*EvtPDL::getMeanMass(phi->getId())/s_new; | |
169 | ||
170 | //L = log(s/(electMass*electMass) | |
171 | double L = 2.*log(wcm_new/electMass); | |
172 | ||
173 | // W(x0) is actually 2*alpha/pi times the following | |
174 | double W = (L-1.)*(1. - x0 +0.5*x0*x0); | |
175 | ||
176 | //Born cross section is actually 12*pi*pi*Gammaee/EvtPDL::getMeanMass(phi->getId()) times the following | |
177 | //(we'd need the full W(x0) as well) | |
178 | cs_Born = W/s_new; | |
179 | } | |
180 | ||
181 | f = cs_Born*f_col; | |
182 | ||
183 | //if fmax was set properly, f should NEVER be larger than fmax | |
184 | if (f > fmax && fmax > 0.){ | |
185 | report(INFO,"EvtGen") << "EvtVectorIsr finds a problem with fmax, the maximum weight setting\n" | |
186 | << "fmax is the third decay argument in the .dec file. VectorIsr attempts to set it reasonably if it wasn't provided\n" | |
187 | << "To determine a more appropriate value, build GeneratorQAApp, and set the third argument for this decay <0.\n" | |
188 | << "If you haven't been providing the first 2 arguments, set them to be 1. 1.). The program will report\n" | |
189 | << "the largest weight it finds. You should set fmax to be slightly larger.\n" | |
190 | << "Alternatively try the following values for various vector particles: " | |
191 | << "phi->1.15 J/psi-psi(4415)->0.105\n" | |
192 | << "The current value of f and fmax for " << EvtPDL::name(phi->getId()) << " are " << f << " " << fmax << "\n" | |
193 | << "Will now assert\n"; | |
194 | assert(0); | |
195 | } | |
196 | ||
197 | ||
198 | if (fmax > 0.) { | |
199 | fran = fmax*EvtRandom::Flat(0.0,1.0); | |
200 | } | |
201 | ||
202 | else { | |
203 | //determine max weight for this vector particle mass | |
204 | if (f>largest_f) { | |
205 | largest_f = f; | |
206 | report(INFO,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " " | |
207 | << "vector_mass " | |
208 | << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f | |
209 | << ". f_col cs_B = " << f_col << " " << cs_Born | |
210 | << std::endl; | |
211 | } | |
212 | if (m%10000 == 0) { | |
213 | report(INFO,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " " | |
214 | << "vector_mass " | |
215 | << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f | |
216 | << ". f_col cs_B = " << f_col << " " << cs_Born | |
217 | << std::endl; | |
218 | } | |
219 | ||
220 | f_col = 0.; | |
221 | f = 0.; | |
222 | //determine max weight for this vector particle mass | |
223 | } | |
224 | ||
225 | if (m > 100000){ | |
226 | ||
227 | if (fmax > 0.) report(INFO,"EvtGen") << "EvtVectorIsr is having problems. Check the fmax value - the 3rd argument in the .dec file\n" | |
228 | << "Recommended values for various vector particles: " | |
229 | << "phi->1.15 J/psi-psi(4415)->0.105 " | |
230 | << "Upsilon(1S,2S,3S)->0.14\n"; | |
231 | assert(0); | |
232 | } | |
233 | }//while (fran > f) | |
234 | ||
235 | }//if (firstorder) | |
236 | ||
237 | //Compute parameters for boost to/from the system after colinear radiation | |
238 | ||
239 | double bet_l; | |
240 | double gam_l; | |
241 | double betgam_l; | |
242 | ||
243 | double csfrmn_new; | |
244 | double csbkmn_new; | |
245 | ||
246 | if (firstorder){ | |
247 | bet_l = 0.; | |
248 | gam_l = 1.; | |
249 | betgam_l = 0.; | |
250 | csfrmn_new = csfrmn; | |
251 | csbkmn_new = csbkmn; | |
252 | } else { | |
253 | double xx = e02/e01; | |
254 | double sq_xx = sqrt(xx); | |
255 | bet_l = (1.-xx)/(1.+xx); | |
256 | gam_l = (1.+xx)/(2.*sq_xx); | |
257 | betgam_l = (1.-xx)/(2.*sq_xx); | |
258 | ||
259 | //Boost photon cos_theta limits in lab to limits in the system after colinear rad | |
260 | csfrmn_new=(csfrmn - bet_l)/(1. - bet_l*csfrmn); | |
261 | csbkmn_new=(csbkmn - bet_l)/(1. - bet_l*csbkmn); | |
262 | } | |
263 | ||
264 | // //generate kinematics according to Bonneau-Martin article | |
265 | // //Nucl. Phys. B27 (1971) 381-397 | |
266 | ||
267 | // For backward compatibility with .dec files before SP5, the backward cos limit for | |
268 | //the ISR photon is actually given as *minus* the actual limit. Sorry, this wouldn't be | |
269 | //my choice. -Joe | |
270 | ||
271 | //gamma momentum in the vpho restframe *after* soft colinear radiation | |
272 | double pg = (s_new - phi->mass()*phi->mass())/(2.*wcm_new); | |
273 | ||
274 | ||
275 | //calculate the beta of incoming electrons after colinear rad in the frame where e= and e- have equal momentum | |
276 | double beta=electMass/ebeam; //electMass/Ebeam = 1/gamma | |
277 | beta=sqrt(1. - beta*beta); //sqrt (1 - (1/gamma)**2) | |
278 | ||
279 | double ymax=log((1.+beta*csfrmn_new)/(1.-beta*csfrmn_new)); | |
280 | double ymin=log((1.-beta*csbkmn_new)/(1.+beta*csbkmn_new)); | |
281 | ||
282 | // photon theta distributed as 2*beta/(1-beta**2*cos(theta)**2) | |
283 | double y=(ymax-ymin)*EvtRandom::Flat(0.0,1.0) + ymin; | |
284 | double cs=exp(y); | |
285 | cs=(cs - 1.)/(cs + 1.)/beta; | |
286 | double sn=sqrt(1-cs*cs); | |
287 | ||
288 | double fi=EvtRandom::Flat(EvtConst::twoPi); | |
289 | ||
290 | //four-vector for the phi | |
291 | double phi_p0 = sqrt(phi->mass()*phi->mass()+pg*pg); | |
292 | double phi_p3 = -pg*cs; | |
293 | ||
294 | ||
295 | //boost back to frame before colinear radiation. | |
296 | EvtVector4R p4phi(gam_l*phi_p0 + betgam_l*phi_p3, | |
297 | -pg*sn*cos(fi), | |
298 | -pg*sn*sin(fi), | |
299 | betgam_l*phi_p0 + gam_l*phi_p3); | |
300 | ||
301 | double isr_p0 = pg; | |
302 | double isr_p3 = -phi_p3; | |
303 | EvtVector4R p4gamma(gam_l*isr_p0 + betgam_l*isr_p3, | |
304 | -p4phi.get(1), | |
305 | -p4phi.get(2), | |
306 | betgam_l*isr_p0 + gam_l*isr_p3); | |
307 | ||
308 | ||
309 | //four-vectors of the collinear photons | |
310 | if (!firstorder) { | |
311 | p4softg1.set(0, eb-e02); p4softg1.set(3, e02-eb); | |
312 | p4softg2.set(0, eb-e01); p4softg2.set(3, eb-e01); | |
313 | } | |
314 | ||
315 | //save momenta for particles | |
316 | phi->init( getDaug(0),p4phi); | |
317 | gamma->init( getDaug(1),p4gamma); | |
318 | ||
319 | ||
320 | //add the two colinear photons as vphoton daughters | |
321 | EvtPhotonParticle *softg1=new EvtPhotonParticle;; | |
322 | EvtPhotonParticle *softg2=new EvtPhotonParticle;; | |
323 | softg1->init(gammaId,p4softg1); | |
324 | softg2->init(gammaId,p4softg2); | |
325 | softg1->addDaug(p); | |
326 | softg2->addDaug(p); | |
327 | ||
328 | //try setting the spin density matrix of the phi | |
329 | //get polarization vector for phi in its parents restframe. | |
330 | EvtVector4C phi0=phi->epsParent(0); | |
331 | EvtVector4C phi1=phi->epsParent(1); | |
332 | EvtVector4C phi2=phi->epsParent(2); | |
333 | ||
334 | //get polarization vector for a photon in its parents restframe. | |
335 | EvtVector4C gamma0=gamma->epsParentPhoton(0); | |
336 | EvtVector4C gamma1=gamma->epsParentPhoton(1); | |
337 | ||
338 | EvtComplex r1p=phi0*gamma0; | |
339 | EvtComplex r2p=phi1*gamma0; | |
340 | EvtComplex r3p=phi2*gamma0; | |
341 | ||
342 | ||
343 | EvtComplex r1m=phi0*gamma1; | |
344 | EvtComplex r2m=phi1*gamma1; | |
345 | EvtComplex r3m=phi2*gamma1; | |
346 | ||
347 | EvtComplex rho33=r3p*conj(r3p)+r3m*conj(r3m); | |
348 | EvtComplex rho22=r2p*conj(r2p)+r2m*conj(r2m); | |
349 | EvtComplex rho11=r1p*conj(r1p)+r1m*conj(r1m); | |
350 | ||
351 | EvtComplex rho13=r3p*conj(r1p)+r3m*conj(r1m); | |
352 | EvtComplex rho12=r2p*conj(r1p)+r2m*conj(r1m); | |
353 | EvtComplex rho23=r3p*conj(r2p)+r3m*conj(r2m); | |
354 | ||
355 | EvtComplex rho31=conj(rho13); | |
356 | EvtComplex rho32=conj(rho23); | |
357 | EvtComplex rho21=conj(rho12); | |
358 | ||
359 | ||
360 | EvtSpinDensity rho; | |
361 | rho.setDim(3); | |
362 | ||
363 | rho.set(0,0,rho11); | |
364 | rho.set(0,1,rho12); | |
365 | rho.set(0,2,rho13); | |
366 | rho.set(1,0,rho21); | |
367 | rho.set(1,1,rho22); | |
368 | rho.set(1,2,rho23); | |
369 | rho.set(2,0,rho31); | |
370 | rho.set(2,1,rho32); | |
371 | rho.set(2,2,rho33); | |
372 | ||
373 | setDaughterSpinDensity(0); | |
374 | phi->setSpinDensityForward(rho); | |
375 | ||
376 | return ; | |
377 | } | |
378 | ||
379 | double EvtVectorIsr::ckhrad1(double xx, double a, double b){ | |
380 | //port of AfkQed/ckhrad.F function ckhrad1 | |
381 | double yy = xx*xx; | |
382 | double zz = 1. - 2*xx + yy; | |
383 | return 0.5* (1. + yy + zz/(a-1.) + 0.25*b*( -0.5*(1. + 3*yy)*log(xx)) - zz ); | |
384 | } | |
385 | ||
386 | void EvtVectorIsr::ckhrad(const double& e_beam,const double& q2_min,double& e01,double& e02,double& f){ | |
387 | //port of AfkQed/ckhrad.F subroutine ckhrad | |
388 | const double adp = 1. / 137.0359895 / EvtConst::pi; | |
389 | const double pi2 = EvtConst::pi*EvtConst::pi; | |
390 | // const double dme = 0.00051099906; | |
391 | const double dme = EvtPDL::getMeanMass(EvtPDL::getId("e-")); | |
392 | ||
393 | double r1=EvtRandom::Flat();//Generates Flat from 0 - 1 | |
394 | double r2=EvtRandom::Flat(); | |
395 | ||
396 | double sss = 4.*e_beam*e_beam; | |
397 | double biglog = log(sss/(dme*dme)); | |
398 | double beta = 2.*adp*(biglog - 1.); | |
399 | double betae_lab = beta; | |
400 | double p3 = adp*(pi2/3. - 0.5); | |
401 | double p12 = adp*adp * (11./8. - 2.*pi2/3.); | |
402 | double coefener = 1. + 0.75*betae_lab + p3; | |
403 | double coef1 = coefener + 0.125*pi2*beta*beta; | |
404 | double coef2 = p12* biglog*biglog; | |
405 | double facts = coef1 + coef2; | |
406 | ||
407 | double y1_min = 0; | |
408 | double e1min = 0.25 * q2_min/e_beam; | |
409 | double y1_max = pow( 1. - e1min/e_beam, 0.5*beta ); | |
410 | double y1 = y1_min +r1 *(y1_max - y1_min); | |
411 | e01 = e_beam *(1. - pow(y1, 2./beta) ); | |
412 | ||
413 | double y2_min = 0.; | |
414 | double e2min = 0.25 * q2_min/e01; | |
415 | double y2_max = pow( 1. - e2min/e_beam, 0.5*beta); | |
416 | double y2 = y2_min +r2 *(y2_max - y2_min); | |
417 | e02 = e_beam *(1. - pow(y2, 2./beta) ); | |
418 | ||
419 | ||
420 | double xx1 = e01/e_beam; | |
421 | double xx2 = e02/e_beam; | |
422 | ||
423 | f = y1_max * y2_max * ckhrad1(xx1,biglog,betae_lab) * ckhrad1(xx2,biglog,betae_lab) * facts; | |
424 | ||
425 | return; | |
426 | } | |
427 |