--- /dev/null
+//--------------------------------------------------------------------------
+//
+// Environment:
+// This software is part of the EvtGen package developed jointly
+// for the BaBar and CLEO collaborations. If you use all or part
+// of it, please give an appropriate acknowledgement.
+//
+// Copyright Information: See EvtGen/COPYRIGHT
+// Copyright (C) 1998 Caltech, UCSB
+//
+// Module: EvtVectorIsr.cc
+//
+// Description:
+// This is a special decay model to generate e+e- -> phi gamma + soft gammas
+// using soft collinear ISR calculation from AfkQed
+// This is implemented as a decay of the VPHO.
+//
+// Modification history:
+//
+// Joe Izen Oct, 2005 Soft Colinear Photons (secondary ISR) ported from AfkQed
+// Joe Izen Dec 16, 2002 Fix cos_theta distribution - prevents boom at cos_theta=+/-1
+// RYD/Adriano June 16, 1998 Module created
+//
+//------------------------------------------------------------------------
+//
+#include "EvtGenBase/EvtPatches.hh"
+#include <stdlib.h>
+
+#include <math.h>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+
+#include "EvtGenBase/EvtParticle.hh"
+#include "EvtGenBase/EvtPhotonParticle.hh"
+#include "EvtGenBase/EvtRandom.hh"
+#include "EvtGenBase/EvtPDL.hh"
+#include "EvtGenBase/EvtAbsLineShape.hh"
+#include "EvtGenModels/EvtVectorIsr.hh"
+#include "EvtGenBase/EvtReport.hh"
+#include "EvtGenBase/EvtConst.hh"
+#include "EvtGenBase/EvtAbsLineShape.hh"
+#include <string>
+#include "EvtGenBase/EvtVector4C.hh"
+
+EvtVectorIsr::~EvtVectorIsr() {}
+
+std::string EvtVectorIsr::getName(){
+
+ return "VECTORISR";
+}
+
+EvtDecayBase* EvtVectorIsr::clone(){
+
+ return new EvtVectorIsr;
+}
+
+void EvtVectorIsr::init(){
+
+ // check that there are 2 arguments
+
+ checkNDaug(2);
+
+ checkSpinParent(EvtSpinType::VECTOR);
+ checkSpinDaughter(0,EvtSpinType::VECTOR);
+ checkSpinDaughter(1,EvtSpinType::PHOTON);
+
+ int narg = getNArg();
+ if ( narg > 4 ) checkNArg(4);
+
+ csfrmn=1.;
+ csbkmn=1.;
+ fmax=1.2;
+ firstorder=false;
+
+ if ( narg > 0 ) csfrmn=getArg(0);
+ if ( narg > 1 ) csbkmn=getArg(1);
+ if ( narg > 2 ) fmax=getArg(2);
+ if ( narg > 3 ) firstorder=true;
+}
+
+
+void EvtVectorIsr::initProbMax(){
+
+ noProbMax();
+}
+
+void EvtVectorIsr::decay( EvtParticle *p ){
+
+ //the elctron mass
+ double electMass=EvtPDL::getMeanMass(EvtPDL::getId("e-"));
+
+ static EvtId gammaId=EvtPDL::getId("gamma");
+
+ EvtParticle *phi;
+ EvtParticle *gamma;
+
+ //4-mom of the two colinear photons to the decay of the vphoton
+ EvtVector4R p4softg1(0.,0.,0.,0.);
+ EvtVector4R p4softg2(0.,0.,0.,0.);
+
+
+ //get pointers to the daughters set
+ //get masses/initial phase space - will overwrite the
+ //p4s below to get the kinematic distributions correct
+ p->initializePhaseSpace(getNDaug(),getDaugs());
+ phi=p->getDaug(0);
+ gamma=p->getDaug(1);
+
+ //Generate soft colinear photons and the electron and positron energies after emission.
+ //based on method of AfkQed and notes of Vladimir Druzhinin.
+ //
+ //function ckhrad(eb,q2m,r1,r2,e01,e02,f_col)
+ //eb: energy of incoming electrons in CM frame
+ //q2m: minimum invariant mass of the virtual photon after soft colinear photon emission
+ //returned arguments
+ //e01,e02: energies of e+ and e- after soft colinear photon emission
+ //fcol: weighting factor for Born cross section for use in an accept/reject test.
+
+
+ double wcm=p->mass();
+ double eb=0.5*wcm;
+
+ //TO guarantee the collinear photons are softer than the ISR photon, require q2m > m*wcm
+ double q2m=phi->mass()*wcm;
+ double f_col(0.);
+ double e01(0.);
+ double e02(0.);
+ double ebeam=eb;
+ double wcm_new = wcm;
+ double s_new = wcm*wcm;
+
+ double fran = 1.;
+ double f = 0;
+ int m = 0;
+ double largest_f=0;//only used when determining max weight for this vector particle mass
+
+ if (!firstorder){
+ while (fran > f){
+ m++;
+
+ int n=0;
+ while (f_col == 0.){
+ n++;
+ ckhrad(eb,q2m,e01,e02,f_col);
+ if (n > 10000){
+ report(INFO,"EvtGen") << "EvtVectorIsr is having problems. Called ckhrad 10000 times.\n";
+ assert(0);
+ }
+ }
+
+ //Effective beam energy after soft photon emission (neglecting electron mass)
+ ebeam = sqrt(e01*e02);
+ wcm_new = 2*ebeam;
+ s_new = wcm_new*wcm_new;
+
+ //The Vector mass should never be greater than wcm_new
+ if (phi->mass() > wcm_new){
+ report(INFO,"EvtGen") << "EvtVectorIsr finds Vector mass="<<phi->mass()<<" > Weff=" << wcm_new<<". Should not happen\n";
+ assert(0);
+ }
+
+ //Determine Born cross section @ wcm_new for e+e- -> gamma V. We aren't interested in the absolute normalization
+ //Just the functional dependence. Assuming a narrow resonance when determining cs_Born
+ double cs_Born = 1.;
+ if (EvtPDL::getMaxRange(phi->getId()) > 0.) {
+ double x0 = 1 - EvtPDL::getMeanMass(phi->getId())*EvtPDL::getMeanMass(phi->getId())/s_new;
+
+ //L = log(s/(electMass*electMass)
+ double L = 2.*log(wcm_new/electMass);
+
+ // W(x0) is actually 2*alpha/pi times the following
+ double W = (L-1.)*(1. - x0 +0.5*x0*x0);
+
+ //Born cross section is actually 12*pi*pi*Gammaee/EvtPDL::getMeanMass(phi->getId()) times the following
+ //(we'd need the full W(x0) as well)
+ cs_Born = W/s_new;
+ }
+
+ f = cs_Born*f_col;
+
+ //if fmax was set properly, f should NEVER be larger than fmax
+ if (f > fmax && fmax > 0.){
+ report(INFO,"EvtGen") << "EvtVectorIsr finds a problem with fmax, the maximum weight setting\n"
+ << "fmax is the third decay argument in the .dec file. VectorIsr attempts to set it reasonably if it wasn't provided\n"
+ << "To determine a more appropriate value, build GeneratorQAApp, and set the third argument for this decay <0.\n"
+ << "If you haven't been providing the first 2 arguments, set them to be 1. 1.). The program will report\n"
+ << "the largest weight it finds. You should set fmax to be slightly larger.\n"
+ << "Alternatively try the following values for various vector particles: "
+ << "phi->1.15 J/psi-psi(4415)->0.105\n"
+ << "The current value of f and fmax for " << EvtPDL::name(phi->getId()) << " are " << f << " " << fmax << "\n"
+ << "Will now assert\n";
+ assert(0);
+ }
+
+
+ if (fmax > 0.) {
+ fran = fmax*EvtRandom::Flat(0.0,1.0);
+ }
+
+ else {
+ //determine max weight for this vector particle mass
+ if (f>largest_f) {
+ largest_f = f;
+ report(INFO,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " "
+ << "vector_mass "
+ << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f
+ << ". f_col cs_B = " << f_col << " " << cs_Born
+ << std::endl;
+ }
+ if (m%10000 == 0) {
+ report(INFO,"EvtGen") << m << " " << EvtPDL::name(phi->getId()) << " "
+ << "vector_mass "
+ << " " << EvtPDL::getMeanMass(phi->getId()) << " fmax should be at least " << largest_f
+ << ". f_col cs_B = " << f_col << " " << cs_Born
+ << std::endl;
+ }
+
+ f_col = 0.;
+ f = 0.;
+ //determine max weight for this vector particle mass
+ }
+
+ if (m > 100000){
+
+ if (fmax > 0.) report(INFO,"EvtGen") << "EvtVectorIsr is having problems. Check the fmax value - the 3rd argument in the .dec file\n"
+ << "Recommended values for various vector particles: "
+ << "phi->1.15 J/psi-psi(4415)->0.105 "
+ << "Upsilon(1S,2S,3S)->0.14\n";
+ assert(0);
+ }
+ }//while (fran > f)
+
+ }//if (firstorder)
+
+ //Compute parameters for boost to/from the system after colinear radiation
+
+ double bet_l;
+ double gam_l;
+ double betgam_l;
+
+ double csfrmn_new;
+ double csbkmn_new;
+
+ if (firstorder){
+ bet_l = 0.;
+ gam_l = 1.;
+ betgam_l = 0.;
+ csfrmn_new = csfrmn;
+ csbkmn_new = csbkmn;
+ } else {
+ double xx = e02/e01;
+ double sq_xx = sqrt(xx);
+ bet_l = (1.-xx)/(1.+xx);
+ gam_l = (1.+xx)/(2.*sq_xx);
+ betgam_l = (1.-xx)/(2.*sq_xx);
+
+ //Boost photon cos_theta limits in lab to limits in the system after colinear rad
+ csfrmn_new=(csfrmn - bet_l)/(1. - bet_l*csfrmn);
+ csbkmn_new=(csbkmn - bet_l)/(1. - bet_l*csbkmn);
+ }
+
+// //generate kinematics according to Bonneau-Martin article
+// //Nucl. Phys. B27 (1971) 381-397
+
+ // For backward compatibility with .dec files before SP5, the backward cos limit for
+ //the ISR photon is actually given as *minus* the actual limit. Sorry, this wouldn't be
+ //my choice. -Joe
+
+ //gamma momentum in the vpho restframe *after* soft colinear radiation
+ double pg = (s_new - phi->mass()*phi->mass())/(2.*wcm_new);
+
+
+ //calculate the beta of incoming electrons after colinear rad in the frame where e= and e- have equal momentum
+ double beta=electMass/ebeam; //electMass/Ebeam = 1/gamma
+ beta=sqrt(1. - beta*beta); //sqrt (1 - (1/gamma)**2)
+
+ double ymax=log((1.+beta*csfrmn_new)/(1.-beta*csfrmn_new));
+ double ymin=log((1.-beta*csbkmn_new)/(1.+beta*csbkmn_new));
+
+ // photon theta distributed as 2*beta/(1-beta**2*cos(theta)**2)
+ double y=(ymax-ymin)*EvtRandom::Flat(0.0,1.0) + ymin;
+ double cs=exp(y);
+ cs=(cs - 1.)/(cs + 1.)/beta;
+ double sn=sqrt(1-cs*cs);
+
+ double fi=EvtRandom::Flat(EvtConst::twoPi);
+
+ //four-vector for the phi
+ double phi_p0 = sqrt(phi->mass()*phi->mass()+pg*pg);
+ double phi_p3 = -pg*cs;
+
+
+ //boost back to frame before colinear radiation.
+ EvtVector4R p4phi(gam_l*phi_p0 + betgam_l*phi_p3,
+ -pg*sn*cos(fi),
+ -pg*sn*sin(fi),
+ betgam_l*phi_p0 + gam_l*phi_p3);
+
+ double isr_p0 = pg;
+ double isr_p3 = -phi_p3;
+ EvtVector4R p4gamma(gam_l*isr_p0 + betgam_l*isr_p3,
+ -p4phi.get(1),
+ -p4phi.get(2),
+ betgam_l*isr_p0 + gam_l*isr_p3);
+
+
+ //four-vectors of the collinear photons
+ if (!firstorder) {
+ p4softg1.set(0, eb-e02); p4softg1.set(3, e02-eb);
+ p4softg2.set(0, eb-e01); p4softg2.set(3, eb-e01);
+ }
+
+ //save momenta for particles
+ phi->init( getDaug(0),p4phi);
+ gamma->init( getDaug(1),p4gamma);
+
+
+ //add the two colinear photons as vphoton daughters
+ EvtPhotonParticle *softg1=new EvtPhotonParticle;;
+ EvtPhotonParticle *softg2=new EvtPhotonParticle;;
+ softg1->init(gammaId,p4softg1);
+ softg2->init(gammaId,p4softg2);
+ softg1->addDaug(p);
+ softg2->addDaug(p);
+
+ //try setting the spin density matrix of the phi
+ //get polarization vector for phi in its parents restframe.
+ EvtVector4C phi0=phi->epsParent(0);
+ EvtVector4C phi1=phi->epsParent(1);
+ EvtVector4C phi2=phi->epsParent(2);
+
+ //get polarization vector for a photon in its parents restframe.
+ EvtVector4C gamma0=gamma->epsParentPhoton(0);
+ EvtVector4C gamma1=gamma->epsParentPhoton(1);
+
+ EvtComplex r1p=phi0*gamma0;
+ EvtComplex r2p=phi1*gamma0;
+ EvtComplex r3p=phi2*gamma0;
+
+
+ EvtComplex r1m=phi0*gamma1;
+ EvtComplex r2m=phi1*gamma1;
+ EvtComplex r3m=phi2*gamma1;
+
+ EvtComplex rho33=r3p*conj(r3p)+r3m*conj(r3m);
+ EvtComplex rho22=r2p*conj(r2p)+r2m*conj(r2m);
+ EvtComplex rho11=r1p*conj(r1p)+r1m*conj(r1m);
+
+ EvtComplex rho13=r3p*conj(r1p)+r3m*conj(r1m);
+ EvtComplex rho12=r2p*conj(r1p)+r2m*conj(r1m);
+ EvtComplex rho23=r3p*conj(r2p)+r3m*conj(r2m);
+
+ EvtComplex rho31=conj(rho13);
+ EvtComplex rho32=conj(rho23);
+ EvtComplex rho21=conj(rho12);
+
+
+ EvtSpinDensity rho;
+ rho.setDim(3);
+
+ rho.set(0,0,rho11);
+ rho.set(0,1,rho12);
+ rho.set(0,2,rho13);
+ rho.set(1,0,rho21);
+ rho.set(1,1,rho22);
+ rho.set(1,2,rho23);
+ rho.set(2,0,rho31);
+ rho.set(2,1,rho32);
+ rho.set(2,2,rho33);
+
+ setDaughterSpinDensity(0);
+ phi->setSpinDensityForward(rho);
+
+ return ;
+}
+
+double EvtVectorIsr::ckhrad1(double xx, double a, double b){
+ //port of AfkQed/ckhrad.F function ckhrad1
+ double yy = xx*xx;
+ double zz = 1. - 2*xx + yy;
+ return 0.5* (1. + yy + zz/(a-1.) + 0.25*b*( -0.5*(1. + 3*yy)*log(xx)) - zz );
+}
+
+void EvtVectorIsr::ckhrad(const double& e_beam,const double& q2_min,double& e01,double& e02,double& f){
+ //port of AfkQed/ckhrad.F subroutine ckhrad
+ const double adp = 1. / 137.0359895 / EvtConst::pi;
+ const double pi2 = EvtConst::pi*EvtConst::pi;
+ // const double dme = 0.00051099906;
+ const double dme = EvtPDL::getMeanMass(EvtPDL::getId("e-"));
+
+ double r1=EvtRandom::Flat();//Generates Flat from 0 - 1
+ double r2=EvtRandom::Flat();
+
+ double sss = 4.*e_beam*e_beam;
+ double biglog = log(sss/(dme*dme));
+ double beta = 2.*adp*(biglog - 1.);
+ double betae_lab = beta;
+ double p3 = adp*(pi2/3. - 0.5);
+ double p12 = adp*adp * (11./8. - 2.*pi2/3.);
+ double coefener = 1. + 0.75*betae_lab + p3;
+ double coef1 = coefener + 0.125*pi2*beta*beta;
+ double coef2 = p12* biglog*biglog;
+ double facts = coef1 + coef2;
+
+ double y1_min = 0;
+ double e1min = 0.25 * q2_min/e_beam;
+ double y1_max = pow( 1. - e1min/e_beam, 0.5*beta );
+ double y1 = y1_min +r1 *(y1_max - y1_min);
+ e01 = e_beam *(1. - pow(y1, 2./beta) );
+
+ double y2_min = 0.;
+ double e2min = 0.25 * q2_min/e01;
+ double y2_max = pow( 1. - e2min/e_beam, 0.5*beta);
+ double y2 = y2_min +r2 *(y2_max - y2_min);
+ e02 = e_beam *(1. - pow(y2, 2./beta) );
+
+
+ double xx1 = e01/e_beam;
+ double xx2 = e02/e_beam;
+
+ f = y1_max * y2_max * ckhrad1(xx1,biglog,betae_lab) * ckhrad1(xx2,biglog,betae_lab) * facts;
+
+ return;
+ }
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff