STARLIGHT/starlight/src/.svn/text-base/gammagammasingle.cpp.svn-base

   1 ///////////////////////////////////////////////////////////////////////////
   2 //
   3 //    Copyright 2010
   4 //
   5 //    This file is part of starlight.
   6 //
   7 //    starlight is free software: you can redistribute it and/or modify
   8 //    it under the terms of the GNU General Public License as published by
   9 //    the Free Software Foundation, either version 3 of the License, or
  10 //    (at your option) any later version.
  11 //
  12 //    starlight is distributed in the hope that it will be useful,
  13 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
  14 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
  15 //    GNU General Public License for more details.
  16 //
  17 //    You should have received a copy of the GNU General Public License
  18 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
  19 //
  20 ///////////////////////////////////////////////////////////////////////////
  21 //
  22 // File and Version Information:
  23 // $Rev::                             $: revision of last commit
  24 // $Author::                          $: author of last commit
  25 // $Date::                            $: date of last commit
  26 //
  27 // Description:
  28 //
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33 #include <iostream>
  34 #include <fstream>
  35 #include <cmath>
  36 #include <vector>
  37
  38
  39 #include "starlightconstants.h"
  40 #include "gammagammasingle.h"
  41 #include "starlightconfig.h"
  42
  43 using namespace std;
  44
  45
  46 //______________________________________________________________________________
  47 Gammagammasingle::Gammagammasingle(beamBeamSystem& bbsystem)
  48 : eventChannel(bbsystem)
  49 #ifdef ENABLE_PYTHIA
  50 ,_pyDecayer()
  51 #endif
  52 {
  53
  54 #ifdef ENABLE_PYTHIA
  55     _pyDecayer.init();
  56 #endif
  57
  58   //Storing inputparameters into protected members for use
  59   _GGsingInputnumw=inputParametersInstance.nmbWBins();
  60   _GGsingInputnumy=inputParametersInstance.nmbRapidityBins();
  61   _GGsingInputpidtest=inputParametersInstance.prodParticleType();
  62   _GGsingInputGamma_em=inputParametersInstance.beamLorentzGamma();
  63   cout<<"SINGLE MESON pid test: "<<_GGsingInputpidtest<<endl;
  64   //reading in luminosity tables
  65   read();
  66   //Now calculating crosssection
  67   singleCrossSection();
  68 }
  69
  70
  71 //______________________________________________________________________________
  72 Gammagammasingle::~Gammagammasingle()
  73 { }
  74
  75
  76 //______________________________________________________________________________
  77 void Gammagammasingle::singleCrossSection()
  78 {
  79   //This function carries out a delta function cross-section calculation. For reference, see STAR Note 243, Eq. 8
  80   //Multiply all _Farray[] by _f_max
  81   double _sigmaSum=0.,remainw=0.;//_remainwd=0.;
  82   int ivalw =0;//_ivalwd;
  83   //calculate the differential cross section and place in the sigma table
  84   _wdelta=getMass();
  85   for(int i=0;i<_GGsingInputnumw;i++){
  86     for(int j=0;j<_GGsingInputnumy;j++){
  87       // Eq. 1 of starnote 347
  88       _sigmax[i][j]=(getSpin()*2.+1.)*4*starlightConstants::pi*starlightConstants::pi*getWidth()/
  89         (getMass()*getMass()*getMass())*_f_max*_Farray[i][j]*starlightConstants::hbarc*starlightConstants::hbarc/100.;
  90     }
  91   }
  92   //find the index, i,for the value of w just less than the mass because we want to use the value from the sigma table that has w=mass
  93
  94   for(int i=0;i<_GGsingInputnumw;i++){
  95     if(getMass()>_Warray[i]) ivalw=i;
  96   }
  97
  98   remainw = (getMass()-_Warray[ivalw])/(_Warray[ivalw+1]-_Warray[ivalw+1]-_Warray[ivalw]);
  99   _ivalwd = ivalw;
 100   _remainwd = remainw;
 101   //if we are interested rho pairs at threshold, the just set sigma to 100nb
 102   switch(_GGsingInputpidtest){
 103   case starlightConstants::ZOVERZ03:
 104     _sigmaSum =0.;
 105     for(int j=0;j<_GGsingInputnumy-1;j++){
 106                         _sigmaSum = _sigmaSum +(_Yarray[j+1]-_Yarray[j])*
 107                           100.0E-9*(.1/getMass())*((1.-remainw)*_f_max*
 108                                                    (_Farray[ivalw][j]+_Farray[ivalw][j])/2.+remainw*_f_max*
 109                                                    (_Farray[ivalw+1][j]+_Farray[ivalw+1][j+1])/2.);
 110     }
 111     break;
 112   default:
 113     //Sum to find the total cross-section
 114     _sigmaSum =0.;
 115     for(int j =0;j<_GGsingInputnumy-1;j++){
 116                         _sigmaSum = _sigmaSum+
 117                           (_Yarray[j+1]-_Yarray[j])*((1.-remainw)*
 118                                                    (_sigmax[ivalw][j]+_sigmax[ivalw][j+1])/2.+remainw*
 119                                                    (_sigmax[ivalw+1][j]+_sigmax[ivalw+1][j+1])/2.);
 120     }
 121   }
 122   if(_sigmaSum > 0.1) cout <<"The total cross-section is: "<<_sigmaSum<<" barns."<<endl;
 123   else if(_sigmaSum > 0.0001)cout <<"The total cross-section is: "<<_sigmaSum*1000<<" mb."<<endl;
 124   else cout <<"The total cross-section is: "<<_sigmaSum*1000000<<" ub."<<endl;
 125
 126
 127   return;
 128 }
 129
 130
 131 //______________________________________________________________________________
 132 void Gammagammasingle::pickw(double &w)
 133 {
 134   //This function picks a w for the 2-photon calculation.
 135   double sgf=0.,signorm=0.,x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.;
 136   int ivalw=0;
 137
 138   double * _sigofw;
 139   double * sgfint;
 140   _sigofw = new double[starlightLimits::MAXWBINS];
 141   sgfint = new double[starlightLimits::MAXYBINS];
 142
 143   if(_wdelta != 0){
 144     w=_wdelta;
 145     ivalw=_ivalwd;
 146     remainw=_remainwd;
 147   }
 148   else{
 149     //deal with the case where sigma is an array
 150     //_sigofw is simga integrated over y using a linear interpolation
 151     //sigint is the integral of sgfint, normalized
 152
 153     //integrate sigma down to a function of just w
 154     for(int i=0;i<_GGsingInputnumw;i++){
 155       _sigofw[i]=0.;
 156       for(int j=0;j<_GGsingInputnumy-1;j++){
 157         _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.;
 158       }
 159     }
 160     //calculate the unnormalized sgfint array
 161     sgfint[0]=0.;
 162     for(int i=0;i<_GGsingInputnumw-1;i++){
 163       sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.;
 164       sgfint[i+1]=sgfint[i]+sgf;
 165     }
 166     //normalize sgfint array
 167     signorm=sgfint[_GGsingInputnumw-1];
 168
 169     for(int i=0;i<_GGsingInputnumw;i++){
 170       sgfint[i]=sgfint[i]/signorm;
 171     }
 172     //pick a random number
 173     x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 174     //compare x and sgfint to find the ivalue which is just less than the random number x
 175     for(int i=0;i<_GGsingInputnumw;i++){
 176       if(x > sgfint[i]) ivalw=i;
 177     }
 178     //remainder above ivalw
 179     remainarea = x - sgfint[ivalw];
 180
 181     //figure out what point corresponds to the excess area in remainarea
 182     c = -remainarea*signorm/(_Warray[ivalw+1]-_Warray[ivalw]);
 183     b = _sigofw[ivalw];
 184     a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.;
 185     if(a==0.){
 186       remainw = -c/b;
 187     }
 188     else{
 189       remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a);
 190     }
 191     _ivalwd = ivalw;
 192     _remainwd = remainw;
 193     //calculate the w value
 194     w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw;
 195     }
 196
 197   delete[] _sigofw;
 198   delete[] sgfint;
 199 }
 200
 201
 202 //______________________________________________________________________________
 203 void Gammagammasingle::picky(double &y)
 204 {
 205   double * sigofy;
 206   double * sgfint;
 207   sigofy = new double[starlightLimits::MAXYBINS];
 208   sgfint = new double[starlightLimits::MAXYBINS];
 209
 210   double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.;
 211   int ivalw=0,ivaly=0;
 212
 213   ivalw=_ivalwd;
 214   remainw=_remainwd;
 215   //average over two colums to get y array
 216   for(int j=0;j<_GGsingInputnumy;j++){
 217     sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw;
 218   }
 219   //calculate the unnormalized sgfint
 220
 221   sgfint[0]=0.;
 222   for(int j=0;j<_GGsingInputnumy-1;j++){
 223     sgf = (sigofy[j+1]+sigofy[j])/2.;
 224     sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]);
 225   }
 226
 227   //normalize the sgfint array
 228   signorm = sgfint[_GGsingInputnumy-1];
 229
 230   for(int j=0;j<_GGsingInputnumy;j++){
 231     sgfint[j]=sgfint[j]/signorm;
 232   }
 233   //pick a random number
 234   x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 235   //compare x and sgfint to find the ivalue which is just less then the random number x
 236   for(int i=0;i<_GGsingInputnumy;i++){
 237     if(x > sgfint[i])
 238       ivaly = i;
 239   }
 240   //remainder above ivaly
 241   remainarea = x - sgfint[ivaly];
 242   //figure what point corresponds to the leftover area in remainarea
 243   c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]);
 244   b = sigofy[ivaly];
 245   a = (sigofy[ivaly+1]-sigofy[ivaly])/2.;
 246   if(a==0.){
 247     remainy = -c/b;
 248   }
 249   else{
 250     remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a);
 251   }
 252   //calculate the y value
 253   y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy;
 254   delete[] sigofy;
 255   delete[] sgfint;
 256 }
 257
 258
 259 //______________________________________________________________________________
 260 void Gammagammasingle::parentMomentum(double w,double y,double &E,double &px,double &py,double &pz)
 261 {
 262   //this function calculates px,py,pz,and E given w and y
 263   double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.;
 264
 265   //E1 and E2 are for the 2 photons in the CM frame
 266   E1 = w*exp(y)/2.;
 267   E2 = w*exp(-y)/2.;
 268   //pz = E1-E2;
 269   //calculate px and py
 270   //to get x and y components-- phi is random between 0 and 2*pi
 271   anglepp1 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 272   anglepp2 = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 273
 274   pp1 = pp(E1);
 275   pp2 = pp(E2);
 276   px = pp1*cos(2.*starlightConstants::pi*anglepp1)+pp2*cos(2.*starlightConstants::pi*anglepp2);
 277   py = pp1*sin(2.*starlightConstants::pi*anglepp1)+pp2*sin(2.*starlightConstants::pi*anglepp2);
 278   //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle
 279   pt = sqrt(px*px+py*py);
 280   //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz
 281   E = sqrt(w*w+pt*pt)*cosh(y);
 282   pz= sqrt(w*w+pt*pt)*sinh(y);
 283   signpx = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 284   //pick the z direction
 285   if(signpx > 0.5)
 286     pz = -pz;
 287 }
 288
 289
 290 //______________________________________________________________________________
 291 double Gammagammasingle::pp(double E)
 292 {
 293   //  will probably have to pass in beambeamsys? that way we can do beam1.formFactor(t) or beam2..., careful with the way sergey did it for asymmetry
 294   //  returns on random draw from pp(E) distribution
 295
 296   double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
 297   double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
 298   int satisfy =0;
 299
 300   ereds = (E/_GGsingInputGamma_em)*(E/_GGsingInputGamma_em);
 301   Cm = sqrt(3.)*E/_GGsingInputGamma_em;
 302   //the amplitude of the p_t spectrum at the maximum
 303   singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
 304   //Doing this once and then storing it as a double, which we square later...SYMMETRY?using beam1 for now.
 305   Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
 306
 307   //pick a test value pp, and find the amplitude there
 308   x = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 309   pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); //Will use nucleus #1, there should be two for symmetry//nextline
 310   singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
 311   test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
 312
 313   while(satisfy==0){
 314     u = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 315     if(u*Coef <= test){
 316       satisfy =1;
 317     }
 318     else{
 319       x =randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 320       pp = 5*starlightConstants::hbarc/_bbs.beam1().nuclearRadius()*x;
 321       singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);//Symmetry
 322       test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
 323     }
 324   }
 325   return pp;
 326 }
 327
 328
 329 //______________________________________________________________________________
 330 void Gammagammasingle::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double /*E*/,double W,double px0,double py0,double pz0,double &px1,double &py1,double &pz1,double &px2,double &py2,double &pz2,int &iFbadevent)
 331 {
 332   //     This routine decays a particle into two particles of mass mdec,
 333   //     taking spin into account
 334
 335   double mdec=0.,E1=0.,E2=0.;
 336   double pmag,ytest=0.;
 337   double phi,theta,xtest,dndtheta,Ecm;
 338   double  betax,betay,betaz;
 339
 340   //    set the mass of the daughter particles
 341   switch(_GGsingInputpidtest){
 342   case starlightConstants::ZOVERZ03:
 343   case starlightConstants::F2:
 344     mdec = starlightConstants::pionChargedMass;
 345     break;
 346   case starlightConstants::F2PRIME:
 347     //  decays 50% to K+/K-, 50% to K_0's
 348     ytest = randyInstance.Rndom();
 349     if(ytest >= 0.5){
 350       mdec = starlightConstants::kaonChargedMass;
 351     }
 352     else{
 353       mdec = 0.493677;
 354     }
 355     break;
 356   default :
 357     cout<<"No default mass selected for single photon-photon particle, expect errant results"<<endl;
 358   }
 359
 360   //Calculating the momentum's magnitude
 361   //add switch for rho pairs at threshold and everything else.
 362   switch(_GGsingInputpidtest){
 363   case starlightConstants::ZOVERZ03:    //the rho pairs produced at threshold
 364     pmag = sqrt(getMass()*getMass()/4. - mdec*mdec);
 365     break;
 366   default :
 367     if(W < 2*mdec){
 368       cout<<" ERROR: W="<<W<<endl;
 369       iFbadevent = 1;
 370       return;
 371     }
 372     pmag = sqrt(W*W/4. - mdec*mdec);
 373   }
 374   //     pick an orientation, based on the spin
 375   //      phi has a flat distribution in 2*pi
 376   phi = randyInstance.Rndom()*2.*starlightConstants::pi; //(random()/(RAND_MAX+1.0))* 2.*starlightConstants::pi;
 377
 378   //     find theta, the angle between one of the outgoing particles and
 379   //    the beamline, in the frame of the two photons
 380   //this will depend on spin, F2,F2' and z/z03 all have spin 2, all other photonphoton-single mesons are handled by jetset
 381   //Applies to spin2 mesons.
 382  L300td:
 383   theta = starlightConstants::pi*randyInstance.Rndom();
 384   xtest = randyInstance.Rndom();
 385   dndtheta = sin(theta)*sin(theta)*sin(theta)*sin(theta)*sin(theta);
 386   if(xtest > dndtheta)
 387     goto L300td;
 388
 389   //     compute unboosted momenta
 390   px1 = sin(theta)*cos(phi)*pmag;
 391   py1 = sin(theta)*sin(phi)*pmag;
 392   pz1 = cos(theta)*pmag;
 393   px2 = -px1;
 394   py2 = -py1;
 395   pz2 = -pz1;
 396   //        compute energies
 397   //Changed mass to W 11/9/2000 SRK
 398   Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
 399   E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
 400   E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
 401
 402   //     Lorentz transform into the lab frame
 403   // betax,betay,betaz are the boost of the complete system
 404   betax = -(px0/Ecm);
 405   betay = -(py0/Ecm);
 406   betaz = -(pz0/Ecm);
 407
 408   transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
 409   transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
 410
 411
 412   if(iFbadevent == 1)
 413     return;
 414
 415   //       change particle id from that of parent to that of daughters
 416
 417   switch(_GGsingInputpidtest){
 418     //These decay into a pi+ pi- pair
 419   case starlightConstants::ZOVERZ03:
 420   case starlightConstants::F2:
 421     ipid=starlightConstants::PION;
 422     break;
 423   case starlightConstants::F2PRIME:
 424     if( ytest >= 0.5 )
 425       {
 426         //Decays 50/50 into k+ k- or k_s k_l
 427         ipid=starlightConstants::KAONCHARGE;
 428       }
 429     else
 430       {
 431         ipid=starlightConstants::KAONNEUTRAL;
 432       }
 433     break;
 434   default:
 435     cout<<"Rethink the daughter particles"<<endl;
 436   }
 437 }
 438
 439
 440 //______________________________________________________________________________
 441 starlightConstants::event Gammagammasingle::produceEvent(int &/*ievent*/)
 442 {
 443   // Not in use anymore, default event struct returned
 444   return starlightConstants::event();
 445 }
 446
 447
 448 //______________________________________________________________________________
 449 // fix it ... lost functionality
 450 //starlightConstants::event Gammagammasingle::produceEvent(int &ievent)
 451 upcEvent Gammagammasingle::produceEvent()
 452 {
 453   //     cout << "NOT IMPLEMENTED!" << endl;
 454
 455   //     return upcEvent();
 456
 457   //    returns the vector with the decay particles inside.
 458   //    onedecayparticle single;
 459   starlightConstants::event single;
 460   double comenergy = 0.;
 461   double rapidity = 0.;
 462   double parentE = 0.;
 463   double parentmomx=0.,parentmomy=0.,parentmomz=0.;
 464
 465   //this function decays particles and writes events to a file
 466   //zeroing out the event structure
 467   single._numberOfTracks=0;
 468   for(int i=0;i<4;i++){
 469     single.px[i]=0.;
 470     single.py[i]=0.;
 471     single.pz[i]=0.;
 472     single._fsParticle[i]=starlightConstants::UNKNOWN;
 473     single._charge[i]=0;
 474   }
 475
 476   pickw(comenergy);
 477   picky(rapidity);
 478   parentMomentum(comenergy,rapidity,parentE,parentmomx,parentmomy,parentmomz);
 479
 480
 481   if(_GGsingInputpidtest != starlightConstants::F2 && _GGsingInputpidtest != starlightConstants::F2PRIME)
 482   {
 483 #ifdef ENABLE_PYTHIA
 484     starlightParticle particle(parentmomx,parentmomy,parentmomz, parentE, getMass(),_GGsingInputpidtest , 0);
 485
 486     _pyDecayer.addParticle(particle);
 487
 488     return _pyDecayer.execute();
 489 #endif
 490   }
 491
 492
 493   int ievent = 0;
 494   int iFbadevent=0;
 495   starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
 496   double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
 497   double px3=0.,px4=0.,py3=0.,py4=0.,pz3=0.,pz4=0.;
 498   //  double theta=0.,phi=0.;//angles from jetset
 499   double xtest=0.,ztest=0.;
 500   switch(_GGsingInputpidtest){
 501   case starlightConstants::ZOVERZ03:
 502     //Decays into two pairs.
 503     parentmomx=parentmomx/2.;
 504     parentmomy=parentmomy/2.;
 505     parentmomz=parentmomz/2.;
 506     //Pair #1
 507     twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 508     //Pair #2
 509     twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px3,py3,pz3,px4,py4,pz4,iFbadevent);
 510     //Now add them to vectors to be written out later.
 511
 512     single._numberOfTracks=4;//number of tracks per event
 513     if (iFbadevent==0){
 514       xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 515       ztest = randyInstance.Rndom();
 516       //Assigning charges randomly.
 517       if (xtest<0.5){
 518         single._charge[0]=1;//q1=1;
 519         single._charge[1]=-1;//q2=-1;
 520       }
 521       else{
 522         single._charge[0]=1;//q1=-1;
 523         single._charge[1]=-1;//q2=1;
 524       }
 525       if (ztest<0.5){
 526         single._charge[2]=1;//q3=1;
 527         single._charge[3]=-1;//q4=-1;
 528       }
 529       else{
 530         single._charge[2]=-1;//q3=-1;
 531         single._charge[3]=1;//q4=1;
 532       }
 533       //Track #1
 534       single.px[0]=px1;
 535       single.py[0]=py1;
 536       single.pz[0]=pz1;
 537       single._fsParticle[0]=ipid;
 538       //Track #2
 539       single.px[1]=px2;
 540       single.py[1]=py2;
 541       single.pz[1]=pz2;
 542       single._fsParticle[1]=ipid;
 543       //Track #3
 544       single.px[2]=px3;
 545       single.py[2]=py3;
 546       single.pz[2]=pz3;
 547       single._fsParticle[2]=ipid;
 548       //Track #4
 549       single.px[3]=px4;
 550       single.py[3]=py4;
 551       single.pz[3]=pz4;
 552       single._fsParticle[3]=ipid;
 553
 554       ievent=ievent+1;
 555     }
 556
 557     break;
 558   case starlightConstants::F2:
 559   case starlightConstants::F2PRIME:
 560     twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 561
 562     single._numberOfTracks=2;
 563     if (iFbadevent==0){
 564       xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 565       if (xtest<0.5){
 566         single._charge[0]=1;//q1=1;
 567         single._charge[1]=-1;//q2=-1;
 568       }
 569       else{
 570         single._charge[0]=-1;//q1=-1;
 571         single._charge[1]=1;//q2=1;
 572       }
 573       //Track #1
 574       single.px[0]=px1;
 575       single.py[0]=py1;
 576       single.pz[0]=pz1;
 577       single._fsParticle[0]=ipid*single._charge[0];
 578       //Track #2
 579       single.px[1]=px2;
 580       single.py[1]=py2;
 581       single.pz[1]=pz2;
 582       single._fsParticle[1]=ipid*single._charge[1];
 583       ievent=ievent+1;
 584     }
 585     break;
 586   default:
 587     break;
 588   }
 589
 590   return upcEvent(single);
 591 }
 592
 593
 594 //______________________________________________________________________________
 595 void Gammagammasingle::thephi(double W,double px,double py,double pz,double E,double &theta,double &phi)
 596 {
 597   //     This subroutine calculates angles for channels decayed by jetset.
 598   //    subroutine thephi(W,px,py,pz,E,theta,phi)
 599   E = sqrt (W*W+px*px+py*py+pz*pz);
 600
 601   theta = acos(pz/sqrt(px*px+py*py+pz*pz));
 602   phi = acos(px/sqrt(px*px+py*py));
 603
 604   if ((px == 0)  && (py == 0))
 605     phi = 0.;
 606   if (py < 0)
 607     phi = 2*starlightConstants::pi - phi;
 608 }
 609
 610
 611 //______________________________________________________________________________
 612 double Gammagammasingle::getMass()
 613 {
 614   using namespace starlightConstants;
 615   double singlemass=0.;
 616   switch(_GGsingInputpidtest){
 617   case starlightConstants::ETA:
 618     singlemass= etaMass;
 619     break;
 620   case starlightConstants::ETAPRIME:
 621     singlemass=etaPrimeMass;
 622     break;
 623   case starlightConstants::ETAC:
 624     singlemass=etaCMass;
 625     break;
 626   case starlightConstants::F0:
 627     singlemass=f0Mass;
 628     break;
 629   case starlightConstants::F2:
 630     singlemass=f2Mass;
 631     break;
 632   case starlightConstants::A2:
 633     singlemass=a2Mass;
 634     break;
 635   case starlightConstants::F2PRIME:
 636     singlemass=f2PrimeMass;
 637     break;
 638   case starlightConstants::ZOVERZ03:
 639     singlemass=1.540;
 640     break;
 641   default:
 642     cout<<"Not a recognized single particle, Gammagammasingle::getmass(), mass = 0."<<endl;
 643   }
 644   return singlemass;
 645 }
 646
 647
 648 //______________________________________________________________________________
 649 double Gammagammasingle::getWidth()
 650 {
 651   double singlewidth=0.;
 652   switch(_GGsingInputpidtest){
 653   case starlightConstants::ETA:
 654     singlewidth=1.E-6;
 655     break;
 656   case starlightConstants::ETAPRIME:
 657     singlewidth=5.E-6;
 658     break;
 659   case starlightConstants::ETAC:
 660     singlewidth=6.4E-6;
 661     break;
 662   case starlightConstants::F0:
 663     singlewidth=0.56E-6;
 664     break;
 665   case starlightConstants::F2:
 666     singlewidth=2.6E-6;
 667     break;
 668   case starlightConstants::A2:
 669     singlewidth=1.04E-6;
 670     break;
 671   case starlightConstants::F2PRIME:
 672     singlewidth=0.1E-6;
 673     break;
 674   case starlightConstants::ZOVERZ03:
 675     singlewidth=0.1E-6;
 676     break;
 677   default:
 678     cout<<"Not a recognized single particle, Gammagammasingle::getwidth(), width = 0."<<endl;
 679   }
 680   return singlewidth;
 681 }
 682
 683
 684 //______________________________________________________________________________
 685 double Gammagammasingle::getSpin()
 686 {
 687   double singlespin=0.5;
 688   switch(_GGsingInputpidtest){
 689   case starlightConstants::ETA:
 690     singlespin=0.0;
 691     break;
 692   case starlightConstants::ETAPRIME:
 693     singlespin=0.0;
 694     break;
 695   case starlightConstants::ETAC:
 696     singlespin=0.0;
 697     break;
 698   case starlightConstants::F0:
 699     singlespin=0.0;
 700     break;
 701   case starlightConstants::F2:
 702     singlespin=2.0;
 703     break;
 704   case starlightConstants::A2:
 705     singlespin=2.0;
 706     break;
 707   case starlightConstants::F2PRIME:
 708     singlespin=2.0;
 709     break;
 710   case starlightConstants::ZOVERZ03:
 711     singlespin=2.0;
 712     break;
 713   default:
 714     cout<<"Not a recognized single particle, Gammagammasingle::getspin(), spin = 0."<<endl;
 715   }
 716   return singlespin;
 717 }
 718
 719
 720