EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
gammaavm.cpp
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file gammaavm.cpp
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 //
21 //
22 // File and Version Information:
23 // $Rev:: 277 $: revision of last commit
24 // $Author:: jnystrand $: author of last commit
25 // $Date:: 2016-09-14 10:55:55 +0100 #$: date of last commit
26 //
27 // Description:
28 // Added incoherent t2-> pt2 selection. Following pp selection scheme
29 //
30 //
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 #include "e_narrowResonanceCrossSection.h"
46 #include "e_wideResonanceCrossSection.h"
47 
48 using namespace std;
49 
50 
51 //______________________________________________________________________________
52 Gammaavectormeson::Gammaavectormeson(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem):eventChannel(inputParametersInstance, bbsystem), _phaseSpaceGen(0)
53 {
54  _VMNPT=inputParametersInstance.nmbPtBinsInterference();
55  _VMWmax=inputParametersInstance.maxW();
56  _VMWmin=inputParametersInstance.minW();
57  //_VMYmax=inputParametersInstance.maxRapidity();
58  //_VMYmin=-1.*_VMYmax;
59  _VMnumw=inputParametersInstance.nmbWBins();
60  //_VMnumy=inputParametersInstance.nmbRapidityBins();
61  _VMnumega=inputParametersInstance.nmbEnergyBins();
62  _VMnumQ2=inputParametersInstance.nmbGammaQ2Bins();
63  _VMgamma_em=inputParametersInstance.beamLorentzGamma();
64  _VMinterferencemode=inputParametersInstance.interferenceEnabled();
65  _VMbslope=0.;//Will define in wide/narrow constructor
66  _bslopeDef=inputParametersInstance.bslopeDefinition();
67  _bslopeVal=inputParametersInstance.bslopeValue();
68  _pEnergy= inputParametersInstance.protonEnergy();
69  _beamNucleus = inputParametersInstance.targetBeamA();
70  // electron energy in CMS frame
71  _eEnergy= inputParametersInstance.electronEnergy();
72  _VMpidtest=inputParametersInstance.prodParticleType();
73  _VMptmax=inputParametersInstance.maxPtInterference();
74  _VMdpt=inputParametersInstance.ptBinWidthInterference();
75  _ProductionMode=inputParametersInstance.productionMode();
76  _targetMaxPhotonEnergy=inputParametersInstance.targetMaxPhotonEnergy();
77  _targetMinPhotonEnergy=inputParametersInstance.targetMinPhotonEnergy();
78  // Now saving the photon energy limits
79  _cmsMaxPhotonEnergy=inputParametersInstance.cmsMaxPhotonEnergy();
80  _cmsMinPhotonEnergy=inputParametersInstance.cmsMinPhotonEnergy();
81  _beamLorentzGamma = inputParametersInstance.beamLorentzGamma();
82  _targetBeamLorentzGamma = inputParametersInstance.targetBeamLorentzGamma();
83  _rap_CM=inputParametersInstance.rap_CM();
84  _targetRadius = inputParametersInstance.targetRadius();
85  //Turn on/off backwards production
86  _backwardsProduction = inputParametersInstance.backwardsProduction();
87 
88  N0 = 0; N1 = 0; N2 = 0;
89  if (_VMpidtest == starlightConstants::FOURPRONG){
90  // create n-body phase-spage generator
91  _phaseSpaceGen = new nBodyPhaseSpaceGen(_randy);
92  }
93  if(_ProductionMode == 12 || _ProductionMode == 13) // Narrow and wide coherent photon-pommeron in eSTARlight
94  _dummy_pncs = new photonNucleusCrossSection(inputParametersInstance, bbsystem);
95  //
96  const double r_04_00 = 0.674;
97  const double cos_delta = 0.925;
98  for( int ii = 0; ii < 100; ++ii){ //epsilon 0-1
99  double epsilon = 0.01*ii;
100  const double R = (1./epsilon)*r_04_00/(1.-r_04_00);
101  for(int jj = 0; jj < 200; ++jj){ //psi 0 - 2pi
102  double psi = jj*starlightConstants::pi/100.;
103  double max_bin;
104  for( int kk = 0; kk < 200; ++kk){ //temp
105  double theta = kk*starlightConstants::pi/100.;
106  // Fin max
107  double this_test = std::pow(std::sin(theta),2.)*(1+epsilon*cos(2.*psi)) + 2.*epsilon*R*std::pow(std::cos(theta),2.)
108  +std::sqrt(2.*epsilon*(1+epsilon))*cos_delta*std::sin(2.*theta)*std::cos(psi);
109  if(this_test > max_bin)
110  max_bin = this_test;
111  }
112  _angular_max[ii][jj] = max_bin;
113  }
114  }
115 }
116 
117 
118 //______________________________________________________________________________
119 Gammaavectormeson::~Gammaavectormeson()
120 {
121  if (_phaseSpaceGen)
122  delete _phaseSpaceGen;
123  if (_dummy_pncs)
124  delete _dummy_pncs;
125 }
126 
127 
128 //______________________________________________________________________________
129 void Gammaavectormeson::pickwy(double &W, double &Y)
130 {
131  double dW, dY, xw,xy,xtest;
132  int IW,IY;
133 
134  dW = (_VMWmax-_VMWmin)/double(_VMnumw);
135  dY = (_VMYmax-_VMYmin)/double(_VMnumy);
136 
137  L201pwy:
138 
139  xw = _randy.Rndom();
140  W = _VMWmin + xw*(_VMWmax-_VMWmin);
141 
142  if (W < 2 * starlightConstants::pionChargedMass)
143  goto L201pwy;
144 
145  IW = int((W-_VMWmin)/dW);
146  xy = _randy.Rndom();
147  Y = _VMYmin + xy*(_VMYmax-_VMYmin);
148  IY = int((Y-_VMYmin)/dY);
149  xtest = _randy.Rndom();
150 
151  if( xtest > _Farray[IW][IY] )
152  goto L201pwy;
153 
154  N0++;
155  // Determine the target nucleus
156  // For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2
157  _TargetBeam=2; // Always true for eX
158 }
159 
160 
161 
162 //______________________________________________________________________________
163 void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
164  double W,
165  double px0, double py0, double pz0,
166  double spin_element,
167  double& px1, double& py1, double& pz1,
168  double& px2, double& py2, double& pz2,
169  int& iFbadevent)
170 {
171  // This routine decays a particle into two particles of mass mdec,
172  // taking spin into account
173  double pmag;
174  double phi,theta,Ecm;
175  double betax,betay,betaz;
176  double mdec1=0.0,mdec2=0.0;
177  double E1=0.0,E2=0.0;
178 
179  // set the mass of the daughter particles
180  mdec1=getDaughterMass(ipid);
181 
182  if(_VMpidtest == starlightConstants::OMEGA_pi0gamma){
183  mdec2=starlightConstants::pionNeutralMass;
184  if(W < mdec2){
185  cout<<" ERROR: W="<<W<<endl;
186  iFbadevent = 1;
187  return;
188  }
189  // calculate the magnitude of the momenta
190  pmag = (W*W - mdec2*mdec2)/(2.0*W);
191  }else{
192  mdec2=mdec1;
193  if(W < 2*mdec1){
194  cout<<" ERROR: W="<<W<<endl;
195  iFbadevent = 1;
196  return;
197  }
198  // calculate the magnitude of the momenta
199  pmag = sqrt(W*W/4. - mdec2*mdec2);
200  }
201 
202  // pick an orientation, based on the spin
203  // phi has a flat distribution in 2*pi
204  phi = _randy.Rndom()*2.*starlightConstants::pi;
205 
206  // find theta, the angle between one of the outgoing particles and
207  // the beamline, in the frame of the two photons
208  //cout<<"spin element: "<<spin_element<<endl;
209  theta=getTheta(ipid, spin_element);
210 
211  // compute unboosted momenta
212  px1 = sin(theta)*cos(phi)*pmag;
213  py1 = sin(theta)*sin(phi)*pmag;
214  pz1 = cos(theta)*pmag;
215  px2 = -px1;
216  py2 = -py1;
217  pz2 = -pz1;
218 
219  Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
220  E1 = sqrt(mdec1*mdec1+px1*px1+py1*py1+pz1*pz1);
221  E2 = sqrt(mdec2*mdec2+px2*px2+py2*py2+pz2*pz2);
222 
223  //cout<<"Ecm: "<<Ecm<<endl;
224  //cout<<"W: "<<W<<endl;
225 
226  betax = -(px0/Ecm);
227  betay = -(py0/Ecm);
228  betaz = -(pz0/Ecm);
229 
230  transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
231  transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
232 
233  if(iFbadevent == 1)
234  return;
235 
236 }
237 
238 
239 //______________________________________________________________________________
240 // decays a particle into four particles with isotropic angular distribution
241 bool Gammaavectormeson::fourBodyDecay
242 (starlightConstants::particleTypeEnum& ipid,
243  const double , // E (unused)
244  const double W, // mass of produced particle
245  const double* p, // momentum of produced particle; expected to have size 3
246  lorentzVector* decayVecs, // array of Lorentz vectors of daughter particles; expected to have size 4
247  int& iFbadevent)
248 {
249  const double parentMass = W;
250 
251  // set the mass of the daughter particles
252  const double daughterMass = getDaughterMass(ipid);
253  if (parentMass < 4 * daughterMass){
254  cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
255  iFbadevent = 1;
256  return false;
257  }
258 
259  // construct parent four-vector
260  const double parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
261  + parentMass * parentMass);
262  const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
263 
264  // setup n-body phase-space generator
265  assert(_phaseSpaceGen);
266  static bool firstCall = true;
267  if (firstCall) {
268  const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
269  _phaseSpaceGen->setDecay(4, m);
270  // estimate maximum phase-space weight
271  _phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
272  firstCall = false;
273  }
274 
275  // generate phase-space event
276  if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
277  return false;
278 
279  // set Lorentzvectors of decay daughters
280  for (unsigned int i = 0; i < 4; ++i)
281  decayVecs[i] = _phaseSpaceGen->daughter(i);
282  return true;
283 }
284 
285 //______________________________________________________________________________
286 void Gammaavectormeson::pi0Decay(double& px_pi0, double& py_pi0, double& pz_pi0,
287  double& e_g1, double& px_g1, double& py_g1, double& pz_g1,
288  double& e_g2, double& px_g2, double& py_g2, double& pz_g2,
289  int& iFbadevent)
290 {
291  double pmag;
292  double phi,theta,Ecm;
293  double betax,betay,betaz;
294  double E1=0.0,E2=0.0;
295 
296  // This routine decays a pi0 into two isotropically produced photons
297  double m_pi0=starlightConstants::pionNeutralMass;
298 
299  // calculate the magnitude of the momenta
300  pmag = sqrt(m_pi0*m_pi0/4.);
301 
302  // pick an orientation, based on the spin
303  // phi has a flat distribution in 2*pi
304  phi = _randy.Rndom()*2.*starlightConstants::pi;
305 
306  // find theta, the angle between one of the outgoing photons and
307  // the beamline, in the frame of the pi0
308  theta=getTheta(starlightConstants::PION, 1.0/3.0);
309 
310  // compute unboosted momenta
311  px_g1 = sin(theta)*cos(phi)*pmag;
312  py_g1 = sin(theta)*sin(phi)*pmag;
313  pz_g1 = cos(theta)*pmag;
314  px_g2 = -px_g1;
315  py_g2 = -py_g1;
316  pz_g2 = -pz_g1;
317 
318  Ecm = sqrt(m_pi0*m_pi0+px_pi0*px_pi0+py_pi0*py_pi0+pz_pi0*pz_pi0);
319  E1 = sqrt(px_g1*px_g1+py_g1*py_g1+pz_g1*pz_g1);
320  E2 = sqrt(px_g2*px_g2+py_g2*py_g2+pz_g2*pz_g2);
321 
322  betax = -(px_pi0/Ecm);
323  betay = -(py_pi0/Ecm);
324  betaz = -(pz_pi0/Ecm);
325 
326  transform (betax,betay,betaz,E1,px_g1,py_g1,pz_g1,iFbadevent);
327  transform (betax,betay,betaz,E2,px_g2,py_g2,pz_g2,iFbadevent);
328 
329  e_g1=E1;
330  e_g2=E2;
331 
332  if(iFbadevent == 1)
333  return;
334 }
335 
336 
337 //______________________________________________________________________________
338 double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
339 {
340  //This will return the daughter particles mass, and the final particles outputed id...
341  double mdec=0.;
342 
343  switch(_VMpidtest){
344  case starlightConstants::RHO:
345  case starlightConstants::RHOZEUS:
346  case starlightConstants::FOURPRONG:
347  case starlightConstants::OMEGA:
348  mdec = starlightConstants::pionChargedMass;
349  ipid = starlightConstants::PION;
350  break;
351  case starlightConstants::OMEGA_pi0gamma:
352  mdec = 0.0;
353  ipid = starlightConstants::PHOTON;
354  break;
355  case starlightConstants::PHI:
356  mdec = starlightConstants::kaonChargedMass;
357  ipid = starlightConstants::KAONCHARGE;
358  break;
359  case starlightConstants::JPSI:
360  mdec = starlightConstants::mel;
361  ipid = starlightConstants::ELECTRON;
362  break;
363  case starlightConstants::JPSI_ee:
364  mdec = starlightConstants::mel;
365  ipid = starlightConstants::ELECTRON;
366  break;
367  case starlightConstants::JPSI_mumu:
368  mdec = starlightConstants::muonMass;
369  ipid = starlightConstants::MUON;
370  break;
371  case starlightConstants::JPSI2S_ee:
372  mdec = starlightConstants::mel;
373  ipid = starlightConstants::ELECTRON;
374  break;
375  case starlightConstants::JPSI2S_mumu:
376  mdec = starlightConstants::muonMass;
377  ipid = starlightConstants::MUON;
378  break;
379 
380  case starlightConstants::JPSI2S:
381  case starlightConstants::UPSILON:
382  case starlightConstants::UPSILON2S:
383  case starlightConstants::UPSILON3S:
384  mdec = starlightConstants::muonMass;
385  ipid = starlightConstants::MUON;
386  break;
387  case starlightConstants::UPSILON_ee:
388  case starlightConstants::UPSILON2S_ee:
389  case starlightConstants::UPSILON3S_ee:
390  mdec = starlightConstants::mel;
391  ipid = starlightConstants::ELECTRON;
392  break;
393  case starlightConstants::UPSILON_mumu:
394  case starlightConstants::UPSILON2S_mumu:
395  case starlightConstants::UPSILON3S_mumu:
396  mdec = starlightConstants::muonMass;
397  ipid = starlightConstants::MUON;
398  break;
399  default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
400  }
401 
402  return mdec;
403 }
404 
405 
406 
407 //______________________________________________________________________________
408 double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid, double r_04_00)
409 {
410  //This depends on the decay angular distribution
411  //Valid for rho, phi, omega.
412  double theta=0.;
413  double xtest=1.;
414  double dndtheta=0.;
415 
416  while(xtest > dndtheta){
417 
418  theta = starlightConstants::pi*_randy.Rndom();
419  xtest = _randy.Rndom();
420  // Follow distribution for helicity +/-1 with finite Q2
421  // Physics Letters B 449, 328; The European Physical Journal C - Particles and Fields 13, 37;
422  // The European Physical Journal C - Particles and Fields 6, 603
423 
424 
425  switch(ipid){
426 
427  case starlightConstants::MUON:
428  case starlightConstants::ELECTRON:
429  //primarily for upsilon/j/psi. VM->ee/mumu
430  dndtheta = sin(theta)*(1 + r_04_00+( 1-3.*r_04_00 )*cos(theta)*cos(theta));
431  break;
432 
433  case starlightConstants::PION:
434  case starlightConstants::KAONCHARGE:
435  //rhos etc
436  dndtheta= sin(theta)*(1 - r_04_00+( 3.*r_04_00-1 )*cos(theta)*cos(theta));
437  break;
438 
439  default: if(!_backwardsProduction) cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
440  }//end of switch
441 
442  // Assume unpolarized vector-mesons for now
443  if(_backwardsProduction){
444  dndtheta = sin(theta);
445  }
446  }
447  return theta;
448 
449 }
450 
451 
452 //______________________________________________________________________________
453 double Gammaavectormeson::getSpinMatrixElement(double W, double Q2, double epsilon)
454 {
455  double m2 = 0.6*(W*W);
456  double R = starlightConstants::pi/2.*m2/Q2 - std::pow(m2,3./2.)/sqrt(Q2)/(Q2+m2) - m2/Q2*atan(sqrt(m2/Q2));
458  double R1 = starlightConstants::pi/2.*m2/Q2;
459  if(R1>1.0E10)R=0.0;
461  R = (Q2+m2)*(Q2+m2)*R*R/m2;
462  return epsilon*R/(1+epsilon*R);
463 }
464 
465 
466 //______________________________________________________________________________
467 double Gammaavectormeson::getWidth()
468 {
469  return _width;
470 }
471 
472 
473 //______________________________________________________________________________
474 double Gammaavectormeson::getMass()
475 {
476  return _mass;
477 }
478 
479 
480 //______________________________________________________________________________
481 double Gammaavectormeson::getSpin()
482 {
483  return 1.0; //VM spins are the same
484 }
485 
486 //______________________________________________________________________________
487 void Gammaavectormeson::momenta(double W,double Egam,double Q2, double gamma_pz, double gamma_pt, //input conditions
488  double &Y,double &E,double &px,double &py,double &pz, //return vm
489  double &t_px, double &t_py, double &t_pz, double &t_E, //return pomeron
490  double &e_phi,int &tcheck) //return electron (angle already known by Q2)
491 {
492  // This subroutine calculates momentum and energy of vector meson for electroproduction (eSTARlight)
493  // given W and photon 4-vector, without interference. No intereference in asymetric eX collisions
494 
495  double Epom,Pom_pz,tmin,pt2,phi1,phi2;
496  double px1,py1;
497  double xt,xtest,ytest;
498  double t2;
499 
500  double target_px, target_py, target_pz, target_E;
501 
502  target_E = _beamNucleus*_pEnergy;
503  target_px = 0.0;
504  target_py = 0.0;
505  target_pz = -_beamNucleus*sqrt(_pEnergy*_pEnergy - pow(starlightConstants::protonMass,2.) );
506  phi1 = 2.*starlightConstants::pi*_randy.Rndom();
507  px1 = gamma_pt*cos(phi1);
508  py1 = gamma_pt*sin(phi1);
509  int isbadevent = 0;
510  double betax_cm = ((px1+target_px)/(Egam+target_E));
511  double betay_cm = ((py1+target_py)/(Egam+target_E));
512  double betaz_cm = ((gamma_pz+target_pz)/(Egam+target_E));
513  transform (betax_cm,betay_cm,betaz_cm,target_E,target_px,target_py,target_pz,isbadevent);
514  transform (betax_cm,betay_cm,betaz_cm,Egam, px1, py1, gamma_pz, isbadevent);
515 
516  e_phi = starlightConstants::pi+phi1;
517  // Pomeron pz is != than its energy in eSTARlight, in order to conserve energy/momentum of scattered
518  // target
519  //Pom_pz = 0.5*(W*W-Q2)/(Egam + gamma_pz);
520  Epom = 0.5*(W*W+Q2)/(Egam + target_E);
521 
522  L522vm:
523  while( e_phi > 2.*starlightConstants::pi ) e_phi-= 2.*starlightConstants::pi;
524  //
525  if ( ( _bbs.targetBeam().A()==1 )
526  || (_ProductionMode == 4) ) {
527  if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
528  // Use dipole form factor for light VM
529  L613vm:
530  xtest = 2.0*_randy.Rndom();
531  double ttest = xtest*xtest;
532  ytest = _randy.Rndom();
533  double t0 = 1./2.23;
534  double yprob = xtest*_bbs.electronBeam().dipoleFormFactor(ttest,t0)*_bbs.electronBeam().dipoleFormFactor(ttest,t0);
535  if( ytest > yprob ) goto L613vm;
536  t2 = ttest;
537  pt2 = xtest;
538  }else{
539  //Use dsig/dt= exp(-_VMbslope*t) for heavy VM
540  double bslope_tdist = _VMbslope;
541  double Wgammap = 0.0;
542  switch(_bslopeDef){
543  case 0:
544  //This is the default, as before
545  bslope_tdist = _VMbslope;
546  break;
547  case 1:
548  //User defined value of bslope. BSLOPE_VALUE default is 4.0 if not set.
549  bslope_tdist = _bslopeVal;
550  if( N0 <= 1 )cout<<" ATTENTION: Using user defined value of bslope = "<<_bslopeVal<<endl;
551  break;
552  case 2:
553  //This is Wgammap dependence of b from H1 (Eur. Phys. J. C 46 (2006) 585)
554  Wgammap = sqrt(4.*Egam*_pEnergy);
555  bslope_tdist = 4.63 + 4.*0.164*log(Wgammap/90.0);
556  if( N0 <= 1 )cout<<" ATTENTION: Using energy dependent value of bslope!"<<endl;
557  break;
558  default:
559  cout<<" Undefined setting for BSLOPE_DEFINITION "<<endl;
560  }
561 
562  xtest = _randy.Rndom();
563  // t2 = (-1./_VMbslope)*log(xtest);
564  t2 = (-1./bslope_tdist)*log(xtest);
565  pt2 = sqrt(1.*t2);
566  }
567  } else {
568  // >> Check tmin
569  tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
570 
571  if(tmin > 0.5){
572  cout<<" WARNING: tmin= "<<tmin<<endl;
573  cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl;
574  cout<<" Will pick a new W,Y "<<endl;
575  tcheck = 1;
576  return;
577  }
578  L663vm:
579  xt = _randy.Rndom();
580  if( _bbs.targetBeam().A() != 1){
581  pt2 = 8.*xt*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
582  }
583  else{
584  std::cout<<"Can't find the electron for eX"<<std::endl;
585  }
586 
587  xtest = _randy.Rndom();
588  t2 = tmin + pt2*pt2;
589 
590  double comp=0.0;
591  if( _bbs.targetBeam().A() != 1){
592  comp = _bbs.targetBeam().formFactor(t2)*_bbs.targetBeam().formFactor(t2)*pt2;
593  }
594  else
595  std::cout<<"Can't find the electron for eX"<<std::endl;
596  if( xtest > comp ) goto L663vm;
597 
598  }
599  phi2 = 2.*starlightConstants::pi*_randy.Rndom();
600 
601  //
602  t_px = pt2*cos(phi2);
603  t_py = pt2*sin(phi2);
604  // Used to return the pomeron pz to generator
605  // Compute scattered target kinematics p_(initial ion) = p_(pomeron) + p_(final ion)
606  double newion_E = target_E-Epom;
607  double newion_px = target_px - t_px;
608  double newion_py = target_py - t_py;
609  double newion_pz_squared = newion_E*newion_E - newion_px*newion_px - newion_py*newion_py - pow(_beamNucleus*starlightConstants::protonMass,2.);
610  if(newion_pz_squared<0)goto L522vm;
611  double newion_pz = -sqrt( newion_pz_squared );
612  Pom_pz = target_pz - newion_pz;
613  t_pz = Pom_pz;
614  // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
615  px = px1 + t_px;
616  py = py1 + t_py;
617 
618  t_E = Epom;
619  // Finally V.M. energy, pz and rapidity from photon + pommeron.
620  E = Egam + t_E;
621  pz = gamma_pz + t_pz;
622  //transform back to e-ion CM frame
623  transform (-betax_cm,-betay_cm,-betaz_cm,target_E,target_px,target_py,target_pz,isbadevent);
624  transform (-betax_cm,-betay_cm,-betaz_cm,newion_E,newion_px,newion_py,newion_pz,isbadevent);
625  transform (-betax_cm,-betay_cm,-betaz_cm,Egam, px1, py1, gamma_pz, isbadevent);
626  transform (-betax_cm,-betay_cm,-betaz_cm,E, px, py, pz, isbadevent);
627  t_px = target_px-newion_px;
628  t_py = target_py-newion_py;
629  t_pz = target_pz-newion_pz;
630  t_E = target_E-newion_E;
631  Y = 0.5*std::log( (E+fabs(pz))/(E-fabs(pz)) );
632  // cout << "pz_final = " << pz << endl;
633 
634 
635  // Backwards Production... updated by Zachary Sweger on 9/21/2021
636  if (_backwardsProduction) {
637  double is_proton_E = target_E;
638  double is_proton_px = target_px;
639  double is_proton_py = target_py;
640  double is_proton_pz = target_pz;
641  double is_gamma_E = Egam;
642  double is_gamma_px = px1;
643  double is_gamma_py = py1;
644  double is_gamma_pz = gamma_pz;
645  double is_tot_E = is_proton_E + is_gamma_E;
646  double is_tot_px = is_proton_px + is_gamma_px;
647  double is_tot_py = is_proton_py + is_gamma_py;
648  double is_tot_pz = is_proton_pz + is_gamma_pz;
649 
650  double fs_proton_E = newion_E;
651  double fs_proton_px = newion_px;
652  double fs_proton_py = newion_py;
653  double fs_proton_pz = newion_pz;
654 
655  double fs_vm_E = E;
656  double fs_vm_px = px;
657  double fs_vm_py = py;
658  double fs_vm_pz = pz;
659 
660  int isbadevent = 0;
661  double betax_cm = (is_tot_px/is_tot_E);
662  double betay_cm = (is_tot_py/is_tot_E);
663  double betaz_cm = (is_tot_pz/is_tot_E);
664  transform (betax_cm,betay_cm,betaz_cm,is_proton_E,is_proton_px,is_proton_py,is_proton_pz,isbadevent);
665  transform (betax_cm,betay_cm,betaz_cm,is_gamma_E, is_gamma_px, is_gamma_py, is_gamma_pz, isbadevent);
666  transform (betax_cm,betay_cm,betaz_cm,fs_proton_E,fs_proton_px,fs_proton_py,fs_proton_pz,isbadevent);
667  transform (betax_cm,betay_cm,betaz_cm,fs_vm_E, fs_vm_px, fs_vm_py, fs_vm_pz, isbadevent);
668 
669  double u_fs_proton_px = -1.0*fs_proton_px;
670  double u_fs_proton_py = -1.0*fs_proton_py;
671  double u_fs_proton_pz = -1.0*fs_proton_pz;
672  double u_fs_proton_E = fs_proton_E;
673 
674  double u_fs_vm_px = -1.0*fs_vm_px;
675  double u_fs_vm_py = -1.0*fs_vm_py;
676  double u_fs_vm_pz = -1.0*fs_vm_pz;
677  double u_fs_vm_E = fs_vm_E;
678 
679  betax_cm = -1.0*betax_cm;
680  betay_cm = -1.0*betay_cm;
681  betaz_cm = -1.0*betaz_cm;
682  transform (betax_cm,betay_cm,betaz_cm,is_proton_E, is_proton_px, is_proton_py, is_proton_pz, isbadevent);
683  transform (betax_cm,betay_cm,betaz_cm,is_gamma_E, is_gamma_px, is_gamma_py, is_gamma_pz, isbadevent);
684  transform (betax_cm,betay_cm,betaz_cm,u_fs_proton_E,u_fs_proton_px,u_fs_proton_py,u_fs_proton_pz,isbadevent);
685  transform (betax_cm,betay_cm,betaz_cm,u_fs_vm_E, u_fs_vm_px, u_fs_vm_py, u_fs_vm_pz, isbadevent);
686 
687  px=u_fs_vm_px;
688  py=u_fs_vm_py;
689  pz=u_fs_vm_pz;
690  E =u_fs_vm_E;
691 
692  t_px=-u_fs_proton_px;
693  t_py=-u_fs_proton_py;
694  t_pz=is_proton_pz-u_fs_proton_pz;
695  t_E =is_proton_E-u_fs_proton_E;
696  //cout<<"Energy Violation: "<<(E+u_fs_proton_E)-(is_tot_E)<<endl;
697  //cout<<"Pz Violation: "<<(pz+u_fs_proton_pz)-(is_tot_pz)<<endl;
698  //cout<<"Px Violation: "<<(px+u_fs_proton_px)-(is_tot_px)<<endl;
699  //cout<<"Py Violation: "<<(py+u_fs_proton_py)-(is_tot_py)<<endl;
700  //cout<<"proton mass: "<< sqrt(u_fs_proton_E*u_fs_proton_E-u_fs_proton_pz*u_fs_proton_pz-u_fs_proton_py*u_fs_proton_py-u_fs_proton_px*u_fs_proton_px)<<endl;
701  //cout << "VM MASS = " << sqrt(E*E-px*px - py*py-pz*pz) << endl;
702 
703  //Calculate the rapidity
704  Y = 0.5*std::log( (E+fabs(pz))/(E-fabs(pz)) );
705  if(Y!=Y){
706  //FIXME the following should be upated if not using omegas
707  if(_VMpidtest==starlightConstants::OMEGA_pi0gamma || _VMpidtest==starlightConstants::OMEGA) E=sqrt(starlightConstants::OmegaMass*starlightConstants::OmegaMass+px*px+py*py+pz*pz);
708  else if(_VMpidtest==starlightConstants::RHO) E=sqrt(starlightConstants::rho0Mass*starlightConstants::rho0Mass+px*px+py*py+pz*pz);
709  //cout<<"Energy Violation: "<<(E+u_fs_proton_E)-(is_tot_E)<<endl;
710  }
711  }
712 
713 }
714 
715 
716 //______________________________________________________________________________
717 double Gammaavectormeson::pTgamma(double E)
718 {
719  // returns on random draw from pp(E) distribution
720  double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
721  double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
722  int satisfy =0;
723 
724  ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
725  //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
726  Cm = sqrt(3.)*E/_VMgamma_em;
727  // If E is very small, the drawing of a pT below is extremely slow.
728  // ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E.
729  // Should have no observable consequences (JN, SRK 11-Sep-2014)
730  if( E < 0.0005 )return Cm;
731 
732  //the amplitude of the p_t spectrum at the maximum
733 
734  if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A() != 1){
735  if( _ProductionMode == 2 || _ProductionMode ==3 ){
736  singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
737  }else{
738  singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
739  }
740  } else if( _bbs.targetBeam().A()==1 && _bbs.electronBeam().A() != 1 ){
741  if( _ProductionMode == 2 || _ProductionMode ==3){
742  singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
743  }else{
744  singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
745  }
746  } else if (_TargetBeam == 1) {
747  singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
748  } else {
749  singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
750  }
751 
752  Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
753 
754  //pick a test value pp, and find the amplitude there
755  x = _randy.Rndom();
756 
757  if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A() != 1){
758  if( _ProductionMode == 2 || _ProductionMode ==3){
759  pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius();
760  singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
761  }else{
762  pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
763  singleformfactorpp1=_bbs.targetBeam().formFactor(pp*pp+ereds);
764  }
765  } else if( _bbs.targetBeam().A()==1 && _bbs.electronBeam().A() != 1 ){
766  if( _ProductionMode == 2 || _ProductionMode ==3){
767  pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
768  singleformfactorpp1=_bbs.targetBeam().formFactor(pp*pp+ereds);
769  }else{
770  pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius();
771  singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
772  }
773  } else if (_TargetBeam == 1) {
774  pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
775  singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
776  } else {
777  pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius();
778  singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
779  }
780 
781  test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
782 
783  while(satisfy==0){
784  u = _randy.Rndom();
785  if(u*Coef <= test)
786  {
787  satisfy =1;
788  }
789  else{
790  x =_randy.Rndom();
791  if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A() != 1){
792  if( _ProductionMode == 2 || _ProductionMode ==3){
793  pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius();
794  singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
795  }else{
796  pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
797  singleformfactorpp2=_bbs.targetBeam().formFactor(pp*pp+ereds);
798  }
799  } else if( _bbs.targetBeam().A()==1 && _bbs.electronBeam().A() != 1 ){
800  if( _ProductionMode == 2 || _ProductionMode ==3){
801  pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
802  singleformfactorpp2=_bbs.targetBeam().formFactor(pp*pp+ereds);
803  }else{
804  pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius();
805  singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
806  }
807  } else if (_TargetBeam == 1) {
808  pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
809  singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
810  } else {
811  pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius();
812  singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
813  }
814  test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
815  }
816  }
817 
818  return pp;
819 }
820 
821 
822 //______________________________________________________________________________
823 void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
824  int&) // tcheck (unused)
825 {
826  // This function calculates momentum and energy of vector meson
827  // given W and Y, including interference.
828  // It gets the pt distribution from a lookup table.
829  double dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
830  int IY=0,j=0;
831 
832  dY = (_VMYmax-_VMYmin)/double(_VMnumy);
833 
834  // Y is already fixed; choose a pt
835  // Follow the approach in pickwy
836  // in _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
837  // Changed, now works -y to +y.
838  IY=int((Y-_VMYmin)/dY);
839  if (IY > (_VMnumy)-1){
840  IY=(_VMnumy)-1;
841  }
842 
843  yleft=(Y-_VMYmin)-(IY)*dY;
844 
845  yfract=yleft*dY;
846 
847  xpt=_randy.Rndom();
848  for(j=0;j<_VMNPT+1;j++){
849  if (xpt < _fptarray[IY][j]) goto L60;
850  }
851  L60:
852 
853  // now do linear interpolation - start with extremes
854  if (j == 0){
855  pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
856  goto L80;
857  }
858  if (j == _VMNPT){
859  pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
860  goto L80;
861  }
862 
863  // we're in the middle
864  ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
865  pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
866 
867  // at an extreme in y?
868  if (IY == (_VMnumy/2)-1){
869  pt=pt1;
870  goto L120;
871  }
872  L80:
873 
874  // interpolate in y repeat for next fractional y bin
875  for(j=0;j<_VMNPT+1;j++){
876  if (xpt < _fptarray[IY+1][j]) goto L90;
877  }
878  L90:
879 
880  // now do linear interpolation - start with extremes
881  if (j == 0){
882  pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
883  goto L100;
884  }
885  if (j == _VMNPT){
886  pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
887  goto L100;
888  }
889 
890  // we're in the middle
891  ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
892  pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
893  L100:
894 
895  // now interpolate in y
896  pt=yfract*pt2+(1-yfract)*pt1;
897  L120:
898 
899  // we have a pt
900  theta=2.*starlightConstants::pi*_randy.Rndom();
901  px=pt*cos(theta);
902  py=pt*sin(theta);
903 
904  E = sqrt(W*W+pt*pt)*cosh(Y);
905  pz = sqrt(W*W+pt*pt)*sinh(Y);
906  // randomly choose to make pz negative 50% of the time
907  if(_randy.Rndom()>=0.5) pz = -pz;
908 }
909 
910 
911 
912 //______________________________________________________________________________
913 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
914 {
915  double pT = sqrt(px*px + py*py);
916  double p = sqrt(pz*pz + pT*pT);
917  double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
918  return eta;
919 }
920 
921 //______________________________________________________________________________
922 Gammaanarrowvm::Gammaanarrowvm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
923 {
924  cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
925  read();
926  cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
927  narrowResonanceCrossSection sigma(input, bbsystem);
928  sigma.crossSectionCalculation(_bwnormsave);
929  setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
930  _VMbslope=sigma.slopeParameter();
931 }
932 
933 
934 //______________________________________________________________________________
935 Gammaanarrowvm::~Gammaanarrowvm()
936 { }
937 
938 
939 //______________________________________________________________________________
940 Gammaaincoherentvm::Gammaaincoherentvm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
941 {
942  cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
943  read();
944  cout<<"Creating and calculating crosssection. Gammaaincoherentvm()"<<endl;
945  incoherentVMCrossSection sigma(input, bbsystem);
946  sigma.crossSectionCalculation(_bwnormsave);
947  setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
948  _VMbslope=sigma.slopeParameter();
949 }
950 
951 
952 //______________________________________________________________________________
953 Gammaaincoherentvm::~Gammaaincoherentvm()
954 { }
955 
956 
957 //______________________________________________________________________________
958 Gammaawidevm::Gammaawidevm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
959 {
960  cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
961  read();
962  cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
963  wideResonanceCrossSection sigma(input, bbsystem);
964  sigma.crossSectionCalculation(_bwnormsave);
965  setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
966  _VMbslope=sigma.slopeParameter();
967 }
968 
969 
970 //______________________________________________________________________________
971 Gammaawidevm::~Gammaawidevm()
972 { }
973 
974 
975 //______________________________________________________________________________
976 e_Gammaanarrowvm::e_Gammaanarrowvm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
977 {
978  cout<<"Reading in luminosity tables. e_Gammaanarrowvm()"<<endl;
979  e_read();
980  cout<<"Creating and calculating crosssection. e_Gammaanarrowvm()"<<endl;
981  e_narrowResonanceCrossSection sigma(input, bbsystem);
982  sigma.crossSectionCalculation(_bwnormsave);
983  setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
984  _VMbslope=sigma.slopeParameter();
985 }
986 
987 
988 //______________________________________________________________________________
989 e_Gammaanarrowvm::~e_Gammaanarrowvm()
990 { }
991 
992 
993 //______________________________________________________________________________
994 e_Gammaawidevm::e_Gammaawidevm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
995 {
996  cout<<"Reading in luminosity tables. e_Gammaawidevm()"<<endl;
997  e_read();
998  cout<<"Creating and calculating crosssection. e_Gammaawidevm()"<<endl;
999  e_wideResonanceCrossSection sigma(input, bbsystem);
1000  sigma.crossSectionCalculation(_bwnormsave);
1001  setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
1002  _VMbslope=sigma.slopeParameter();
1003 }
1004 
1005 
1006 //______________________________________________________________________________
1007 e_Gammaawidevm::~e_Gammaawidevm()
1008 { }
1009 
1010 
1011 //______________________________________________________________________________
1012 void Gammaavectormeson::pickwEgamq2(double &W, double &cmsEgamma, double &targetEgamma,
1013  double &Q2, double &gamma_pz, double &gamma_pt,//photon in target frame
1014  double &E_prime, double &theta_e //electron
1015  )
1016 {
1017  double dW, dEgamma;
1018  double xw,xEgamma, xQ2, xtest, q2test; // btest;
1019  int IW,IGamma, IQ2;
1020  // ---------
1021  // int egamma_draws = 0, cms_egamma_draws =0, q2_draws =0 ;
1022  // ---------
1023  dW = (_VMWmax-_VMWmin)/double(_VMnumw);
1024  // Following used for debugging/timing for event generation
1025  //std::chrono::steady_clock::time_point begin_evt = std::chrono::steady_clock::now();
1026  bool pick_state = false;
1027  dEgamma = std::log(_targetMaxPhotonEnergy/_targetMinPhotonEnergy);
1028  while( pick_state == false ){
1029  xw = _randy.Rndom();
1030  W = _VMWmin + xw*(_VMWmax-_VMWmin);
1031  double w_test = _randy.Rndom();
1032  if (W < 2 * starlightConstants::pionChargedMass)
1033  continue;
1034  IW = int((W-_VMWmin)/dW);
1035  if (_BWarray[IW] < w_test)
1036  continue;
1037  xEgamma = _randy.Rndom();
1038  //
1039  targetEgamma = std::exp(std::log(_targetMinPhotonEnergy) + xEgamma*(dEgamma));
1040  IGamma = int(_VMnumega*xEgamma);
1041  // Holds Q2 and integrated Q2 dependence. Array is saved in target frame
1042  std::pair< double, std::vector<double> > this_energy = _g_EQ2array->operator[](IGamma);
1043  double intgrated_q2 = this_energy.first;
1044  // Sample single differential photon flux to obtain photon energy
1045  xtest = _randy.Rndom();
1046  if( xtest > intgrated_q2 ){
1047  //egamma_draws+=1;
1048  continue;
1049  }
1050  N0++;
1051  // btest = _randy.Rndom();
1052  // Load double differential photon flux table for selected energy
1053  std::vector<double> photon_flux = this_energy.second;
1054  double VMQ2min = photon_flux[0];
1055  double VMQ2max = photon_flux[1];
1056  //
1057  double ratio = std::log(VMQ2max/VMQ2min);
1058  double ln_min = std::log(VMQ2min);
1059 
1060  xQ2 = _randy.Rndom();
1061  Q2 = std::exp(ln_min+xQ2*ratio);
1062  IQ2 = int(_VMnumQ2*xQ2);
1063  // Load from look-up table. Use linear interpolation to evaluate at Q2
1064  double Q2_bin_0_1 = std::exp(ln_min+1*ratio/_VMnumQ2) - (std::exp(ln_min+0*ratio/_VMnumQ2));
1065  double x_1 = std::exp(ln_min+IQ2*ratio/_VMnumQ2);
1066  double x_2 = std::exp(ln_min+(1+IQ2)*ratio/_VMnumQ2);
1067  double x_3 = std::exp(ln_min+(2+IQ2)*ratio/_VMnumQ2);
1068  double scale_max = 0.0001;
1069  for(int ii = 0; ii < _VMnumQ2; ii++){
1070  double x_2_ii = std::exp(ln_min+(1+ii)*ratio/_VMnumQ2);
1071  double x_1_ii = std::exp(ln_min+ii*ratio/_VMnumQ2);
1072  if((photon_flux[ii+2]*(x_2_ii-x_1_ii)/Q2_bin_0_1 )>scale_max){
1073  scale_max = (photon_flux[ii+2]*(x_2_ii-x_1_ii)/Q2_bin_0_1);
1074  }
1075  }
1076  double y_1 = photon_flux[IQ2+2] *(x_2-x_1)/Q2_bin_0_1/scale_max;
1077  double y_2 = photon_flux[IQ2+3] *(x_3-x_2)/Q2_bin_0_1/scale_max;
1078  double m = (y_2 - y_1)/(x_2 - x_1);
1079  double c = y_1-m*x_1;
1080  double y = m*Q2+c;
1081  q2test = _randy.Rndom();
1082  if( y < q2test ){
1083  //q2_draws++;
1084  continue;
1085  }
1086  // -- Generate electron and photon in Target frame
1087  E_prime = _eEnergy - targetEgamma;
1088  if(Q2>1E-6){
1089  double cos_theta_e = 1. - Q2/(2.*_eEnergy*E_prime);
1090  theta_e = acos(cos_theta_e);
1091  }
1092  else {theta_e = sqrt(Q2/(_eEnergy*E_prime));}//updated using small angle approximation to avoid rounding to 0
1093  double beam_y = acosh(_targetBeamLorentzGamma)+_rap_CM;
1094  gamma_pt = E_prime*sin(theta_e);
1095 
1096  double pz_squared = targetEgamma*targetEgamma + Q2 - gamma_pt*gamma_pt;
1097  if( pz_squared < 0 )
1098  continue;
1099  double temp_pz = sqrt(pz_squared);
1100  // Now boost to CMS frame to check validity
1101  gamma_pz = temp_pz*cosh(beam_y) - targetEgamma*sinh(beam_y);
1102  cmsEgamma = targetEgamma*cosh(beam_y) - temp_pz*sinh(beam_y);
1103  // Simple checkl, should not be needed but used for safety
1104  if( cmsEgamma < _cmsMinPhotonEnergy || 2.*targetEgamma/(Q2+W*W) < _targetRadius){ //This cut is roughly RA = 0.8 fm the radius of proton and 1 eV^{-1} = 1.97 x 10 ^{-7} m
1105  continue;
1106  }
1107  xtest = _randy.Rndom();
1108  // Test against photonuclear cross section gamma+X --> VM+X
1109  if( _f_WYarray[IGamma][IQ2] < xtest ){
1110  continue;
1111  }
1112  pick_state = true;
1113  }
1114  return;
1115 }
1116 
1117 
1118 //______________________________________________________________________________
1119 eXEvent Gammaavectormeson::e_produceEvent()
1120 {
1121  // The new event type
1122  eXEvent event;
1123 
1124  int iFbadevent=0;
1125  int tcheck=0;
1126  starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
1127  starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
1128  // at present 4 prong decay is not implemented
1129  double comenergy = 0.;
1130  double rapidity = 0.;
1131  double Q2 = 0;
1132  double E = 0.;
1133  double momx=0.,momy=0.,momz=0.;
1134  double targetEgamma = 0, cmsEgamma = 0 ;
1135  double gamma_pz = 0 , gamma_pt = 0, e_theta = 0;
1136  double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
1137  double e_E=0., e_phi=0;
1138  double t_px =0, t_py=0., t_pz=0, t_E;
1139  bool accepted = false;
1140  do{
1141  pickwEgamq2(comenergy,cmsEgamma, targetEgamma,
1142  Q2, gamma_pz, gamma_pt, //photon infor in CMS frame
1143  e_E, e_theta); //electron info in target frame
1144  //
1145  momenta(comenergy,cmsEgamma, Q2, gamma_pz, gamma_pt, //input
1146  rapidity, E, momx, momy, momz, //VM
1147  t_px, t_py, t_pz, t_E, //pomeron
1148  e_phi,tcheck); //electron
1149  //
1150  // inelasticity: used for angular distributions
1151  double col_y = 1. - (e_E/_eEnergy)*std::pow(std::cos(e_theta/2.),2.);
1152  double col_polarization = (1 - col_y)/(1-col_y+col_y*col_y/2.);
1153  double spin_element = getSpinMatrixElement(comenergy, Q2, col_polarization);
1154  _nmbAttempts++;
1155 
1156  vmpid = ipid;
1157  // Two body dedcay in eSTARlight includes the angular corrections due to finite virtuality
1158  twoBodyDecay(ipid,comenergy,momx,momy,momz,spin_element,
1159  px1,py1,pz1,px2,py2,pz2,iFbadevent);
1160  double pt1chk = sqrt(px1*px1+py1*py1);
1161  double pt2chk = sqrt(px2*px2+py2*py2);
1162  double eta1 = pseudoRapidity(px1, py1, pz1);
1163  double eta2 = pseudoRapidity(px2, py2, pz2);
1164 
1165 
1166  if(_ptCutEnabled && !_etaCutEnabled){
1167  if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
1168  accepted = true;
1169  _nmbAccepted++;
1170  }
1171  }
1172  else if(!_ptCutEnabled && _etaCutEnabled){
1173  if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
1174  accepted = true;
1175  _nmbAccepted++;
1176  }
1177  }
1178  else if(_ptCutEnabled && _etaCutEnabled){
1179  if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
1180  if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
1181  accepted = true;
1182  _nmbAccepted++;
1183  }
1184  }
1185  }
1186  else if(!_ptCutEnabled && !_etaCutEnabled)
1187  _nmbAccepted++;
1188  }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
1189  if (iFbadevent==0&&tcheck==0) {
1190  int q1=0,q2=0;
1191  int ipid1,ipid2=0;
1192 
1193  double xtest = _randy.Rndom();
1194  if (xtest<0.5)
1195  {
1196  q1=1;
1197  q2=-1;
1198  }
1199  else {
1200  q1=-1;
1201  q2=1;
1202  }
1203 
1204  if ( ipid == 11 || ipid == 13 ){
1205  ipid1 = -q1*ipid;
1206  ipid2 = -q2*ipid;
1207  } else if (ipid == 22){ // omega --> pi0 + gamma
1208  ipid1 = 22;
1209  ipid2 = 111;
1210  } else {
1211  ipid1 = q1*ipid;
1212  ipid2 = q2*ipid;
1213  }
1214 
1215  // - Outgoing electron - target frame
1216  double e_ptot = sqrt(e_E*e_E - starlightConstants::mel*starlightConstants::mel);
1217  double e_px = e_ptot*sin(e_theta)*cos(e_phi);
1218  double e_py = e_ptot*sin(e_theta)*sin(e_phi);
1219  double e_pz = e_ptot*cos(e_theta);
1220  lorentzVector electron(e_px, e_py, e_pz, e_E);
1221  event.addSourceElectron(electron);
1222  // - Generated photon - target frame
1223  double gamma_x = gamma_pt*cos(e_phi+starlightConstants::pi);
1224  double gamma_y = gamma_pt*sin(e_phi+starlightConstants::pi);
1225  lorentzVector gamma(gamma_x,gamma_y,gamma_pz,cmsEgamma);
1226  vector3 boostVector(0, 0, tanh(_rap_CM));
1227  (gamma).Boost(boostVector);
1228  event.addGamma(gamma, targetEgamma, Q2);
1229  // - Saving V.M. daughters
1230  double md = getDaughterMass(vmpid);
1231  bool isOmegaNeutralDecay = (vmpid == starlightConstants::PHOTON);
1232  if(isOmegaNeutralDecay){
1233  //double mpi0 = starlightConstants::pionNeutralMass;
1234  double e_gamma1,px_gamma1,py_gamma1,pz_gamma1,e_gamma2,px_gamma2,py_gamma2,pz_gamma2;
1235  double Ed1 = sqrt(px1*px1+py1*py1+pz1*pz1);
1236  starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
1237  event.addParticle(particle1);
1238  //
1239  pi0Decay(px2,py2,pz2,e_gamma1,px_gamma1,py_gamma1,pz_gamma1,e_gamma2,px_gamma2,py_gamma2,pz_gamma2,iFbadevent);
1240  starlightParticle particle2(px_gamma1, py_gamma1, pz_gamma1, e_gamma1, starlightConstants::UNKNOWN, ipid1, q2);
1241  starlightParticle particle3(px_gamma2, py_gamma2, pz_gamma2, e_gamma2, starlightConstants::UNKNOWN, ipid1, q2);
1242  event.addParticle(particle2);
1243  event.addParticle(particle3);
1244  //double invmass = sqrt(mpi0*mpi0 + 2.0*(Ed1*Ed2 - (px1*px2+py1*py2+pz1*pz2)));
1245  } else {
1246  double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1);
1247  starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
1248  event.addParticle(particle1);
1249  //
1250  double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2);
1251  starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
1252  event.addParticle(particle2);
1253  }
1254 
1255  // - Scattered target and transfered momenta at target vertex
1256  double target_pz = - _beamNucleus*sqrt(_pEnergy*_pEnergy - pow(starlightConstants::protonMass,2.) ) - t_pz;
1257  //Sign of t_px in following equation changed to fix sign error and conserve p_z. Change made by Spencer Klein based on a bug report from Ya-Ping Xie. Nov. 14, 2019
1258  lorentzVector target(-t_px, -t_py, target_pz, _beamNucleus*_pEnergy - t_E);
1259  double t_var = t_E*t_E - t_px*t_px - t_py*t_py - t_pz*t_pz;
1260  event.addScatteredTarget(target, t_var);
1261  }
1262  return event;
1263 
1264 }
1265 string Gammaavectormeson::gammaTableParse(int ii, int jj)
1266 {
1267  ostringstream tag1, tag2;
1268  tag1<<ii;
1269  tag2<<jj;
1270  string to_ret = tag1.str()+","+tag2.str();
1271  return to_ret;
1272 }