Line data Source code
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 0 : EvtVectorIsr::~EvtVectorIsr() {}
48 :
49 : std::string EvtVectorIsr::getName(){
50 :
51 0 : return "VECTORISR";
52 : }
53 :
54 : EvtDecayBase* EvtVectorIsr::clone(){
55 :
56 0 : return new EvtVectorIsr;
57 0 : }
58 :
59 : void EvtVectorIsr::init(){
60 :
61 : // check that there are 2 arguments
62 :
63 0 : checkNDaug(2);
64 :
65 0 : checkSpinParent(EvtSpinType::VECTOR);
66 0 : checkSpinDaughter(0,EvtSpinType::VECTOR);
67 0 : checkSpinDaughter(1,EvtSpinType::PHOTON);
68 :
69 0 : int narg = getNArg();
70 0 : if ( narg > 4 ) checkNArg(4);
71 :
72 0 : csfrmn=1.;
73 0 : csbkmn=1.;
74 0 : fmax=1.2;
75 0 : firstorder=false;
76 :
77 0 : if ( narg > 0 ) csfrmn=getArg(0);
78 0 : if ( narg > 1 ) csbkmn=getArg(1);
79 0 : if ( narg > 2 ) fmax=getArg(2);
80 0 : if ( narg > 3 ) firstorder=true;
81 0 : }
82 :
83 :
84 : void EvtVectorIsr::initProbMax(){
85 :
86 0 : noProbMax();
87 0 : }
88 :
89 : void EvtVectorIsr::decay( EvtParticle *p ){
90 :
91 : //the elctron mass
92 0 : double electMass=EvtPDL::getMeanMass(EvtPDL::getId("e-"));
93 :
94 0 : 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 0 : EvtVector4R p4softg1(0.,0.,0.,0.);
101 0 : 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 0 : p->initializePhaseSpace(getNDaug(),getDaugs());
108 0 : phi=p->getDaug(0);
109 0 : 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 0 : double wcm=p->mass();
123 0 : double eb=0.5*wcm;
124 :
125 : //TO guarantee the collinear photons are softer than the ISR photon, require q2m > m*wcm
126 0 : double q2m=phi->mass()*wcm;
127 0 : double f_col(0.);
128 0 : double e01(0.);
129 0 : double e02(0.);
130 0 : double ebeam=eb;
131 : double wcm_new = wcm;
132 0 : 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 0 : if (!firstorder){
140 0 : while (fran > f){
141 0 : m++;
142 :
143 : int n=0;
144 0 : while (f_col == 0.){
145 0 : n++;
146 0 : ckhrad(eb,q2m,e01,e02,f_col);
147 0 : if (n > 10000){
148 0 : report(Severity::Info,"EvtGen") << "EvtVectorIsr is having problems. Called ckhrad 10000 times.\n";
149 0 : assert(0);
150 : }
151 : }
152 :
153 : //Effective beam energy after soft photon emission (neglecting electron mass)
154 0 : ebeam = sqrt(e01*e02);
155 0 : wcm_new = 2*ebeam;
156 0 : s_new = wcm_new*wcm_new;
157 :
158 : //The Vector mass should never be greater than wcm_new
159 0 : if (phi->mass() > wcm_new){
160 0 : report(Severity::Info,"EvtGen") << "EvtVectorIsr finds Vector mass="<<phi->mass()<<" > Weff=" << wcm_new<<". Should not happen\n";
161 0 : 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 0 : if (EvtPDL::getMaxRange(phi->getId()) > 0.) {
168 0 : double x0 = 1 - EvtPDL::getMeanMass(phi->getId())*EvtPDL::getMeanMass(phi->getId())/s_new;
169 :
170 : //L = log(s/(electMass*electMass)
171 0 : double L = 2.*log(wcm_new/electMass);
172 :
173 : // W(x0) is actually 2*alpha/pi times the following
174 0 : 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 0 : cs_Born = W/s_new;
179 0 : }
180 :
181 0 : f = cs_Born*f_col;
182 :
183 : //if fmax was set properly, f should NEVER be larger than fmax
184 0 : if (f > fmax && fmax > 0.){
185 0 : report(Severity::Info,"EvtGen") << "EvtVectorIsr finds a problem with fmax, the maximum weight setting\n"
186 0 : << "fmax is the third decay argument in the .dec file. VectorIsr attempts to set it reasonably if it wasn't provided\n"
187 0 : << "To determine a more appropriate value, build GeneratorQAApp, and set the third argument for this decay <0.\n"
188 0 : << "If you haven't been providing the first 2 arguments, set them to be 1. 1.). The program will report\n"
189 0 : << "the largest weight it finds. You should set fmax to be slightly larger.\n"
190 0 : << "Alternatively try the following values for various vector particles: "
191 0 : << "phi->1.15 J/psi-psi(4415)->0.105\n"
192 0 : << "The current value of f and fmax for " << EvtPDL::name(phi->getId()) << " are " << f << " " << fmax << "\n"
193 0 : << "Will now assert\n";
194 0 : assert(0);
195 : }
196 :
197 :
198 0 : if (fmax > 0.) {
199 0 : fran = fmax*EvtRandom::Flat(0.0,1.0);
200 0 : }
201 :
202 : else {
203 : //determine max weight for this vector particle mass
204 0 : if (f>largest_f) {
205 : largest_f = f;
206 0 : report(Severity::Info,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " "
207 0 : << "vector_mass "
208 0 : << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f
209 0 : << ". f_col cs_B = " << f_col << " " << cs_Born
210 0 : << std::endl;
211 0 : }
212 0 : if (m%10000 == 0) {
213 0 : report(Severity::Info,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " "
214 0 : << "vector_mass "
215 0 : << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f
216 0 : << ". f_col cs_B = " << f_col << " " << cs_Born
217 0 : << std::endl;
218 0 : }
219 :
220 0 : f_col = 0.;
221 : f = 0.;
222 : //determine max weight for this vector particle mass
223 : }
224 :
225 0 : if (m > 100000){
226 :
227 0 : if (fmax > 0.) report(Severity::Info,"EvtGen") << "EvtVectorIsr is having problems. Check the fmax value - the 3rd argument in the .dec file\n"
228 0 : << "Recommended values for various vector particles: "
229 0 : << "phi->1.15 J/psi-psi(4415)->0.105 "
230 0 : << "Upsilon(1S,2S,3S)->0.14\n";
231 0 : 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 0 : if (firstorder){
247 : bet_l = 0.;
248 : gam_l = 1.;
249 : betgam_l = 0.;
250 0 : csfrmn_new = csfrmn;
251 0 : csbkmn_new = csbkmn;
252 0 : } else {
253 0 : double xx = e02/e01;
254 0 : double sq_xx = sqrt(xx);
255 0 : bet_l = (1.-xx)/(1.+xx);
256 0 : gam_l = (1.+xx)/(2.*sq_xx);
257 0 : 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 0 : csfrmn_new=(csfrmn - bet_l)/(1. - bet_l*csfrmn);
261 0 : 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 0 : 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 0 : double beta=electMass/ebeam; //electMass/Ebeam = 1/gamma
277 0 : beta=sqrt(1. - beta*beta); //sqrt (1 - (1/gamma)**2)
278 :
279 0 : double ymax=log((1.+beta*csfrmn_new)/(1.-beta*csfrmn_new));
280 0 : 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 0 : double y=(ymax-ymin)*EvtRandom::Flat(0.0,1.0) + ymin;
284 0 : double cs=exp(y);
285 0 : cs=(cs - 1.)/(cs + 1.)/beta;
286 0 : double sn=sqrt(1-cs*cs);
287 :
288 0 : double fi=EvtRandom::Flat(EvtConst::twoPi);
289 :
290 : //four-vector for the phi
291 0 : double phi_p0 = sqrt(phi->mass()*phi->mass()+pg*pg);
292 0 : double phi_p3 = -pg*cs;
293 :
294 :
295 : //boost back to frame before colinear radiation.
296 0 : EvtVector4R p4phi(gam_l*phi_p0 + betgam_l*phi_p3,
297 0 : -pg*sn*cos(fi),
298 0 : -pg*sn*sin(fi),
299 0 : betgam_l*phi_p0 + gam_l*phi_p3);
300 :
301 : double isr_p0 = pg;
302 0 : double isr_p3 = -phi_p3;
303 0 : EvtVector4R p4gamma(gam_l*isr_p0 + betgam_l*isr_p3,
304 0 : -p4phi.get(1),
305 0 : -p4phi.get(2),
306 0 : betgam_l*isr_p0 + gam_l*isr_p3);
307 :
308 :
309 : //four-vectors of the collinear photons
310 0 : if (!firstorder) {
311 0 : p4softg1.set(0, eb-e02); p4softg1.set(3, e02-eb);
312 0 : p4softg2.set(0, eb-e01); p4softg2.set(3, eb-e01);
313 0 : }
314 :
315 : //save momenta for particles
316 0 : phi->init( getDaug(0),p4phi);
317 0 : gamma->init( getDaug(1),p4gamma);
318 :
319 :
320 : //add the two colinear photons as vphoton daughters
321 0 : EvtPhotonParticle *softg1=new EvtPhotonParticle;;
322 0 : EvtPhotonParticle *softg2=new EvtPhotonParticle;;
323 0 : softg1->init(gammaId,p4softg1);
324 0 : softg2->init(gammaId,p4softg2);
325 0 : softg1->addDaug(p);
326 0 : 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 0 : EvtVector4C phi0=phi->epsParent(0);
331 0 : EvtVector4C phi1=phi->epsParent(1);
332 0 : EvtVector4C phi2=phi->epsParent(2);
333 :
334 : //get polarization vector for a photon in its parents restframe.
335 0 : EvtVector4C gamma0=gamma->epsParentPhoton(0);
336 0 : EvtVector4C gamma1=gamma->epsParentPhoton(1);
337 :
338 0 : EvtComplex r1p=phi0*gamma0;
339 0 : EvtComplex r2p=phi1*gamma0;
340 0 : EvtComplex r3p=phi2*gamma0;
341 :
342 :
343 0 : EvtComplex r1m=phi0*gamma1;
344 0 : EvtComplex r2m=phi1*gamma1;
345 0 : EvtComplex r3m=phi2*gamma1;
346 :
347 0 : EvtComplex rho33=r3p*conj(r3p)+r3m*conj(r3m);
348 0 : EvtComplex rho22=r2p*conj(r2p)+r2m*conj(r2m);
349 0 : EvtComplex rho11=r1p*conj(r1p)+r1m*conj(r1m);
350 :
351 0 : EvtComplex rho13=r3p*conj(r1p)+r3m*conj(r1m);
352 0 : EvtComplex rho12=r2p*conj(r1p)+r2m*conj(r1m);
353 0 : EvtComplex rho23=r3p*conj(r2p)+r3m*conj(r2m);
354 :
355 0 : EvtComplex rho31=conj(rho13);
356 0 : EvtComplex rho32=conj(rho23);
357 0 : EvtComplex rho21=conj(rho12);
358 :
359 :
360 0 : EvtSpinDensity rho;
361 0 : rho.setDim(3);
362 :
363 0 : rho.set(0,0,rho11);
364 0 : rho.set(0,1,rho12);
365 0 : rho.set(0,2,rho13);
366 0 : rho.set(1,0,rho21);
367 0 : rho.set(1,1,rho22);
368 0 : rho.set(1,2,rho23);
369 0 : rho.set(2,0,rho31);
370 0 : rho.set(2,1,rho32);
371 0 : rho.set(2,2,rho33);
372 :
373 0 : setDaughterSpinDensity(0);
374 0 : phi->setSpinDensityForward(rho);
375 :
376 : return ;
377 0 : }
378 :
379 : double EvtVectorIsr::ckhrad1(double xx, double a, double b){
380 : //port of AfkQed/ckhrad.F function ckhrad1
381 0 : double yy = xx*xx;
382 0 : double zz = 1. - 2*xx + yy;
383 0 : 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 0 : const double adp = 1. / 137.0359895 / EvtConst::pi;
389 0 : const double pi2 = EvtConst::pi*EvtConst::pi;
390 : // const double dme = 0.00051099906;
391 0 : const double dme = EvtPDL::getMeanMass(EvtPDL::getId("e-"));
392 :
393 0 : double r1=EvtRandom::Flat();//Generates Flat from 0 - 1
394 0 : double r2=EvtRandom::Flat();
395 :
396 0 : double sss = 4.*e_beam*e_beam;
397 0 : double biglog = log(sss/(dme*dme));
398 0 : double beta = 2.*adp*(biglog - 1.);
399 : double betae_lab = beta;
400 0 : double p3 = adp*(pi2/3. - 0.5);
401 0 : double p12 = adp*adp * (11./8. - 2.*pi2/3.);
402 0 : double coefener = 1. + 0.75*betae_lab + p3;
403 0 : double coef1 = coefener + 0.125*pi2*beta*beta;
404 0 : double coef2 = p12* biglog*biglog;
405 0 : double facts = coef1 + coef2;
406 :
407 : double y1_min = 0;
408 0 : double e1min = 0.25 * q2_min/e_beam;
409 0 : double y1_max = pow( 1. - e1min/e_beam, 0.5*beta );
410 0 : double y1 = y1_min +r1 *(y1_max - y1_min);
411 0 : e01 = e_beam *(1. - pow(y1, 2./beta) );
412 :
413 : double y2_min = 0.;
414 0 : double e2min = 0.25 * q2_min/e01;
415 0 : double y2_max = pow( 1. - e2min/e_beam, 0.5*beta);
416 0 : double y2 = y2_min +r2 *(y2_max - y2_min);
417 0 : e02 = e_beam *(1. - pow(y2, 2./beta) );
418 :
419 :
420 0 : double xx1 = e01/e_beam;
421 0 : double xx2 = e02/e_beam;
422 :
423 0 : f = y1_max * y2_max * ckhrad1(xx1,biglog,betae_lab) * ckhrad1(xx2,biglog,betae_lab) * facts;
424 :
425 : return;
426 0 : }
427 :
|
