STARLIGHT/starlight/src/.svn/text-base/gammaavm.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 //    Added incoherent t2-> pt2 selection.  Following pp selection scheme
  29 //
  30 //
  31 ///////////////////////////////////////////////////////////////////////////
  32
  33
  34 #include <iostream>
  35 #include <fstream>
  36 #include <cassert>
  37 #include <cmath>
  38
  39 #include "gammaavm.h"
  40 #include "photonNucleusCrossSection.h"
  41 #include "wideResonanceCrossSection.h"
  42 #include "narrowResonanceCrossSection.h"
  43 #include "incoherentVMCrossSection.h"
  44
  45 using namespace std;
  46
  47
  48 //______________________________________________________________________________
  49 Gammaavectormeson::Gammaavectormeson(beamBeamSystem& bbsystem):eventChannel(bbsystem), _phaseSpaceGen(0)  //:readLuminosity(input),_bbs(bbsystem)
  50 {
  51         _VMNPT=inputParametersInstance.nmbPtBinsInterference();
  52         _VMWmax=inputParametersInstance.maxW();
  53         _VMWmin=inputParametersInstance.minW();
  54         _VMYmax=inputParametersInstance.maxRapidity();
  55         _VMYmin=-1.*_VMYmax;
  56         _VMnumw=inputParametersInstance.nmbWBins();
  57         _VMnumy=inputParametersInstance.nmbRapidityBins();
  58         _VMgamma_em=inputParametersInstance.beamLorentzGamma();
  59         _VMinterferencemode=inputParametersInstance.interferenceEnabled();
  60         _VMbslope=0.;//Will define in wide/narrow constructor
  61         _VMpidtest=inputParametersInstance.prodParticleType();
  62         _VMptmax=inputParametersInstance.maxPtInterference();
  63         _VMdpt=inputParametersInstance.ptBinWidthInterference();
  64         _VMCoherence=inputParametersInstance.coherentProduction();
  65         _VMCoherenceFactor=inputParametersInstance.coherentProduction();
  66         _ProductionMode=inputParametersInstance.productionMode();
  67
  68         switch(_VMpidtest){
  69         case starlightConstants::RHO:
  70         case starlightConstants::RHOZEUS:
  71                 _width=0.1507;
  72                 _mass=0.7685;
  73                 break;
  74         case starlightConstants::FOURPRONG:
  75                 // create n-body phase-space generator instance
  76                 _phaseSpaceGen = new nBodyPhaseSpaceGen();
  77                 _width = 0.360;
  78                 _mass  = 1.350;
  79                 break;
  80         case starlightConstants::OMEGA:
  81                 _width=0.00843;
  82                 _mass=0.78194;
  83                 break;
  84         case starlightConstants::PHI:
  85                 _width=0.00443;
  86                 _mass=1.019413;
  87                 break;
  88         case starlightConstants::JPSI:
  89         case starlightConstants::JPSI_ee:
  90         case starlightConstants::JPSI_mumu:
  91                 _width=0.000091;
  92                 _mass=3.09692;
  93                 break;
  94         case starlightConstants::JPSI2S:
  95         case starlightConstants::JPSI2S_ee:
  96         case starlightConstants::JPSI2S_mumu:
  97                 _width=0.000337;
  98                 _mass=3.686093;
  99                 break;
 100         case starlightConstants::UPSILON:
 101         case starlightConstants::UPSILON_ee:
 102         case starlightConstants::UPSILON_mumu:
 103                 _width=0.00005402;
 104                 _mass=9.46030;
 105                 break;
 106         case starlightConstants::UPSILON2S:
 107         case starlightConstants::UPSILON2S_ee:
 108         case starlightConstants::UPSILON2S_mumu:
 109                 _width=0.00003198;
 110                 _mass=10.02326;
 111                 break;
 112         case starlightConstants::UPSILON3S:
 113         case starlightConstants::UPSILON3S_ee:
 114         case starlightConstants::UPSILON3S_mumu:
 115                 _width=0.00002032;
 116                 _mass=10.3552;
 117                 break;
 118         default: cout<<"No proper vector meson defined, gammaavectormeson::gammaavectormeson"<<endl;
 119         }
 120 }
 121
 122
 123 //______________________________________________________________________________
 124 Gammaavectormeson::~Gammaavectormeson()
 125 {
 126         if (_phaseSpaceGen)
 127                 delete _phaseSpaceGen;
 128 }
 129
 130
 131 //______________________________________________________________________________
 132 void Gammaavectormeson::pickwy(double &W, double &Y)
 133 {
 134         double dW, dY, xw,xy,xtest;
 135         int  IW,IY;
 136
 137         dW = (_VMWmax-_VMWmin)/double(_VMnumw);
 138         dY = (_VMYmax-_VMYmin)/double(_VMnumy);
 139
 140  L201pwy:
 141
 142         xw = randyInstance.Rndom();// random()/(RAND_MAX+1.0);
 143         W = _VMWmin + xw*(_VMWmax-_VMWmin);
 144
 145         if (W < 2 * starlightConstants::pionChargedMass)
 146                 goto L201pwy;
 147
 148         IW = int((W-_VMWmin)/dW); //+ 1;
 149         xy = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 150         Y = _VMYmin + xy*(_VMYmax-_VMYmin);
 151         IY = int((Y-_VMYmin)/dY); //+ 1;
 152         xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 153
 154         if( xtest > _Farray[IW][IY] )
 155                 goto L201pwy;
 156
 157 }
 158
 159
 160 //______________________________________________________________________________
 161 void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
 162                                      double,  // E (unused)
 163                                      double  W,
 164                                      double  px0, double  py0, double  pz0,
 165                                      double& px1, double& py1, double& pz1,
 166                                      double& px2, double& py2, double& pz2,
 167                                      int&    iFbadevent)
 168 {
 169         // This routine decays a particle into two particles of mass mdec,
 170         // taking spin into account
 171
 172         double pmag;
 173         // double anglelep[20001],xtest,ytest=0.,dndtheta;
 174         double phi,theta,Ecm;
 175         double betax,betay,betaz;
 176         double mdec=0.0;
 177         double E1=0.0,E2=0.0;
 178
 179         //    set the mass of the daughter particles
 180         mdec=getDaughterMass(ipid);
 181
 182         //     calculate the magnitude of the momenta
 183         if(W < 2*mdec){
 184                 cout<<" ERROR: W="<<W<<endl;
 185                 iFbadevent = 1;
 186                 return;
 187         }
 188         pmag = sqrt(W*W/4. - mdec*mdec);
 189
 190         //     pick an orientation, based on the spin
 191         //      phi has a flat distribution in 2*pi
 192         phi = randyInstance.Rndom()*2.*starlightConstants::pi;//(random()/(RAND_MAX+1.0))* 2.*pi;
 193
 194         //     find theta, the angle between one of the outgoing particles and
 195         //    the beamline, in the frame of the two photons
 196
 197         theta=getTheta(ipid);
 198
 199         //     compute unboosted momenta
 200         px1 = sin(theta)*cos(phi)*pmag;
 201         py1 = sin(theta)*sin(phi)*pmag;
 202         pz1 = cos(theta)*pmag;
 203         px2 = -px1;
 204         py2 = -py1;
 205         pz2 = -pz1;
 206
 207         Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
 208         E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
 209         E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
 210
 211         betax = -(px0/Ecm);
 212         betay = -(py0/Ecm);
 213         betaz = -(pz0/Ecm);
 214
 215         transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
 216         transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
 217
 218         if(iFbadevent == 1)
 219                 return;
 220
 221 }
 222
 223
 224 //______________________________________________________________________________
 225 // decays a particle into four particles with isotropic angular distribution
 226 bool Gammaavectormeson::fourBodyDecay
 227 (starlightConstants::particleTypeEnum& ipid,
 228  const double                  ,           // E (unused)
 229  const double                  W,          // mass of produced particle
 230  const double*                 p,          // momentum of produced particle; expected to have size 3
 231  lorentzVector*                decayVecs,  // array of Lorentz vectors of daughter particles; expected to have size 4
 232  int&                          iFbadevent)
 233 {
 234         const double parentMass = W;
 235
 236         // set the mass of the daughter particles
 237         const double daughterMass = getDaughterMass(ipid);
 238         if (parentMass < 4 * daughterMass){
 239                 cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
 240                 iFbadevent = 1;
 241                 return false;
 242         }
 243
 244         // construct parent four-vector
 245         const double        parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
 246                                                 + parentMass * parentMass);
 247         const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
 248
 249         // setup n-body phase-space generator
 250         assert(_phaseSpaceGen);
 251         static bool firstCall = true;
 252         if (firstCall) {
 253                 const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
 254                 _phaseSpaceGen->setDecay(4, m);
 255                 // estimate maximum phase-space weight
 256                 _phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
 257                 firstCall = false;
 258         }
 259
 260         // generate phase-space event
 261         if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
 262                 return false;
 263
 264         // set Lorentzvectors of decay daughters
 265         for (unsigned int i = 0; i < 4; ++i)
 266                 decayVecs[i] = _phaseSpaceGen->daughter(i);
 267         return true;
 268 }
 269
 270
 271 //______________________________________________________________________________
 272 double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
 273 {
 274         //This will return the daughter particles mass, and the final particles outputed id...
 275         double ytest=0.,mdec=0.;
 276
 277         switch(_VMpidtest){
 278         case starlightConstants::RHO:
 279         case starlightConstants::RHOZEUS:
 280         case starlightConstants::FOURPRONG:
 281         case starlightConstants::OMEGA:
 282                 mdec = starlightConstants::pionChargedMass;
 283                 ipid = starlightConstants::PION;
 284                 break;
 285         case starlightConstants::PHI:
 286                 mdec = starlightConstants::kaonChargedMass;
 287                 ipid = starlightConstants::KAONCHARGE;
 288                 break;
 289         case starlightConstants::JPSI:
 290                 mdec = starlightConstants::mel;
 291                 ipid = starlightConstants::ELECTRON;
 292                 break;
 293         case starlightConstants::JPSI_ee:
 294                 mdec = starlightConstants::mel;
 295                 ipid = starlightConstants::ELECTRON;
 296                 break;
 297         case starlightConstants::JPSI_mumu:
 298                 mdec = starlightConstants::muonMass;
 299                 ipid = starlightConstants::MUON;
 300                 break;
 301         case starlightConstants::JPSI2S_ee:
 302                 mdec = starlightConstants::mel;
 303                 ipid = starlightConstants::ELECTRON;
 304                 break;
 305         case starlightConstants::JPSI2S_mumu:
 306                 mdec = starlightConstants::muonMass;
 307                 ipid = starlightConstants::MUON;
 308                 break;
 309
 310         case starlightConstants::JPSI2S:
 311         case starlightConstants::UPSILON:
 312         case starlightConstants::UPSILON2S:
 313         case starlightConstants::UPSILON3S:
 314                 //  decays 50% to e+/e-, 50% to mu+/mu-
 315                 ytest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 316
 317                 mdec = starlightConstants::muonMass;
 318                 ipid = starlightConstants::MUON;
 319                 break;
 320         case starlightConstants::UPSILON_ee:
 321         case starlightConstants::UPSILON2S_ee:
 322         case starlightConstants::UPSILON3S_ee:
 323                 mdec = starlightConstants::mel;
 324                 ipid = starlightConstants::ELECTRON;
 325                 break;
 326         case starlightConstants::UPSILON_mumu:
 327         case starlightConstants::UPSILON2S_mumu:
 328         case starlightConstants::UPSILON3S_mumu:
 329                 mdec = starlightConstants::muonMass;
 330                 ipid = starlightConstants::MUON;
 331                 break;
 332         default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
 333         }
 334
 335         return mdec;
 336 }
 337
 338
 339 //______________________________________________________________________________
 340 double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid)
 341 {
 342         //This depends on the decay angular distribution
 343         //Valid for rho, phi, omega.
 344         double theta=0.;
 345         double xtest=0.;
 346         double dndtheta=0.;
 347  L200td:
 348
 349         theta = starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 350         xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 351         //  Follow distribution for helicity +/-1
 352         //  Eq. 19 of J. Breitweg et al., Eur. Phys. J. C2, 247 (1998)
 353         //  SRK 11/14/2000
 354
 355         switch(ipid){
 356
 357         case starlightConstants::MUON:
 358         case starlightConstants::ELECTRON:
 359                 //primarily for upsilon/j/psi.  VM->ee/mumu
 360                 dndtheta = sin(theta)*(1.+((cos(theta))*(cos(theta))));
 361                 break;
 362
 363         case starlightConstants::PION:
 364         case starlightConstants::KAONCHARGE:
 365                 //rhos etc
 366                 dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
 367                 break;
 368
 369         default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
 370         }//end of switch
 371
 372         if(xtest > dndtheta)
 373                 goto L200td;
 374
 375         return theta;
 376
 377 }
 378
 379
 380 //______________________________________________________________________________
 381 double Gammaavectormeson::getWidth()
 382 {
 383         return _width;
 384 }
 385
 386
 387 //______________________________________________________________________________
 388 double Gammaavectormeson::getMass()
 389 {
 390         return _mass;
 391 }
 392
 393
 394 //______________________________________________________________________________
 395 double Gammaavectormeson::getSpin()
 396 {
 397         return 1.0; //VM spins are the same
 398 }
 399
 400
 401 //______________________________________________________________________________
 402 void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)
 403 {
 404         //     This subroutine calculates momentum and energy of vector meson
 405         //     given W and Y,   without interference.  Subroutine vmpt.f handles
 406         //     production with interference
 407
 408         double dW,dY;
 409         double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
 410         double px1,py1,px2,py2;
 411         double pt,xt,xtest,ytest;
 412         //      double photon_spectrum;
 413         double t2;
 414
 415         dW = (_VMWmax-_VMWmin)/double(_VMnumw);
 416         dY  = (_VMYmax-_VMYmin)/double(_VMnumy);
 417
 418         //Find Egam,Epom in CM frame
 419         if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
 420           if( _ProductionMode == 2 ){
 421             Egam = 0.5*W*exp(Y);
 422             Epom = 0.5*W*exp(-Y);
 423           }else{
 424             Egam = 0.5*W*exp(-Y);
 425             Epom = 0.5*W*exp(Y);
 426           }
 427         } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
 428           if( _ProductionMode == 2 ){
 429             Egam = 0.5*W*exp(-Y);
 430             Epom = 0.5*W*exp(Y);
 431           }else{
 432             Egam = 0.5*W*exp(Y);
 433             Epom = 0.5*W*exp(-Y);
 434           }
 435         } else {
 436           Egam = 0.5*W*exp(Y);
 437           Epom = 0.5*W*exp(-Y);
 438         }
 439
 440         // cout<<" Y: "<<Y<<" W: "<<W<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl;
 441         pt1 = pTgamma(Egam);
 442         phi1 = 2.*starlightConstants::pi*randyInstance.Rndom();
 443
 444         //      cout<<" pt1: "<<pt1<<endl;
 445
 446         if( (_bbs.beam1().A()==1 && _bbs.beam2().A()==1) ||
 447             (_ProductionMode == 4) ) {
 448             if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
 449               // Use dipole form factor for light VM
 450               L555vm:
 451               xtest = 2.0*randyInstance.Rndom();
 452               double ttest = xtest*xtest;
 453               ytest = randyInstance.Rndom();
 454               double t0 = 1./2.23;
 455               double yprob = xtest*_bbs.beam1().dipoleFormFactor(ttest,t0)*_bbs.beam1().dipoleFormFactor(ttest,t0);
 456               if( ytest > yprob ) goto L555vm;
 457               t2 = ttest;
 458               pt2 = xtest;
 459             }else{
 460                 //Use dsig/dt= exp(-_VMbslope*t) for heavy VM
 461                 xtest = randyInstance.Rndom();
 462                 t2 = (-1./_VMbslope)*log(xtest);
 463                 pt2 = sqrt(1.*t2);
 464             }
 465         } else {
 466                 //       >> Check tmin
 467                 tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
 468
 469                 if(tmin > 0.5){
 470                         cout<<" WARNING: tmin= "<<tmin<<endl;
 471                         cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl;
 472                         cout<<" Will pick a new W,Y "<<endl;
 473                         tcheck = 1;
 474                         return;
 475                 }
 476  L203vm:
 477                 xt = randyInstance.Rndom();
 478                 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
 479                   if( _ProductionMode == 2 ){
 480                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 481                   }else{
 482                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 483                   }
 484                 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
 485                   if( _ProductionMode == 2 ){
 486                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 487                   }else{
 488                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 489                   }
 490                 } else {
 491                     pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 492                 }
 493
 494                 xtest = randyInstance.Rndom();
 495                 t2 = tmin + pt2*pt2;
 496
 497                 if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
 498                         if(1.0 < _bbs.beam2().formFactor(t2)*pt2)  cout <<"POMERON"<<endl;
 499                         if( xtest > _bbs.beam2().formFactor(t2)*pt2) goto L203vm;
 500                 }
 501                 else{
 502
 503                   double comp=0.0;
 504                   if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
 505                     if( _ProductionMode == 2 ){
 506                       comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
 507                     }else{
 508                       comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
 509                     }
 510                   } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
 511                     if( _ProductionMode == 2 ){
 512                       comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
 513                     }else{
 514                       comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
 515                     }
 516                   } else {
 517                     comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
 518                   }
 519
 520                   if( xtest > comp ) goto L203vm;
 521                 }
 522
 523                 /*
 524                 if(_VMCoherence==0 && (!(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2))){
 525                      //dsig/dt= exp(-_VMbslope*t)
 526                      xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 527                      t2 = (-1./_VMbslope)*log(xtest);
 528                      pt2 = sqrt(1.*t2);
 529                 }
 530                 */
 531
 532         }//else end from pp
 533         phi2 = 2.*starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
 534
 535         //        cout<<" pt2: "<<pt2<<endl;
 536
 537         px1 = pt1*cos(phi1);
 538         py1 = pt1*sin(phi1);
 539         px2 = pt2*cos(phi2);
 540         py2 = pt2*sin(phi2);
 541
 542         // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
 543         px = px1 + px2;
 544         py = py1 + py2;
 545         pt = sqrt( px*px + py*py );
 546
 547         E  = sqrt(W*W+pt*pt)*cosh(Y);
 548         pz = sqrt(W*W+pt*pt)*sinh(Y);
 549
 550         // cout<< " Y = "<<Y<<" W = "<<W<<" Egam = "<<Egam<<" gamma = "<<_VMgamma_em<<endl;
 551
 552         // Randomly choose to make pz negative 50% of the time
 553         if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
 554                 pz = -pz;
 555         }else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
 556           // Don't switch
 557         }
 558         else{
 559                 if (randyInstance.Rndom() >= 0.5) pz = -pz;
 560         }
 561
 562 }
 563
 564 //______________________________________________________________________________
 565 double Gammaavectormeson::pTgamma(double E)
 566 {
 567     // returns on random draw from pp(E) distribution
 568     double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
 569     double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
 570     int satisfy =0;
 571
 572     ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
 573     //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
 574     Cm = sqrt(3.)*E/_VMgamma_em;
 575
 576     //the amplitude of the p_t spectrum at the maximum
 577
 578     if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
 579       if( _ProductionMode == 2 ){
 580          singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
 581       }else{
 582          singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
 583       }
 584     } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
 585       if( _ProductionMode == 2 ){
 586          singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
 587       }else{
 588          singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
 589       }
 590     } else {
 591       singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
 592     }
 593
 594     Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
 595
 596     //pick a test value pp, and find the amplitude there
 597     x = randyInstance.Rndom();
 598
 599     if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
 600       if( _ProductionMode == 2 ){
 601         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 602         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
 603       }else{
 604         pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 605         singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
 606       }
 607     } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
 608       if( _ProductionMode == 2 ){
 609         pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 610         singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
 611       }else{
 612         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 613         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
 614       }
 615     } else {
 616         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 617         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
 618     }
 619
 620     test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
 621
 622     while(satisfy==0){
 623         u = randyInstance.Rndom();
 624         if(u*Coef <= test)
 625         {
 626             satisfy =1;
 627         }
 628         else{
 629             x =randyInstance.Rndom();
 630             if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
 631               if( _ProductionMode == 2 ){
 632                 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 633                 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
 634               }else{
 635                 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 636                 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
 637               }
 638             } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
 639               if( _ProductionMode == 2 ){
 640                 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
 641                 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
 642               }else{
 643                 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 644                 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
 645               }
 646             } else {
 647               pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
 648               singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
 649             }
 650             test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
 651         }
 652     }
 653
 654     return pp;
 655 }
 656
 657
 658 //______________________________________________________________________________
 659 void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
 660                              int&) // tcheck (unused)
 661 {
 662         //    This function calculates momentum and energy of vector meson
 663         //     given W and Y, including interference.
 664         //     It gets the pt distribution from a lookup table.
 665         double dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
 666         int IY=0,j=0;
 667
 668         dW = (_VMWmax-_VMWmin)/double(_VMnumw);
 669         dY  = (_VMYmax-_VMYmin)/double(_VMnumy);
 670
 671         //  Y is already fixed; choose a pt
 672         //  Follow the approavh in pickwy.f
 673         // in   _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
 674
 675         IY=int(fabs(Y)/dY);//+1;
 676         if (IY > (_VMnumy/2)-1){
 677                 IY=(_VMnumy/2)-1;
 678         }
 679
 680         yleft=fabs(Y)-(IY)*dY;
 681         yfract=yleft*dY;
 682
 683         xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
 684
 685         for(j=0;j<_VMNPT+1;j++){
 686                 if (xpt < _fptarray[IY][j]) goto L60;
 687         }
 688  L60:
 689
 690         //  now do linear interpolation - start with extremes
 691
 692         if (j == 0){
 693                 pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
 694                 goto L80;
 695         }
 696         if (j == _VMNPT){
 697                 pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
 698                 goto L80;
 699         }
 700
 701         //  we're in the middle
 702
 703         ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
 704         pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
 705
 706         //  at an extreme in y?
 707         if (IY == (_VMnumy/2)-1){
 708                 pt=pt1;
 709                 goto L120;
 710         }
 711  L80:
 712         //  interpolate in y repeat for next fractional y bin
 713
 714         for(j=0;j<_VMNPT+1;j++){
 715                 if (xpt < _fptarray[IY+1][j]) goto L90;
 716         }
 717  L90:
 718
 719         //  now do linear interpolation - start with extremes
 720
 721         if (j == 0){
 722                 pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
 723                 goto L100;
 724         }
 725         if (j == _VMNPT){
 726                 pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
 727                 goto L100;
 728         }
 729
 730         //  we're in the middle
 731
 732         ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
 733         pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
 734
 735  L100:
 736
 737         //  now interpolate in y
 738
 739         pt=yfract*pt2+(1-yfract)*pt1;
 740
 741  L120:
 742
 743         //  we have a pt
 744
 745         theta=2.*starlightConstants::pi*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
 746         px=pt*cos(theta);
 747         py=pt*sin(theta);
 748
 749         //      I guess W is the mass of the vector meson (not necessarily
 750         //      on-mass-shell), and E is the energy
 751
 752         E  = sqrt(W*W+pt*pt)*cosh(Y);
 753         pz = sqrt(W*W+pt*pt)*sinh(Y);
 754         //      randomly choose to make pz negative 50% of the time
 755         if(randyInstance.Rndom()>=0.5) pz = -pz;
 756 }
 757
 758
 759 //______________________________________________________________________________
 760 starlightConstants::event Gammaavectormeson::produceEvent(int&)
 761 {
 762         // Not used; return default event
 763         return starlightConstants::event();
 764 }
 765
 766
 767 //______________________________________________________________________________
 768 upcEvent Gammaavectormeson::produceEvent()
 769 {
 770         // The new event type
 771         upcEvent event;
 772
 773         int iFbadevent=0;
 774         int tcheck=0;
 775         starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
 776         starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
 777
 778         if (_VMpidtest == starlightConstants::FOURPRONG) {
 779                 double        comenergy = 0;
 780                 double        mom[3]    = {0, 0, 0};
 781                 double        E         = 0;
 782                 lorentzVector decayVecs[4];
 783                 do {
 784                         double rapidity = 0;
 785                         pickwy(comenergy, rapidity);
 786                         if (_VMinterferencemode == 0)
 787                                 momenta(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
 788                         else if (_VMinterferencemode==1)
 789                                 vmpt(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
 790                 } while (!fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent));
 791                 if ((iFbadevent == 0) and (tcheck == 0))
 792                         for (unsigned int i = 0; i < 4; ++i) {
 793                                 starlightParticle daughter(decayVecs[i].GetPx(),
 794                                                            decayVecs[i].GetPy(),
 795                                                            decayVecs[i].GetPz(),
 796                                                            starlightConstants::UNKNOWN,  // energy
 797                                                            starlightConstants::UNKNOWN,  // _mass
 798                                                            ipid,
 799                                                            (i < 2) ? -1 : +1);
 800                                 event.addParticle(daughter);
 801                         }
 802         } else {
 803                 double comenergy = 0.;
 804                 double rapidity = 0.;
 805                 double E = 0.;
 806                 double momx=0.,momy=0.,momz=0.;
 807
 808                 double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
 809                 bool accepted = false;
 810                 //  if(_accCut){
 811                 do{
 812                         pickwy(comenergy,rapidity);
 813
 814                         if (_VMinterferencemode==0){
 815                                 momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
 816                         } else if (_VMinterferencemode==1){
 817                                 vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck);
 818                         }
 819
 820                         // cout << "_ptCutMin: " << _ptCutMin << " _ptCutMax: " << _ptCutMax << " _etaCutMin: " << _etaCutMin << " _etaCutMax: " << _etaCutMax << endl;
 821                         _nmbAttempts++;
 822                         //cout << "n tries: " << _nmbAttempts<< endl;
 823                         vmpid = ipid;
 824                         twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 825                         double pt1chk = sqrt(px1*px1+py1*py1);
 826                         double pt2chk = sqrt(px2*px2+py2*py2);
 827
 828                         //cout << "pt1: " << pt1chk  << " pt2: " << pt2chk << endl;
 829                         double eta1 = pseudoRapidity(px1, py1, pz1);
 830                         double eta2 = pseudoRapidity(px2, py2, pz2);
 831                         //cout << "eta1: " << eta1 << " eta2: " << eta2 << endl;
 832                         if(_ptCutEnabled && !_etaCutEnabled){
 833                                 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
 834                                         accepted = true;
 835                                         _nmbAccepted++;
 836                                 }
 837                         }
 838                         else if(!_ptCutEnabled && _etaCutEnabled){
 839                                 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
 840                                         accepted = true;
 841                                         _nmbAccepted++;
 842                                 }
 843                         }
 844                         else if(_ptCutEnabled && _etaCutEnabled){
 845                                 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
 846                                         if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
 847                                                 accepted = true;
 848                                                 _nmbAccepted++;
 849                                         }
 850                                 }
 851                         }
 852                         else if(!_ptCutEnabled && !_etaCutEnabled)
 853                                 _nmbAccepted++;
 854
 855                 }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
 856                 /*  }else{
 857                     twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 858                     }*/
 859                 if (iFbadevent==0&&tcheck==0) {
 860                         int q1=0,q2=0;
 861                         int ipid1,ipid2=0;
 862
 863                         double xtest = randyInstance.Rndom();
 864                         if (xtest<0.5)
 865                                 {
 866                                         q1=1;
 867                                         q2=-1;
 868                                 }
 869                         else {
 870                                 q1=-1;
 871                                 q2=1;
 872                         }
 873
 874                         if ( ipid == 11 || ipid == 13 ){
 875                           ipid1 = -q1*ipid;
 876                           ipid2 = -q2*ipid;
 877                         } else {
 878                           ipid1 = q1*ipid;
 879                           ipid2 = q2*ipid;
 880                         }
 881                         //     The new stuff
 882                         double md = getDaughterMass(vmpid);
 883                         double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1);
 884                         starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
 885                         event.addParticle(particle1);
 886
 887                         double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2);
 888                         starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
 889                         event.addParticle(particle2);
 890                         //     End of the new stuff
 891
 892                 }
 893         }
 894
 895         return event;
 896
 897 }
 898 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
 899 {
 900         double pT = sqrt(px*px + py*py);
 901         double p = sqrt(pz*pz + pT*pT);
 902         double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
 903         return eta;
 904 }
 905
 906 //______________________________________________________________________________
 907 Gammaanarrowvm::Gammaanarrowvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
 908 {
 909         cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
 910         read();
 911         cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
 912         narrowResonanceCrossSection sigma(bbsystem);
 913         sigma.crossSectionCalculation(_bwnormsave);
 914         _VMbslope=sigma.slopeParameter();
 915 }
 916
 917
 918 //______________________________________________________________________________
 919 Gammaanarrowvm::~Gammaanarrowvm()
 920 { }
 921
 922
 923 //______________________________________________________________________________
 924 Gammaaincoherentvm::Gammaaincoherentvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
 925 {
 926         cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
 927         read();
 928         cout<<"Creating and calculating crosssection. Gammaainkoherentvm()"<<endl;
 929         incoherentVMCrossSection sigma(bbsystem); sigma.crossSectionCalculation(_bwnormsave);
 930         _VMbslope=sigma.slopeParameter();
 931 }
 932
 933
 934 //______________________________________________________________________________
 935 Gammaaincoherentvm::~Gammaaincoherentvm()
 936 { }
 937
 938
 939 //______________________________________________________________________________
 940 Gammaawidevm::Gammaawidevm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
 941 {
 942         cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
 943         read();
 944         cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
 945         wideResonanceCrossSection sigma(bbsystem);
 946         sigma.crossSectionCalculation(_bwnormsave);
 947         _VMbslope=sigma.slopeParameter();
 948 }
 949
 950
 951 //______________________________________________________________________________
 952 Gammaawidevm::~Gammaawidevm()
 953 { }
 954
 955