TEvtGen/EvtGenModels/EvtVectorIsr.cpp

   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