1 //--------------------------------------------------------------------------
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.
8 // Copyright Information: See EvtGen/COPYRIGHT
9 // Copyright (C) 1998 Caltech, UCSB
11 // Module: EvtVectorIsr.cc
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.
18 // Modification history:
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
24 //------------------------------------------------------------------------
26 #include "EvtGenBase/EvtPatches.hh"
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"
45 #include "EvtGenBase/EvtVector4C.hh"
47 EvtVectorIsr::~EvtVectorIsr() {}
49 std::string EvtVectorIsr::getName(){
54 EvtDecayBase* EvtVectorIsr::clone(){
56 return new EvtVectorIsr;
59 void EvtVectorIsr::init(){
61 // check that there are 2 arguments
65 checkSpinParent(EvtSpinType::VECTOR);
66 checkSpinDaughter(0,EvtSpinType::VECTOR);
67 checkSpinDaughter(1,EvtSpinType::PHOTON);
70 if ( narg > 4 ) checkNArg(4);
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;
84 void EvtVectorIsr::initProbMax(){
89 void EvtVectorIsr::decay( EvtParticle *p ){
92 double electMass=EvtPDL::getMeanMass(EvtPDL::getId("e-"));
94 static EvtId gammaId=EvtPDL::getId("gamma");
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.);
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());
111 //Generate soft colinear photons and the electron and positron energies after emission.
112 //based on method of AfkQed and notes of Vladimir Druzhinin.
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
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.
122 double wcm=p->mass();
125 //TO guarantee the collinear photons are softer than the ISR photon, require q2m > m*wcm
126 double q2m=phi->mass()*wcm;
131 double wcm_new = wcm;
132 double s_new = wcm*wcm;
137 double largest_f=0;//only used when determining max weight for this vector particle mass
146 ckhrad(eb,q2m,e01,e02,f_col);
148 report(INFO,"EvtGen") << "EvtVectorIsr is having problems. Called ckhrad 10000 times.\n";
153 //Effective beam energy after soft photon emission (neglecting electron mass)
154 ebeam = sqrt(e01*e02);
156 s_new = wcm_new*wcm_new;
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";
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
167 if (EvtPDL::getMaxRange(phi->getId()) > 0.) {
168 double x0 = 1 - EvtPDL::getMeanMass(phi->getId())*EvtPDL::getMeanMass(phi->getId())/s_new;
170 //L = log(s/(electMass*electMass)
171 double L = 2.*log(wcm_new/electMass);
173 // W(x0) is actually 2*alpha/pi times the following
174 double W = (L-1.)*(1. - x0 +0.5*x0*x0);
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)
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";
199 fran = fmax*EvtRandom::Flat(0.0,1.0);
203 //determine max weight for this vector particle mass
206 report(INFO,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " "
208 << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f
209 << ". f_col cs_B = " << f_col << " " << cs_Born
213 report(INFO,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " "
215 << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f
216 << ". f_col cs_B = " << f_col << " " << cs_Born
222 //determine max weight for this vector particle mass
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";
237 //Compute parameters for boost to/from the system after colinear radiation
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);
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);
264 // //generate kinematics according to Bonneau-Martin article
265 // //Nucl. Phys. B27 (1971) 381-397
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
271 //gamma momentum in the vpho restframe *after* soft colinear radiation
272 double pg = (s_new - phi->mass()*phi->mass())/(2.*wcm_new);
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)
279 double ymax=log((1.+beta*csfrmn_new)/(1.-beta*csfrmn_new));
280 double ymin=log((1.-beta*csbkmn_new)/(1.+beta*csbkmn_new));
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;
285 cs=(cs - 1.)/(cs + 1.)/beta;
286 double sn=sqrt(1-cs*cs);
288 double fi=EvtRandom::Flat(EvtConst::twoPi);
290 //four-vector for the phi
291 double phi_p0 = sqrt(phi->mass()*phi->mass()+pg*pg);
292 double phi_p3 = -pg*cs;
295 //boost back to frame before colinear radiation.
296 EvtVector4R p4phi(gam_l*phi_p0 + betgam_l*phi_p3,
299 betgam_l*phi_p0 + gam_l*phi_p3);
302 double isr_p3 = -phi_p3;
303 EvtVector4R p4gamma(gam_l*isr_p0 + betgam_l*isr_p3,
306 betgam_l*isr_p0 + gam_l*isr_p3);
309 //four-vectors of the collinear photons
311 p4softg1.set(0, eb-e02); p4softg1.set(3, e02-eb);
312 p4softg2.set(0, eb-e01); p4softg2.set(3, eb-e01);
315 //save momenta for particles
316 phi->init( getDaug(0),p4phi);
317 gamma->init( getDaug(1),p4gamma);
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);
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);
334 //get polarization vector for a photon in its parents restframe.
335 EvtVector4C gamma0=gamma->epsParentPhoton(0);
336 EvtVector4C gamma1=gamma->epsParentPhoton(1);
338 EvtComplex r1p=phi0*gamma0;
339 EvtComplex r2p=phi1*gamma0;
340 EvtComplex r3p=phi2*gamma0;
343 EvtComplex r1m=phi0*gamma1;
344 EvtComplex r2m=phi1*gamma1;
345 EvtComplex r3m=phi2*gamma1;
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);
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);
355 EvtComplex rho31=conj(rho13);
356 EvtComplex rho32=conj(rho23);
357 EvtComplex rho21=conj(rho12);
373 setDaughterSpinDensity(0);
374 phi->setSpinDensityForward(rho);
379 double EvtVectorIsr::ckhrad1(double xx, double a, double b){
380 //port of AfkQed/ckhrad.F function ckhrad1
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 );
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-"));
393 double r1=EvtRandom::Flat();//Generates Flat from 0 - 1
394 double r2=EvtRandom::Flat();
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;
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) );
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) );
420 double xx1 = e01/e_beam;
421 double xx2 = e02/e_beam;
423 f = y1_max * y2_max * ckhrad1(xx1,biglog,betae_lab) * ckhrad1(xx2,biglog,betae_lab) * facts;