Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-07-05 08:15:50 

0001 ///////////////////////////////////////////////////////////////////////////
0002 //
0003 //    Copyright 2010
0004 //
0005 //    This file is part of starlight.
0006 //
0007 //    starlight is free software: you can redistribute it and/or modify
0008 //    it under the terms of the GNU General Public License as published by
0009 //    the Free Software Foundation, either version 3 of the License, or
0010 //    (at your option) any later version.
0011 //
0012 //    starlight is distributed in the hope that it will be useful,
0013 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
0014 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
0015 //    GNU General Public License for more details.
0016 //
0017 //    You should have received a copy of the GNU General Public License
0018 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
0019 //
0020 ///////////////////////////////////////////////////////////////////////////
0021 //
0022 // File and Version Information:
0023 // $Rev:: 277                         $: revision of last commit
0024 // $Author:: jnystrand                $: author of last commit
0025 // $Date:: 2016-09-14 10:55:55 +0100 #$: date of last commit
0026 //
0027 // Description:
0028 //    Added incoherent t2-> pt2 selection.  Following pp selection scheme
0029 //
0030 //
0031 ///////////////////////////////////////////////////////////////////////////
0032 
0033 #include <iostream>
0034 #include <fstream>
0035 #include <cassert>
0036 #include <cmath>
0037 
0038 #include "gammaavm.h"
0039 // #include "photonNucleusCrossSection.h"
0040 #include "wideResonanceCrossSection.h"
0041 #include "narrowResonanceCrossSection.h"
0042 #include "incoherentVMCrossSection.h"
0043 //
0044 #include "e_narrowResonanceCrossSection.h"
0045 #include "e_wideResonanceCrossSection.h"
0046 
0047 using namespace std;
0048 
0049 //______________________________________________________________________________
0050 Gammaavectormeson::Gammaavectormeson(const inputParameters &inputParametersInstance, beamBeamSystem &bbsystem) : eventChannel(inputParametersInstance, bbsystem), _phaseSpaceGen(0)
0051 {
0052     _VMNPT = inputParametersInstance.nmbPtBinsInterference();
0053     _VMWmax = inputParametersInstance.maxW();
0054     _VMWmin = inputParametersInstance.minW();
0055     //_VMYmax=inputParametersInstance.maxRapidity();
0056     //_VMYmin=-1.*_VMYmax;
0057     _VMnumw = inputParametersInstance.nmbWBins();
0058     //_VMnumy=inputParametersInstance.nmbRapidityBins();
0059     _VMnumega = inputParametersInstance.nmbEnergyBins();
0060     _VMnumQ2 = inputParametersInstance.nmbGammaQ2Bins();
0061     _VMgamma_em = inputParametersInstance.beamLorentzGamma();
0062     _VMinterferencemode = inputParametersInstance.interferenceEnabled();
0063     _VMbslope = 0.; // Will define in wide/narrow constructor
0064     _bslopeDef = inputParametersInstance.bslopeDefinition();
0065     _bslopeVal = inputParametersInstance.bslopeValue();
0066     _pEnergy = inputParametersInstance.protonEnergy();
0067     _beamNucleus = inputParametersInstance.targetBeamA();
0068     // electron energy in CMS frame
0069     _eEnergy = inputParametersInstance.electronEnergy();
0070     _VMpidtest = inputParametersInstance.prodParticleType();
0071     _VMptmax = inputParametersInstance.maxPtInterference();
0072     _VMdpt = inputParametersInstance.ptBinWidthInterference();
0073     _ProductionMode = inputParametersInstance.productionMode();
0074     _targetMaxPhotonEnergy = inputParametersInstance.targetMaxPhotonEnergy();
0075     _targetMinPhotonEnergy = inputParametersInstance.targetMinPhotonEnergy();
0076     // Now saving the photon energy limits
0077     _cmsMaxPhotonEnergy = inputParametersInstance.cmsMaxPhotonEnergy();
0078     _cmsMinPhotonEnergy = inputParametersInstance.cmsMinPhotonEnergy();
0079     _beamLorentzGamma = inputParametersInstance.beamLorentzGamma();
0080     _targetBeamLorentzGamma = inputParametersInstance.targetBeamLorentzGamma();
0081     _rap_CM = inputParametersInstance.rap_CM();
0082     _targetRadius = inputParametersInstance.targetRadius();
0083     // Turn on/off backwards production
0084     _backwardsProduction = inputParametersInstance.backwardsProduction();
0085 
0086     N0 = 0;
0087     N1 = 0;
0088     N2 = 0;
0089     if (_VMpidtest == starlightConstants::FOURPRONG)
0090     {
0091         // create n-body phase-spage generator
0092         _phaseSpaceGen = new nBodyPhaseSpaceGen(_randy);
0093     }
0094     if (_ProductionMode == 12 || _ProductionMode == 13) // Narrow and wide coherent photon-pommeron in eSTARlight
0095         _dummy_pncs = new photonNucleusCrossSection(inputParametersInstance, bbsystem);
0096     //
0097     const double r_04_00 = 0.674;
0098     const double cos_delta = 0.925;
0099     for (int ii = 0; ii < 100; ++ii)
0100     { // epsilon 0-1
0101         double epsilon = 0.01 * ii;
0102         const double R = (1. / epsilon) * r_04_00 / (1. - r_04_00);
0103         for (int jj = 0; jj < 200; ++jj)
0104         { // psi 0 - 2pi
0105             double psi = jj * starlightConstants::pi / 100.;
0106             double max_bin;
0107             for (int kk = 0; kk < 200; ++kk)
0108             { // temp
0109                 double theta = kk * starlightConstants::pi / 100.;
0110                 // Fin max
0111                 double this_test = std::pow(std::sin(theta), 2.) * (1 + epsilon * cos(2. * psi)) + 2. * epsilon * R * std::pow(std::cos(theta), 2.) + std::sqrt(2. * epsilon * (1 + epsilon)) * cos_delta * std::sin(2. * theta) * std::cos(psi);
0112                 if (this_test > max_bin)
0113                     max_bin = this_test;
0114             }
0115             _angular_max[ii][jj] = max_bin;
0116         }
0117     }
0118 }
0119 
0120 //______________________________________________________________________________
0121 Gammaavectormeson::~Gammaavectormeson()
0122 {
0123     if (_phaseSpaceGen)
0124         delete _phaseSpaceGen;
0125     if (_dummy_pncs)
0126         delete _dummy_pncs;
0127 }
0128 
0129 //______________________________________________________________________________
0130 void Gammaavectormeson::pickwy(double &W, double &Y)
0131 {
0132     double dW, dY, xw, xy, xtest;
0133     int IW, IY;
0134 
0135     dW = (_VMWmax - _VMWmin) / double(_VMnumw);
0136     dY = (_VMYmax - _VMYmin) / double(_VMnumy);
0137 
0138 L201pwy:
0139 
0140     xw = _randy.Rndom();
0141     W = _VMWmin + xw * (_VMWmax - _VMWmin);
0142 
0143     if (W < 2 * starlightConstants::pionChargedMass)
0144         goto L201pwy;
0145 
0146     IW = int((W - _VMWmin) / dW);
0147     xy = _randy.Rndom();
0148     Y = _VMYmin + xy * (_VMYmax - _VMYmin);
0149     IY = int((Y - _VMYmin) / dY);
0150     xtest = _randy.Rndom();
0151 
0152     if (xtest > _Farray[IW][IY])
0153         goto L201pwy;
0154 
0155     N0++;
0156     // Determine the target nucleus
0157     // For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2
0158     _TargetBeam = 2; // Always true for eX
0159 }
0160 
0161 //______________________________________________________________________________
0162 void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
0163                                      double W,
0164                                      double px0, double py0, double pz0,
0165                                      double spin_element,
0166                                      double &px1, double &py1, double &pz1,
0167                                      double &px2, double &py2, double &pz2,
0168                                      int &iFbadevent)
0169 {
0170     // This routine decays a particle into two particles of mass mdec,
0171     // taking spin into account
0172     double pmag;
0173     double phi, theta, Ecm;
0174     double betax, betay, betaz;
0175     double mdec1 = 0.0, mdec2 = 0.0;
0176     double E1 = 0.0, E2 = 0.0;
0177 
0178     //    set the mass of the daughter particles
0179     mdec1 = getDaughterMass(ipid);
0180 
0181     if (_VMpidtest == starlightConstants::OMEGA_pi0gamma)
0182     {
0183         mdec2 = starlightConstants::pionNeutralMass;
0184         if (W < mdec2)
0185         {
0186             cout << " ERROR: W=" << W << endl;
0187             iFbadevent = 1;
0188             return;
0189         }
0190         //  calculate the magnitude of the momenta
0191         pmag = (W * W - mdec2 * mdec2) / (2.0 * W);
0192     }
0193     else
0194     {
0195         mdec2 = mdec1;
0196         if (W < 2 * mdec1)
0197         {
0198             cout << " ERROR: W=" << W << endl;
0199             iFbadevent = 1;
0200             return;
0201         }
0202         //  calculate the magnitude of the momenta
0203         pmag = sqrt(W * W / 4. - mdec2 * mdec2);
0204     }
0205 
0206     //  pick an orientation, based on the spin
0207     //  phi has a flat distribution in 2*pi
0208     phi = _randy.Rndom() * 2. * starlightConstants::pi;
0209 
0210     //  find theta, the angle between one of the outgoing particles and
0211     //  the beamline, in the frame of the two photons
0212     // cout<<"spin element: "<<spin_element<<endl;
0213     theta = getTheta(ipid, spin_element);
0214 
0215     //  compute unboosted momenta
0216     px1 = sin(theta) * cos(phi) * pmag;
0217     py1 = sin(theta) * sin(phi) * pmag;
0218     pz1 = cos(theta) * pmag;
0219     px2 = -px1;
0220     py2 = -py1;
0221     pz2 = -pz1;
0222 
0223     Ecm = sqrt(W * W + px0 * px0 + py0 * py0 + pz0 * pz0);
0224     E1 = sqrt(mdec1 * mdec1 + px1 * px1 + py1 * py1 + pz1 * pz1);
0225     E2 = sqrt(mdec2 * mdec2 + px2 * px2 + py2 * py2 + pz2 * pz2);
0226 
0227     // cout<<"Ecm: "<<Ecm<<endl;
0228     // cout<<"W: "<<W<<endl;
0229 
0230     betax = -(px0 / Ecm);
0231     betay = -(py0 / Ecm);
0232     betaz = -(pz0 / Ecm);
0233 
0234     transform(betax, betay, betaz, E1, px1, py1, pz1, iFbadevent);
0235     transform(betax, betay, betaz, E2, px2, py2, pz2, iFbadevent);
0236 
0237     if (iFbadevent == 1)
0238         return;
0239 }
0240 
0241 //______________________________________________________________________________
0242 // decays a particle into four particles with isotropic angular distribution
0243 bool Gammaavectormeson::fourBodyDecay(starlightConstants::particleTypeEnum &ipid,
0244                                       const double,             // E (unused)
0245                                       const double W,           // mass of produced particle
0246                                       const double *p,          // momentum of produced particle; expected to have size 3
0247                                       lorentzVector *decayVecs, // array of Lorentz vectors of daughter particles; expected to have size 4
0248                                       int &iFbadevent)
0249 {
0250     const double parentMass = W;
0251 
0252     // set the mass of the daughter particles
0253     const double daughterMass = getDaughterMass(ipid);
0254     if (parentMass < 4 * daughterMass)
0255     {
0256         cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
0257         iFbadevent = 1;
0258         return false;
0259     }
0260 
0261     // construct parent four-vector
0262     const double parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2] + parentMass * parentMass);
0263     const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
0264 
0265     // setup n-body phase-space generator
0266     assert(_phaseSpaceGen);
0267     static bool firstCall = true;
0268     if (firstCall)
0269     {
0270         const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
0271         _phaseSpaceGen->setDecay(4, m);
0272         // estimate maximum phase-space weight
0273         _phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
0274         firstCall = false;
0275     }
0276 
0277     // generate phase-space event
0278     if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
0279         return false;
0280 
0281     // set Lorentzvectors of decay daughters
0282     for (unsigned int i = 0; i < 4; ++i)
0283         decayVecs[i] = _phaseSpaceGen->daughter(i);
0284     return true;
0285 }
0286 
0287 //______________________________________________________________________________
0288 void Gammaavectormeson::pi0Decay(double &px_pi0, double &py_pi0, double &pz_pi0,
0289                                  double &e_g1, double &px_g1, double &py_g1, double &pz_g1,
0290                                  double &e_g2, double &px_g2, double &py_g2, double &pz_g2,
0291                                  int &iFbadevent)
0292 {
0293     double pmag;
0294     double phi, theta, Ecm;
0295     double betax, betay, betaz;
0296     double E1 = 0.0, E2 = 0.0;
0297 
0298     // This routine decays a pi0 into two isotropically produced photons
0299     double m_pi0 = starlightConstants::pionNeutralMass;
0300 
0301     //  calculate the magnitude of the momenta
0302     pmag = sqrt(m_pi0 * m_pi0 / 4.);
0303 
0304     //  pick an orientation, based on the spin
0305     //  phi has a flat distribution in 2*pi
0306     phi = _randy.Rndom() * 2. * starlightConstants::pi;
0307 
0308     //  find theta, the angle between one of the outgoing photons and
0309     //  the beamline, in the frame of the pi0
0310     theta = getTheta(starlightConstants::PION, 1.0 / 3.0);
0311 
0312     //  compute unboosted momenta
0313     px_g1 = sin(theta) * cos(phi) * pmag;
0314     py_g1 = sin(theta) * sin(phi) * pmag;
0315     pz_g1 = cos(theta) * pmag;
0316     px_g2 = -px_g1;
0317     py_g2 = -py_g1;
0318     pz_g2 = -pz_g1;
0319 
0320     Ecm = sqrt(m_pi0 * m_pi0 + px_pi0 * px_pi0 + py_pi0 * py_pi0 + pz_pi0 * pz_pi0);
0321     E1 = sqrt(px_g1 * px_g1 + py_g1 * py_g1 + pz_g1 * pz_g1);
0322     E2 = sqrt(px_g2 * px_g2 + py_g2 * py_g2 + pz_g2 * pz_g2);
0323 
0324     betax = -(px_pi0 / Ecm);
0325     betay = -(py_pi0 / Ecm);
0326     betaz = -(pz_pi0 / Ecm);
0327 
0328     transform(betax, betay, betaz, E1, px_g1, py_g1, pz_g1, iFbadevent);
0329     transform(betax, betay, betaz, E2, px_g2, py_g2, pz_g2, iFbadevent);
0330 
0331     e_g1 = E1;
0332     e_g2 = E2;
0333 
0334     if (iFbadevent == 1)
0335         return;
0336 }
0337 
0338 //______________________________________________________________________________
0339 double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
0340 {
0341     // This will return the daughter particles mass, and the final particles outputed id...
0342     double mdec = 0.;
0343 
0344     switch (_VMpidtest)
0345     {
0346     case starlightConstants::RHO:
0347     case starlightConstants::RHOZEUS:
0348     case starlightConstants::FOURPRONG:
0349     case starlightConstants::OMEGA:
0350         mdec = starlightConstants::pionChargedMass;
0351         ipid = starlightConstants::PION;
0352         break;
0353     case starlightConstants::OMEGA_pi0gamma:
0354         mdec = 0.0;
0355         ipid = starlightConstants::PHOTON;
0356         break;
0357     case starlightConstants::PHI:
0358         mdec = starlightConstants::kaonChargedMass;
0359         ipid = starlightConstants::KAONCHARGE;
0360         break;
0361     case starlightConstants::JPSI:
0362         mdec = starlightConstants::mel;
0363         ipid = starlightConstants::ELECTRON;
0364         break;
0365     case starlightConstants::JPSI_ee:
0366         mdec = starlightConstants::mel;
0367         ipid = starlightConstants::ELECTRON;
0368         break;
0369     case starlightConstants::JPSI_mumu:
0370         mdec = starlightConstants::muonMass;
0371         ipid = starlightConstants::MUON;
0372         break;
0373     case starlightConstants::JPSI2S_ee:
0374         mdec = starlightConstants::mel;
0375         ipid = starlightConstants::ELECTRON;
0376         break;
0377     case starlightConstants::JPSI2S_mumu:
0378         mdec = starlightConstants::muonMass;
0379         ipid = starlightConstants::MUON;
0380         break;
0381 
0382     case starlightConstants::JPSI2S:
0383     case starlightConstants::UPSILON:
0384     case starlightConstants::UPSILON2S:
0385     case starlightConstants::UPSILON3S:
0386         mdec = starlightConstants::muonMass;
0387         ipid = starlightConstants::MUON;
0388         break;
0389     case starlightConstants::UPSILON_ee:
0390     case starlightConstants::UPSILON2S_ee:
0391     case starlightConstants::UPSILON3S_ee:
0392         mdec = starlightConstants::mel;
0393         ipid = starlightConstants::ELECTRON;
0394         break;
0395     case starlightConstants::UPSILON_mumu:
0396     case starlightConstants::UPSILON2S_mumu:
0397     case starlightConstants::UPSILON3S_mumu:
0398         mdec = starlightConstants::muonMass;
0399         ipid = starlightConstants::MUON;
0400         break;
0401     default:
0402         cout << "No daughtermass defined, gammaavectormeson::getdaughtermass" << endl;
0403     }
0404 
0405     return mdec;
0406 }
0407 
0408 //______________________________________________________________________________
0409 double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid, double r_04_00)
0410 {
0411     // This depends on the decay angular distribution
0412     // Valid for rho, phi, omega.
0413     double theta = 0.;
0414     double xtest = 1.;
0415     double dndtheta = 0.;
0416 
0417     while (xtest > dndtheta)
0418     {
0419 
0420         theta = starlightConstants::pi * _randy.Rndom();
0421         xtest = _randy.Rndom();
0422         //  Follow distribution for helicity +/-1 with finite Q2
0423         //  Physics Letters B 449, 328; The European Physical Journal C - Particles and Fields 13, 37;
0424         //  The European Physical Journal C - Particles and Fields 6, 603
0425 
0426         switch (ipid)
0427         {
0428 
0429         case starlightConstants::MUON:
0430         case starlightConstants::ELECTRON:
0431             // primarily for upsilon/j/psi.  VM->ee/mumu
0432             dndtheta = sin(theta) * (1 + r_04_00 + (1 - 3. * r_04_00) * cos(theta) * cos(theta));
0433             break;
0434 
0435         case starlightConstants::PION:
0436         case starlightConstants::KAONCHARGE:
0437             // rhos etc
0438             dndtheta = sin(theta) * (1 - r_04_00 + (3. * r_04_00 - 1) * cos(theta) * cos(theta));
0439             break;
0440 
0441         default:
0442             if (!_backwardsProduction)
0443                 cout << "No proper theta dependence defined, check gammaavectormeson::gettheta" << endl;
0444         } // end of switch
0445 
0446         // Assume unpolarized vector-mesons for now
0447         if (_backwardsProduction)
0448         {
0449             dndtheta = sin(theta);
0450         }
0451     }
0452     return theta;
0453 }
0454 
0455 //______________________________________________________________________________
0456 double Gammaavectormeson::getSpinMatrixElement(double W, double Q2, double epsilon)
0457 {
0458     double m2 = 0.6 * (W * W);
0459     double R = starlightConstants::pi / 2. * m2 / Q2 - std::pow(m2, 3. / 2.) / sqrt(Q2) / (Q2 + m2) - m2 / Q2 * atan(sqrt(m2 / Q2));
0460     //////////needed for problem with double precision
0461     double R1 = starlightConstants::pi / 2. * m2 / Q2;
0462     if (R1 > 1.0E10)
0463         R = 0.0;
0464     //////////
0465     R = (Q2 + m2) * (Q2 + m2) * R * R / m2;
0466     return epsilon * R / (1 + epsilon * R);
0467 }
0468 
0469 //______________________________________________________________________________
0470 double Gammaavectormeson::getWidth()
0471 {
0472     return _width;
0473 }
0474 
0475 //______________________________________________________________________________
0476 double Gammaavectormeson::getMass()
0477 {
0478     return _mass;
0479 }
0480 
0481 //______________________________________________________________________________
0482 double Gammaavectormeson::getSpin()
0483 {
0484     return 1.0; // VM spins are the same
0485 }
0486 
0487 //______________________________________________________________________________
0488 void Gammaavectormeson::momenta(double W, double Egam, double Q2, double gamma_pz, double gamma_pt, // input conditions
0489                                 double &Y, double &E, double &px, double &py, double &pz,           // return vm
0490                                 double &t_px, double &t_py, double &t_pz, double &t_E,              // return pomeron
0491                                 double &e_phi, int &tcheck)                                         // return electron (angle already known by Q2)
0492 {
0493     //     This subroutine calculates momentum and energy of vector meson for electroproduction (eSTARlight)
0494     //     given W and photon 4-vector,   without interference.  No intereference in asymetric eX collisions
0495 
0496     double Epom, Pom_pz, tmin, pt2, phi1, phi2;
0497     double px1, py1;
0498     double xt, xtest, ytest;
0499     double t2;
0500 
0501     double target_px, target_py, target_pz, target_E;
0502 
0503     target_E = _beamNucleus * _pEnergy;
0504     target_px = 0.0;
0505     target_py = 0.0;
0506     target_pz = -_beamNucleus * sqrt(_pEnergy * _pEnergy - pow(starlightConstants::protonMass, 2.));
0507     phi1 = 2. * starlightConstants::pi * _randy.Rndom();
0508     px1 = gamma_pt * cos(phi1);
0509     py1 = gamma_pt * sin(phi1);
0510     int isbadevent = 0;
0511     double betax_cm = ((px1 + target_px) / (Egam + target_E));
0512     double betay_cm = ((py1 + target_py) / (Egam + target_E));
0513     double betaz_cm = ((gamma_pz + target_pz) / (Egam + target_E));
0514     transform(betax_cm, betay_cm, betaz_cm, target_E, target_px, target_py, target_pz, isbadevent);
0515     transform(betax_cm, betay_cm, betaz_cm, Egam, px1, py1, gamma_pz, isbadevent);
0516 
0517     e_phi = starlightConstants::pi + phi1;
0518     // Pomeron pz is != than its energy in eSTARlight, in order to conserve energy/momentum of scattered
0519     // target
0520     // Pom_pz = 0.5*(W*W-Q2)/(Egam + gamma_pz);
0521     Epom = 0.5 * (W * W + Q2) / (Egam + target_E);
0522 
0523 L522vm:
0524     while (e_phi > 2. * starlightConstants::pi)
0525         e_phi -= 2. * starlightConstants::pi;
0526     //
0527     if ((_bbs.targetBeam().A() == 1) || (_ProductionMode == 4))
0528     {
0529         if ((_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA))
0530         {
0531             // Use dipole form factor for light VM
0532         L613vm:
0533             xtest = 2.0 * _randy.Rndom();
0534             double ttest = xtest * xtest;
0535             ytest = _randy.Rndom();
0536             double t0 = 1. / 2.23;
0537             double yprob = xtest * _bbs.electronBeam().dipoleFormFactor(ttest, t0) * _bbs.electronBeam().dipoleFormFactor(ttest, t0);
0538             if (ytest > yprob)
0539                 goto L613vm;
0540             t2 = ttest;
0541             pt2 = xtest;
0542         }
0543         else
0544         {
0545             // Use dsig/dt= exp(-_VMbslope*t) for heavy VM
0546             double bslope_tdist = _VMbslope;
0547             double Wgammap = 0.0;
0548             switch (_bslopeDef)
0549             {
0550             case 0:
0551                 // This is the default, as before
0552                 bslope_tdist = _VMbslope;
0553                 break;
0554             case 1:
0555                 // User defined value of bslope. BSLOPE_VALUE default is 4.0 if not set.
0556                 bslope_tdist = _bslopeVal;
0557                 if (N0 <= 1)
0558                     cout << " ATTENTION: Using user defined value of bslope = " << _bslopeVal << endl;
0559                 break;
0560             case 2:
0561                 // This is Wgammap dependence of b from H1 (Eur. Phys. J. C 46 (2006) 585)
0562                 Wgammap = sqrt(4. * Egam * _pEnergy);
0563                 bslope_tdist = 4.63 + 4. * 0.164 * log(Wgammap / 90.0);
0564                 if (N0 <= 1)
0565                     cout << " ATTENTION: Using energy dependent value of bslope!" << endl;
0566                 break;
0567             default:
0568                 cout << " Undefined setting for BSLOPE_DEFINITION " << endl;
0569             }
0570 
0571             xtest = _randy.Rndom();
0572             // t2 = (-1./_VMbslope)*log(xtest);
0573             t2 = (-1. / bslope_tdist) * log(xtest);
0574             pt2 = sqrt(1. * t2);
0575         }
0576     }
0577     else
0578     {
0579         // >> Check tmin
0580         tmin = ((Epom / _VMgamma_em) * (Epom / _VMgamma_em));
0581 
0582         if (tmin > 0.5)
0583         {
0584             cout << " WARNING: tmin= " << tmin << endl;
0585             cout << " Y = " << Y << " W = " << W << " Epom = " << Epom << " gamma = " << _VMgamma_em << endl;
0586             cout << " Will pick a new W,Y " << endl;
0587             tcheck = 1;
0588             return;
0589         }
0590     L663vm:
0591         xt = _randy.Rndom();
0592         if (_bbs.targetBeam().A() != 1)
0593         {
0594             pt2 = 8. * xt * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0595         }
0596         else
0597         {
0598             std::cout << "Can't find the electron for eX" << std::endl;
0599         }
0600 
0601         xtest = _randy.Rndom();
0602         t2 = tmin + pt2 * pt2;
0603 
0604         double comp = 0.0;
0605         if (_bbs.targetBeam().A() != 1)
0606         {
0607             comp = _bbs.targetBeam().formFactor(t2) * _bbs.targetBeam().formFactor(t2) * pt2;
0608         }
0609         else
0610             std::cout << "Can't find the electron for eX" << std::endl;
0611         if (xtest > comp)
0612             goto L663vm;
0613     }
0614     phi2 = 2. * starlightConstants::pi * _randy.Rndom();
0615 
0616     //
0617     t_px = pt2 * cos(phi2);
0618     t_py = pt2 * sin(phi2);
0619     // Used to return the pomeron pz to generator
0620     // Compute scattered target kinematics p_(initial ion) = p_(pomeron) + p_(final ion)
0621     double newion_E = target_E - Epom;
0622     double newion_px = target_px - t_px;
0623     double newion_py = target_py - t_py;
0624     double newion_pz_squared = newion_E * newion_E - newion_px * newion_px - newion_py * newion_py - pow(_beamNucleus * starlightConstants::protonMass, 2.);
0625     if (newion_pz_squared < 0)
0626         goto L522vm;
0627     double newion_pz = -sqrt(newion_pz_squared);
0628     Pom_pz = target_pz - newion_pz;
0629     t_pz = Pom_pz;
0630     // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
0631     px = px1 + t_px;
0632     py = py1 + t_py;
0633 
0634     t_E = Epom;
0635     // Finally V.M. energy, pz and rapidity from photon + pommeron.
0636     E = Egam + t_E;
0637     pz = gamma_pz + t_pz;
0638     // transform back to e-ion CM frame
0639     transform(-betax_cm, -betay_cm, -betaz_cm, target_E, target_px, target_py, target_pz, isbadevent);
0640     transform(-betax_cm, -betay_cm, -betaz_cm, newion_E, newion_px, newion_py, newion_pz, isbadevent);
0641     transform(-betax_cm, -betay_cm, -betaz_cm, Egam, px1, py1, gamma_pz, isbadevent);
0642     transform(-betax_cm, -betay_cm, -betaz_cm, E, px, py, pz, isbadevent);
0643     t_px = target_px - newion_px;
0644     t_py = target_py - newion_py;
0645     t_pz = target_pz - newion_pz;
0646     t_E = target_E - newion_E;
0647     Y = 0.5 * std::log((E + fabs(pz)) / (E - fabs(pz)));
0648     //  cout << "pz_final = " << pz << endl;
0649 
0650     // Backwards Production... updated by Zachary Sweger on 9/21/2021
0651     if (_backwardsProduction)
0652     {
0653         double is_proton_E = target_E;
0654         double is_proton_px = target_px;
0655         double is_proton_py = target_py;
0656         double is_proton_pz = target_pz;
0657         double is_gamma_E = Egam;
0658         double is_gamma_px = px1;
0659         double is_gamma_py = py1;
0660         double is_gamma_pz = gamma_pz;
0661         double is_tot_E = is_proton_E + is_gamma_E;
0662         double is_tot_px = is_proton_px + is_gamma_px;
0663         double is_tot_py = is_proton_py + is_gamma_py;
0664         double is_tot_pz = is_proton_pz + is_gamma_pz;
0665 
0666         double fs_proton_E = newion_E;
0667         double fs_proton_px = newion_px;
0668         double fs_proton_py = newion_py;
0669         double fs_proton_pz = newion_pz;
0670 
0671         double fs_vm_E = E;
0672         double fs_vm_px = px;
0673         double fs_vm_py = py;
0674         double fs_vm_pz = pz;
0675 
0676         int isbadevent = 0;
0677         double betax_cm = (is_tot_px / is_tot_E);
0678         double betay_cm = (is_tot_py / is_tot_E);
0679         double betaz_cm = (is_tot_pz / is_tot_E);
0680         transform(betax_cm, betay_cm, betaz_cm, is_proton_E, is_proton_px, is_proton_py, is_proton_pz, isbadevent);
0681         transform(betax_cm, betay_cm, betaz_cm, is_gamma_E, is_gamma_px, is_gamma_py, is_gamma_pz, isbadevent);
0682         transform(betax_cm, betay_cm, betaz_cm, fs_proton_E, fs_proton_px, fs_proton_py, fs_proton_pz, isbadevent);
0683         transform(betax_cm, betay_cm, betaz_cm, fs_vm_E, fs_vm_px, fs_vm_py, fs_vm_pz, isbadevent);
0684 
0685         double u_fs_proton_px = -1.0 * fs_proton_px;
0686         double u_fs_proton_py = -1.0 * fs_proton_py;
0687         double u_fs_proton_pz = -1.0 * fs_proton_pz;
0688         double u_fs_proton_E = fs_proton_E;
0689 
0690         double u_fs_vm_px = -1.0 * fs_vm_px;
0691         double u_fs_vm_py = -1.0 * fs_vm_py;
0692         double u_fs_vm_pz = -1.0 * fs_vm_pz;
0693         double u_fs_vm_E = fs_vm_E;
0694 
0695         betax_cm = -1.0 * betax_cm;
0696         betay_cm = -1.0 * betay_cm;
0697         betaz_cm = -1.0 * betaz_cm;
0698         transform(betax_cm, betay_cm, betaz_cm, is_proton_E, is_proton_px, is_proton_py, is_proton_pz, isbadevent);
0699         transform(betax_cm, betay_cm, betaz_cm, is_gamma_E, is_gamma_px, is_gamma_py, is_gamma_pz, isbadevent);
0700         transform(betax_cm, betay_cm, betaz_cm, u_fs_proton_E, u_fs_proton_px, u_fs_proton_py, u_fs_proton_pz, isbadevent);
0701         transform(betax_cm, betay_cm, betaz_cm, u_fs_vm_E, u_fs_vm_px, u_fs_vm_py, u_fs_vm_pz, isbadevent);
0702 
0703         px = u_fs_vm_px;
0704         py = u_fs_vm_py;
0705         pz = u_fs_vm_pz;
0706         E = u_fs_vm_E;
0707 
0708         t_px = -u_fs_proton_px;
0709         t_py = -u_fs_proton_py;
0710         t_pz = is_proton_pz - u_fs_proton_pz;
0711         t_E = is_proton_E - u_fs_proton_E;
0712         // cout<<"Energy Violation: "<<(E+u_fs_proton_E)-(is_tot_E)<<endl;
0713         // cout<<"Pz Violation: "<<(pz+u_fs_proton_pz)-(is_tot_pz)<<endl;
0714         // cout<<"Px Violation: "<<(px+u_fs_proton_px)-(is_tot_px)<<endl;
0715         // cout<<"Py Violation: "<<(py+u_fs_proton_py)-(is_tot_py)<<endl;
0716         // 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;
0717         // cout << "VM MASS = " << sqrt(E*E-px*px - py*py-pz*pz) << endl;
0718 
0719         // Calculate the rapidity
0720         Y = 0.5 * std::log((E + fabs(pz)) / (E - fabs(pz)));
0721         if (Y != Y)
0722         {
0723             // FIXME the following should be upated if not using omegas
0724             if (_VMpidtest == starlightConstants::OMEGA_pi0gamma || _VMpidtest == starlightConstants::OMEGA)
0725                 E = sqrt(starlightConstants::OmegaMass * starlightConstants::OmegaMass + px * px + py * py + pz * pz);
0726             else if (_VMpidtest == starlightConstants::RHO)
0727                 E = sqrt(starlightConstants::rho0Mass * starlightConstants::rho0Mass + px * px + py * py + pz * pz);
0728             // cout<<"Energy Violation: "<<(E+u_fs_proton_E)-(is_tot_E)<<endl;
0729         }
0730     }
0731 }
0732 
0733 //______________________________________________________________________________
0734 double Gammaavectormeson::pTgamma(double E)
0735 {
0736     // returns on random draw from pp(E) distribution
0737     double ereds = 0., Cm = 0., Coef = 0., x = 0., pp = 0., test = 0., u = 0.;
0738     double singleformfactorCm = 0., singleformfactorpp1 = 0., singleformfactorpp2 = 0.;
0739     int satisfy = 0;
0740 
0741     ereds = (E / _VMgamma_em) * (E / _VMgamma_em);
0742     // sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
0743     Cm = sqrt(3.) * E / _VMgamma_em;
0744     // If E is very small, the drawing of a pT below is extremely slow.
0745     // ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E.
0746     // Should have no observable consequences (JN, SRK 11-Sep-2014)
0747     if (E < 0.0005)
0748         return Cm;
0749 
0750     // the amplitude of the p_t spectrum at the maximum
0751 
0752     if (_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() != 1)
0753     {
0754         if (_ProductionMode == 2 || _ProductionMode == 3)
0755         {
0756             singleformfactorCm = _bbs.electronBeam().formFactor(Cm * Cm + ereds);
0757         }
0758         else
0759         {
0760             singleformfactorCm = _bbs.targetBeam().formFactor(Cm * Cm + ereds);
0761         }
0762     }
0763     else if (_bbs.targetBeam().A() == 1 && _bbs.electronBeam().A() != 1)
0764     {
0765         if (_ProductionMode == 2 || _ProductionMode == 3)
0766         {
0767             singleformfactorCm = _bbs.targetBeam().formFactor(Cm * Cm + ereds);
0768         }
0769         else
0770         {
0771             singleformfactorCm = _bbs.electronBeam().formFactor(Cm * Cm + ereds);
0772         }
0773     }
0774     else if (_TargetBeam == 1)
0775     {
0776         singleformfactorCm = _bbs.targetBeam().formFactor(Cm * Cm + ereds);
0777     }
0778     else
0779     {
0780         singleformfactorCm = _bbs.electronBeam().formFactor(Cm * Cm + ereds);
0781     }
0782 
0783     Coef = 3.0 * (singleformfactorCm * singleformfactorCm * Cm * Cm * Cm) / ((2. * (starlightConstants::pi) * (ereds + Cm * Cm)) * (2. * (starlightConstants::pi) * (ereds + Cm * Cm)));
0784 
0785     // pick a test value pp, and find the amplitude there
0786     x = _randy.Rndom();
0787 
0788     if (_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() != 1)
0789     {
0790         if (_ProductionMode == 2 || _ProductionMode == 3)
0791         {
0792             pp = x * 5. * starlightConstants::hbarc / _bbs.electronBeam().nuclearRadius();
0793             singleformfactorpp1 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0794         }
0795         else
0796         {
0797             pp = x * 5. * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0798             singleformfactorpp1 = _bbs.targetBeam().formFactor(pp * pp + ereds);
0799         }
0800     }
0801     else if (_bbs.targetBeam().A() == 1 && _bbs.electronBeam().A() != 1)
0802     {
0803         if (_ProductionMode == 2 || _ProductionMode == 3)
0804         {
0805             pp = x * 5. * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0806             singleformfactorpp1 = _bbs.targetBeam().formFactor(pp * pp + ereds);
0807         }
0808         else
0809         {
0810             pp = x * 5. * starlightConstants::hbarc / _bbs.electronBeam().nuclearRadius();
0811             singleformfactorpp1 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0812         }
0813     }
0814     else if (_TargetBeam == 1)
0815     {
0816         pp = x * 5. * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0817         singleformfactorpp1 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0818     }
0819     else
0820     {
0821         pp = x * 5. * starlightConstants::hbarc / _bbs.electronBeam().nuclearRadius();
0822         singleformfactorpp1 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0823     }
0824 
0825     test = (singleformfactorpp1 * singleformfactorpp1) * pp * pp * pp / ((2. * starlightConstants::pi * (ereds + pp * pp)) * (2. * starlightConstants::pi * (ereds + pp * pp)));
0826 
0827     while (satisfy == 0)
0828     {
0829         u = _randy.Rndom();
0830         if (u * Coef <= test)
0831         {
0832             satisfy = 1;
0833         }
0834         else
0835         {
0836             x = _randy.Rndom();
0837             if (_bbs.electronBeam().A() == 1 && _bbs.targetBeam().A() != 1)
0838             {
0839                 if (_ProductionMode == 2 || _ProductionMode == 3)
0840                 {
0841                     pp = x * 5. * starlightConstants::hbarc / _bbs.electronBeam().nuclearRadius();
0842                     singleformfactorpp2 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0843                 }
0844                 else
0845                 {
0846                     pp = x * 5. * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0847                     singleformfactorpp2 = _bbs.targetBeam().formFactor(pp * pp + ereds);
0848                 }
0849             }
0850             else if (_bbs.targetBeam().A() == 1 && _bbs.electronBeam().A() != 1)
0851             {
0852                 if (_ProductionMode == 2 || _ProductionMode == 3)
0853                 {
0854                     pp = x * 5. * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0855                     singleformfactorpp2 = _bbs.targetBeam().formFactor(pp * pp + ereds);
0856                 }
0857                 else
0858                 {
0859                     pp = x * 5. * starlightConstants::hbarc / _bbs.electronBeam().nuclearRadius();
0860                     singleformfactorpp2 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0861                 }
0862             }
0863             else if (_TargetBeam == 1)
0864             {
0865                 pp = x * 5. * starlightConstants::hbarc / _bbs.targetBeam().nuclearRadius();
0866                 singleformfactorpp2 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0867             }
0868             else
0869             {
0870                 pp = x * 5. * starlightConstants::hbarc / _bbs.electronBeam().nuclearRadius();
0871                 singleformfactorpp2 = _bbs.electronBeam().formFactor(pp * pp + ereds);
0872             }
0873             test = (singleformfactorpp2 * singleformfactorpp2) * pp * pp * pp / (2. * starlightConstants::pi * (ereds + pp * pp) * 2. * starlightConstants::pi * (ereds + pp * pp));
0874         }
0875     }
0876 
0877     return pp;
0878 }
0879 
0880 //______________________________________________________________________________
0881 void Gammaavectormeson::vmpt(double W, double Y, double &E, double &px, double &py, double &pz,
0882                              int &) // tcheck (unused)
0883 {
0884     //    This function calculates momentum and energy of vector meson
0885     //    given W and Y, including interference.
0886     //    It gets the pt distribution from a lookup table.
0887     double dY = 0., yleft = 0., yfract = 0., xpt = 0., pt1 = 0., ptfract = 0., pt = 0., pt2 = 0., theta = 0.;
0888     int IY = 0, j = 0;
0889 
0890     dY = (_VMYmax - _VMYmin) / double(_VMnumy);
0891 
0892     //  Y is already fixed; choose a pt
0893     //  Follow the approach in pickwy
0894     //  in  _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
0895     // Changed,  now works -y to +y.
0896     IY = int((Y - _VMYmin) / dY);
0897     if (IY > (_VMnumy)-1)
0898     {
0899         IY = (_VMnumy)-1;
0900     }
0901 
0902     yleft = (Y - _VMYmin) - (IY)*dY;
0903 
0904     yfract = yleft * dY;
0905 
0906     xpt = _randy.Rndom();
0907     for (j = 0; j < _VMNPT + 1; j++)
0908     {
0909         if (xpt < _fptarray[IY][j])
0910             goto L60;
0911     }
0912 L60:
0913 
0914     //  now do linear interpolation - start with extremes
0915     if (j == 0)
0916     {
0917         pt1 = xpt / _fptarray[IY][j] * _VMdpt / 2.;
0918         goto L80;
0919     }
0920     if (j == _VMNPT)
0921     {
0922         pt1 = (_VMptmax - _VMdpt / 2.) + _VMdpt / 2. * (xpt - _fptarray[IY][j]) / (1. - _fptarray[IY][j]);
0923         goto L80;
0924     }
0925 
0926     //  we're in the middle
0927     ptfract = (xpt - _fptarray[IY][j]) / (_fptarray[IY][j + 1] - _fptarray[IY][j]);
0928     pt1 = (j + 1) * _VMdpt + ptfract * _VMdpt;
0929 
0930     //  at an extreme in y?
0931     if (IY == (_VMnumy / 2) - 1)
0932     {
0933         pt = pt1;
0934         goto L120;
0935     }
0936 L80:
0937 
0938     //  interpolate in y repeat for next fractional y bin
0939     for (j = 0; j < _VMNPT + 1; j++)
0940     {
0941         if (xpt < _fptarray[IY + 1][j])
0942             goto L90;
0943     }
0944 L90:
0945 
0946     //  now do linear interpolation - start with extremes
0947     if (j == 0)
0948     {
0949         pt2 = xpt / _fptarray[IY + 1][j] * _VMdpt / 2.;
0950         goto L100;
0951     }
0952     if (j == _VMNPT)
0953     {
0954         pt2 = (_VMptmax - _VMdpt / 2.) + _VMdpt / 2. * (xpt - _fptarray[IY + 1][j]) / (1. - _fptarray[IY + 1][j]);
0955         goto L100;
0956     }
0957 
0958     //  we're in the middle
0959     ptfract = (xpt - _fptarray[IY + 1][j]) / (_fptarray[IY + 1][j + 1] - _fptarray[IY + 1][j]);
0960     pt2 = (j + 1) * _VMdpt + ptfract * _VMdpt;
0961 L100:
0962 
0963     //  now interpolate in y
0964     pt = yfract * pt2 + (1 - yfract) * pt1;
0965 L120:
0966 
0967     //  we have a pt
0968     theta = 2. * starlightConstants::pi * _randy.Rndom();
0969     px = pt * cos(theta);
0970     py = pt * sin(theta);
0971 
0972     E = sqrt(W * W + pt * pt) * cosh(Y);
0973     pz = sqrt(W * W + pt * pt) * sinh(Y);
0974     //      randomly choose to make pz negative 50% of the time
0975     if (_randy.Rndom() >= 0.5)
0976         pz = -pz;
0977 }
0978 
0979 //______________________________________________________________________________
0980 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
0981 {
0982     double pT = sqrt(px * px + py * py);
0983     double p = sqrt(pz * pz + pT * pT);
0984     double eta = -99.9;
0985     if ((p - pz) != 0)
0986     {
0987         eta = 0.5 * log((p + pz) / (p - pz));
0988     }
0989     return eta;
0990 }
0991 
0992 //______________________________________________________________________________
0993 Gammaanarrowvm::Gammaanarrowvm(const inputParameters &input, beamBeamSystem &bbsystem) : Gammaavectormeson(input, bbsystem)
0994 {
0995     cout << "Reading in luminosity tables. Gammaanarrowvm()" << endl;
0996     read();
0997     cout << "Creating and calculating crosssection. Gammaanarrowvm()" << endl;
0998     narrowResonanceCrossSection sigma(input, bbsystem);
0999     sigma.crossSectionCalculation(_bwnormsave);
1000     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
1001     _VMbslope = sigma.slopeParameter();
1002 }
1003 
1004 //______________________________________________________________________________
1005 Gammaanarrowvm::~Gammaanarrowvm()
1006 {
1007 }
1008 
1009 //______________________________________________________________________________
1010 Gammaaincoherentvm::Gammaaincoherentvm(const inputParameters &input, beamBeamSystem &bbsystem) : Gammaavectormeson(input, bbsystem)
1011 {
1012     cout << "Reading in luminosity tables. Gammaainkoherentvm()" << endl;
1013     read();
1014     cout << "Creating and calculating crosssection. Gammaaincoherentvm()" << endl;
1015     incoherentVMCrossSection sigma(input, bbsystem);
1016     sigma.crossSectionCalculation(_bwnormsave);
1017     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
1018     _VMbslope = sigma.slopeParameter();
1019 }
1020 
1021 //______________________________________________________________________________
1022 Gammaaincoherentvm::~Gammaaincoherentvm()
1023 {
1024 }
1025 
1026 //______________________________________________________________________________
1027 Gammaawidevm::Gammaawidevm(const inputParameters &input, beamBeamSystem &bbsystem) : Gammaavectormeson(input, bbsystem)
1028 {
1029     cout << "Reading in luminosity tables. Gammaawidevm()" << endl;
1030     read();
1031     cout << "Creating and calculating crosssection. Gammaawidevm()" << endl;
1032     wideResonanceCrossSection sigma(input, bbsystem);
1033     sigma.crossSectionCalculation(_bwnormsave);
1034     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
1035     _VMbslope = sigma.slopeParameter();
1036 }
1037 
1038 //______________________________________________________________________________
1039 Gammaawidevm::~Gammaawidevm()
1040 {
1041 }
1042 
1043 //______________________________________________________________________________
1044 e_Gammaanarrowvm::e_Gammaanarrowvm(const inputParameters &input, beamBeamSystem &bbsystem) : Gammaavectormeson(input, bbsystem)
1045 {
1046     cout << "Reading in luminosity tables. e_Gammaanarrowvm()" << endl;
1047     e_read();
1048     cout << "Creating and calculating crosssection. e_Gammaanarrowvm()" << endl;
1049     e_narrowResonanceCrossSection sigma(input, bbsystem);
1050     sigma.crossSectionCalculation(_bwnormsave);
1051     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
1052     _VMbslope = sigma.slopeParameter();
1053 }
1054 
1055 //______________________________________________________________________________
1056 e_Gammaanarrowvm::~e_Gammaanarrowvm()
1057 {
1058 }
1059 
1060 //______________________________________________________________________________
1061 e_Gammaawidevm::e_Gammaawidevm(const inputParameters &input, beamBeamSystem &bbsystem) : Gammaavectormeson(input, bbsystem)
1062 {
1063     cout << "Reading in luminosity tables. e_Gammaawidevm()" << endl;
1064     e_read();
1065     cout << "Creating and calculating crosssection. e_Gammaawidevm()" << endl;
1066     e_wideResonanceCrossSection sigma(input, bbsystem);
1067     sigma.crossSectionCalculation(_bwnormsave);
1068     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
1069     _VMbslope = sigma.slopeParameter();
1070 }
1071 
1072 //______________________________________________________________________________
1073 e_Gammaawidevm::~e_Gammaawidevm()
1074 {
1075 }
1076 
1077 //______________________________________________________________________________
1078 void Gammaavectormeson::pickwEgamq2(double &W, double &cmsEgamma, double &targetEgamma,
1079                                     double &Q2, double &gamma_pz, double &gamma_pt, // photon in target frame
1080                                     double &E_prime, double &theta_e                // electron
1081 )
1082 {
1083     double dW, dEgamma;
1084     double xw, xEgamma, xQ2, xtest, q2test; // btest;
1085     int IW, IGamma, IQ2;
1086     // ---------
1087     //  int egamma_draws = 0, cms_egamma_draws =0, q2_draws =0 ;
1088     // ---------
1089     dW = (_VMWmax - _VMWmin) / double(_VMnumw);
1090     // Following used for debugging/timing for event generation
1091     // std::chrono::steady_clock::time_point begin_evt = std::chrono::steady_clock::now();
1092     bool pick_state = false;
1093     dEgamma = std::log(_targetMaxPhotonEnergy / _targetMinPhotonEnergy);
1094 
1095     while (pick_state == false)
1096     {
1097         xw = _randy.Rndom();
1098         W = _VMWmin + xw * (_VMWmax - _VMWmin);
1099         double w_test = _randy.Rndom();
1100         if (W < 2 * starlightConstants::pionChargedMass)
1101             continue;
1102         IW = int((W - _VMWmin) / dW);
1103         if (_BWarray[IW] < w_test)
1104             continue;
1105         xEgamma = _randy.Rndom();
1106         //
1107         targetEgamma = std::exp(std::log(_targetMinPhotonEnergy) + xEgamma * (dEgamma));
1108         IGamma = int(_VMnumega * xEgamma);
1109         // Holds Q2 and integrated Q2 dependence. Array is saved in target frame
1110         std::pair<double, std::vector<double>> this_energy = _g_EQ2array->operator[](IGamma);
1111         double intgrated_q2 = this_energy.first;
1112         // Sample single differential photon flux to obtain photon energy
1113         xtest = _randy.Rndom();
1114         if (xtest > intgrated_q2)
1115         {
1116             // egamma_draws+=1;
1117             continue;
1118         }
1119         N0++;
1120         //    btest = _randy.Rndom();
1121         // Load double differential photon flux table for selected energy
1122         std::vector<double> photon_flux = this_energy.second;
1123         double VMQ2min = photon_flux[0];
1124         double VMQ2max = photon_flux[1];
1125         //
1126         double ratio = std::log(VMQ2max / VMQ2min);
1127         double ln_min = std::log(VMQ2min);
1128 
1129         xQ2 = _randy.Rndom();
1130         Q2 = std::exp(ln_min + xQ2 * ratio);
1131         IQ2 = int(_VMnumQ2 * xQ2);
1132         // Load from look-up table. Use linear interpolation to evaluate at Q2
1133         double Q2_bin_0_1 = std::exp(ln_min + 1 * ratio / _VMnumQ2) - (std::exp(ln_min + 0 * ratio / _VMnumQ2));
1134         double x_1 = std::exp(ln_min + IQ2 * ratio / _VMnumQ2);
1135         double x_2 = std::exp(ln_min + (1 + IQ2) * ratio / _VMnumQ2);
1136         double x_3 = std::exp(ln_min + (2 + IQ2) * ratio / _VMnumQ2);
1137         double scale_max = 0.0001;
1138         for (int ii = 0; ii < _VMnumQ2; ii++)
1139         {
1140             double x_2_ii = std::exp(ln_min + (1 + ii) * ratio / _VMnumQ2);
1141             double x_1_ii = std::exp(ln_min + ii * ratio / _VMnumQ2);
1142             if ((photon_flux[ii + 2] * (x_2_ii - x_1_ii) / Q2_bin_0_1) > scale_max)
1143             {
1144                 scale_max = (photon_flux[ii + 2] * (x_2_ii - x_1_ii) / Q2_bin_0_1);
1145             }
1146         }
1147         double y_1 = photon_flux[IQ2 + 2] * (x_2 - x_1) / Q2_bin_0_1 / scale_max;
1148         double y_2 = photon_flux[IQ2 + 3] * (x_3 - x_2) / Q2_bin_0_1 / scale_max;
1149         double m = (y_2 - y_1) / (x_2 - x_1);
1150         double c = y_1 - m * x_1;
1151         double y = m * Q2 + c;
1152         q2test = _randy.Rndom();
1153         if (y < q2test)
1154         {
1155             // q2_draws++;
1156             continue;
1157         }
1158         // -- Generate electron and photon in Target frame
1159         E_prime = _eEnergy - targetEgamma;
1160         if (Q2 > 1E-6)
1161         {
1162             double cos_theta_e = 1. - Q2 / (2. * _eEnergy * E_prime);
1163             theta_e = acos(cos_theta_e);
1164         }
1165         else
1166         {
1167             theta_e = sqrt(Q2 / (_eEnergy * E_prime));
1168         } // updated using small angle approximation to avoid rounding to 0
1169         double beam_y = acosh(_targetBeamLorentzGamma) + _rap_CM;
1170         gamma_pt = E_prime * sin(theta_e);
1171 
1172         double pz_squared = targetEgamma * targetEgamma + Q2 - gamma_pt * gamma_pt;
1173         if (pz_squared < 0)
1174             continue;
1175         double temp_pz = sqrt(pz_squared);
1176         // Now boost to CMS frame to check validity
1177         gamma_pz = temp_pz * cosh(beam_y) - targetEgamma * sinh(beam_y);
1178         cmsEgamma = targetEgamma * cosh(beam_y) - temp_pz * sinh(beam_y);
1179         // Simple checkl, should not be needed but used for safety
1180         if (cmsEgamma < _cmsMinPhotonEnergy || 2. * targetEgamma / (Q2 + W * W) < _targetRadius)
1181         { // This cut is roughly RA = 0.8 fm the radius of proton and 1 eV^{-1} = 1.97 x 10 ^{-7} m
1182             continue;
1183         }
1184         xtest = _randy.Rndom();
1185         // Test against photonuclear cross section gamma+X --> VM+X
1186         if (_f_WYarray[IGamma][IQ2] < xtest)
1187         {
1188             continue;
1189         }
1190         pick_state = true;
1191     }
1192     return;
1193 }
1194 
1195 //______________________________________________________________________________
1196 eXEvent Gammaavectormeson::e_produceEvent()
1197 {
1198     // The new event type
1199     eXEvent event;
1200 
1201     int iFbadevent = 0;
1202     int tcheck = 0;
1203     starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
1204     starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
1205     // at present 4 prong decay is not implemented
1206     double comenergy = 0.;
1207     double rapidity = 0.;
1208     double Q2 = 0;
1209     double E = 0.;
1210     double momx = 0., momy = 0., momz = 0.;
1211     double targetEgamma = 0, cmsEgamma = 0;
1212     double gamma_pz = 0, gamma_pt = 0, e_theta = 0;
1213     double px2 = 0., px1 = 0., py2 = 0., py1 = 0., pz2 = 0., pz1 = 0.;
1214     double e_E = 0., e_phi = 0;
1215     double t_px = 0, t_py = 0., t_pz = 0, t_E;
1216     bool accepted = false;
1217     if (_VMpidtest == starlightConstants::FOURPRONG)
1218     {
1219 
1220         double comenergy = 0;
1221         double mom[3] = {0, 0, 0};
1222         double E = 0;
1223         lorentzVector decayVecs[4];
1224         double rapidity = 0;
1225         do
1226         {
1227             tcheck = 0;     // reinitialized after every loop cycle - to avoid infinite loop
1228             iFbadevent = 0; // same as above.
1229             pickwEgamq2(comenergy, cmsEgamma, targetEgamma,
1230                         Q2, gamma_pz, gamma_pt,                   // photon infor in CMS frame
1231                         e_E, e_theta);                            // electron info in target frame
1232             momenta(comenergy, cmsEgamma, Q2, gamma_pz, gamma_pt, // input
1233                     rapidity, E, momx, momy, momz,                // VM
1234                     t_px, t_py, t_pz, t_E,                        // pomeron
1235                     e_phi, tcheck);                               // electron
1236             _nmbAttempts++;
1237             accepted = true; // re-initialized after every loop cycle -to avoid infinite loop
1238             mom[0] = momx;
1239             mom[1] = momy;
1240             mom[2] = momz;
1241 
1242             if (tcheck != 0 || !fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent))
1243             { // if either vector meson creation, or further decay into four pions, is impossible
1244                 accepted = false;
1245                 continue; // this skips the etaCut and ptCut checks.
1246             }
1247 
1248             if (_ptCutEnabled)
1249             {
1250                 for (int i = 0; i < 4; i++)
1251                 {
1252                     double pt_chk2 = 0;
1253                     pt_chk2 += pow(decayVecs[i].GetPx(), 2);
1254                     pt_chk2 += pow(decayVecs[i].GetPy(), 2); // compute transverse momentum (squared) for the particle, for ptCut checks.
1255 
1256                     if (sqrt(pt_chk2) < _ptCutMin || sqrt(pt_chk2) > _ptCutMax)
1257                     {
1258                         accepted = false;
1259                         continue; // skips checking other daughter particles
1260                     }
1261                 }
1262             }
1263             if (_etaCutEnabled)
1264             {
1265                 for (int i = 0; i < 4; i++)
1266                 {
1267                     double eta_chk = pseudoRapidity(
1268                         decayVecs[i].GetPx(),
1269                         decayVecs[i].GetPy(),
1270                         decayVecs[i].GetPz()); // computes the pseudorapidity in the laboratory frame.
1271                     if (eta_chk < _etaCutMin || eta_chk > _etaCutMax)
1272                     { // if this particle does not fall into range
1273                         accepted = false;
1274                         continue; // skips checking other daughter particles
1275                     }
1276                 }
1277             }
1278             accepted = true;
1279             _nmbAccepted++;
1280 
1281         } while (!accepted);
1282         double md = getDaughterMass(ipid);
1283         if ((iFbadevent == 0) and (tcheck == 0))
1284         {
1285             // adds daughters as particles into the output event.
1286             for (unsigned int i = 0; i < 4; ++i)
1287             {
1288                 int sign = (i % 2 == 0) ? 1 : -1;
1289                 starlightParticle daughter(decayVecs[i].GetPx(),
1290                                            decayVecs[i].GetPy(),
1291                                            decayVecs[i].GetPz(),
1292                                            sqrt(decayVecs[i].GetPx() * decayVecs[i].GetPx() + decayVecs[i].GetPy() * decayVecs[i].GetPy() + decayVecs[i].GetPz() * decayVecs[i].GetPz() + md * md), // energy
1293                                            md,                                                                                                                                                      // _mass
1294                                            ipid * sign,                                                                                                                                             // make half of the particles pi^+, half pi^-
1295                                            sign);                                                                                                                                                   // charge
1296                 event.addParticle(daughter);
1297             }
1298         }
1299     }
1300 
1301     else
1302     {
1303         do
1304         {
1305             pickwEgamq2(comenergy, cmsEgamma, targetEgamma,
1306                         Q2, gamma_pz, gamma_pt, // photon infor in CMS frame
1307                         e_E, e_theta);          // electron info in target frame
1308             //
1309             momenta(comenergy, cmsEgamma, Q2, gamma_pz, gamma_pt, // input
1310                     rapidity, E, momx, momy, momz,                // VM
1311                     t_px, t_py, t_pz, t_E,                        // pomeron
1312                     e_phi, tcheck);                               // electron
1313             //
1314             // inelasticity: used for angular distributions
1315             double col_y = 1. - (e_E / _eEnergy) * std::pow(std::cos(e_theta / 2.), 2.);
1316             double col_polarization = (1 - col_y) / (1 - col_y + col_y * col_y / 2.);
1317             double spin_element = getSpinMatrixElement(comenergy, Q2, col_polarization);
1318             _nmbAttempts++;
1319 
1320             vmpid = ipid;
1321             // Two body dedcay in eSTARlight includes the angular corrections due to finite virtuality
1322             twoBodyDecay(ipid, comenergy, momx, momy, momz, spin_element,
1323                          px1, py1, pz1, px2, py2, pz2, iFbadevent);
1324             double pt1chk = sqrt(px1 * px1 + py1 * py1);
1325             double pt2chk = sqrt(px2 * px2 + py2 * py2);
1326             double eta1 = pseudoRapidity(px1, py1, pz1);
1327             double eta2 = pseudoRapidity(px2, py2, pz2);
1328 
1329             if (_ptCutEnabled && !_etaCutEnabled)
1330             {
1331                 if (pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax)
1332                 {
1333                     accepted = true;
1334                     _nmbAccepted++;
1335                 }
1336             }
1337             else if (!_ptCutEnabled && _etaCutEnabled)
1338             {
1339                 if (eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax)
1340                 {
1341                     accepted = true;
1342                     _nmbAccepted++;
1343                 }
1344             }
1345             else if (_ptCutEnabled && _etaCutEnabled)
1346             {
1347                 if (pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax)
1348                 {
1349                     if (eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax)
1350                     {
1351                         accepted = true;
1352                         _nmbAccepted++;
1353                     }
1354                 }
1355             }
1356             else if (!_ptCutEnabled && !_etaCutEnabled)
1357                 _nmbAccepted++;
1358         } while ((_ptCutEnabled || _etaCutEnabled) && !accepted);
1359     }
1360     if (iFbadevent == 0 && tcheck == 0)
1361     {
1362         int q1 = 0, q2 = 0;
1363         int ipid1, ipid2 = 0;
1364 
1365         double xtest = _randy.Rndom();
1366         if (xtest < 0.5)
1367         {
1368             q1 = 1;
1369             q2 = -1;
1370         }
1371         else
1372         {
1373             q1 = -1;
1374             q2 = 1;
1375         }
1376 
1377         if (ipid == 11 || ipid == 13)
1378         {
1379             ipid1 = -q1 * ipid;
1380             ipid2 = -q2 * ipid;
1381         }
1382         else if (ipid == 22)
1383         { // omega --> pi0 + gamma
1384             ipid1 = 22;
1385             ipid2 = 111;
1386         }
1387         else
1388         {
1389             ipid1 = q1 * ipid;
1390             ipid2 = q2 * ipid;
1391         }
1392 
1393         // - Outgoing electron - target frame
1394         double e_ptot = sqrt(e_E * e_E - starlightConstants::mel * starlightConstants::mel);
1395         double e_px = e_ptot * sin(e_theta) * cos(e_phi);
1396         double e_py = e_ptot * sin(e_theta) * sin(e_phi);
1397         double e_pz = e_ptot * cos(e_theta);
1398         lorentzVector electron(e_px, e_py, e_pz, e_E);
1399         event.addSourceElectron(electron);
1400         // - Generated photon - target frame
1401         double gamma_x = gamma_pt * cos(e_phi + starlightConstants::pi);
1402         double gamma_y = gamma_pt * sin(e_phi + starlightConstants::pi);
1403         lorentzVector gamma(gamma_x, gamma_y, gamma_pz, cmsEgamma);
1404         vector3 boostVector(0, 0, tanh(_rap_CM));
1405         (gamma).Boost(boostVector);
1406         event.addGamma(gamma, targetEgamma, Q2);
1407         // - Scattered target and transfered momenta at target vertex
1408         double target_pz = -_beamNucleus * sqrt(_pEnergy * _pEnergy - pow(starlightConstants::protonMass, 2.)) - t_pz;
1409         // 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
1410         lorentzVector target(-t_px, -t_py, target_pz, _beamNucleus * _pEnergy - t_E);
1411         double t_var = t_E * t_E - t_px * t_px - t_py * t_py - t_pz * t_pz;
1412         event.addScatteredTarget(target, t_var);
1413         // - Vector meson
1414         // - Vector meson
1415         if (_VMpidtest == starlightConstants::FOURPRONG)
1416         {
1417             return event; // Exit the method if FOURPRONG
1418         }
1419 
1420         // - Saving V.M. daughters
1421         double md = getDaughterMass(vmpid);
1422         bool isOmegaNeutralDecay = (vmpid == starlightConstants::PHOTON);
1423         if (isOmegaNeutralDecay)
1424         {
1425             // double mpi0 = starlightConstants::pionNeutralMass;
1426             double e_gamma1, px_gamma1, py_gamma1, pz_gamma1, e_gamma2, px_gamma2, py_gamma2, pz_gamma2;
1427             double Ed1 = sqrt(px1 * px1 + py1 * py1 + pz1 * pz1);
1428             starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
1429             event.addParticle(particle1);
1430             //
1431             pi0Decay(px2, py2, pz2, e_gamma1, px_gamma1, py_gamma1, pz_gamma1, e_gamma2, px_gamma2, py_gamma2, pz_gamma2, iFbadevent);
1432             starlightParticle particle2(px_gamma1, py_gamma1, pz_gamma1, e_gamma1, starlightConstants::UNKNOWN, ipid1, q2);
1433             starlightParticle particle3(px_gamma2, py_gamma2, pz_gamma2, e_gamma2, starlightConstants::UNKNOWN, ipid1, q2);
1434             event.addParticle(particle2);
1435             event.addParticle(particle3);
1436             // double invmass = sqrt(mpi0*mpi0 + 2.0*(Ed1*Ed2 - (px1*px2+py1*py2+pz1*pz2)));
1437         }
1438         else
1439         {
1440             double Ed1 = sqrt(md * md + px1 * px1 + py1 * py1 + pz1 * pz1);
1441             starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
1442             event.addParticle(particle1);
1443             //
1444             double Ed2 = sqrt(md * md + px2 * px2 + py2 * py2 + pz2 * pz2);
1445             starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
1446             event.addParticle(particle2);
1447         }
1448     }
1449     return event;
1450 }
1451 string Gammaavectormeson::gammaTableParse(int ii, int jj)
1452 {
1453     ostringstream tag1, tag2;
1454     tag1 << ii;
1455     tag2 << jj;
1456     string to_ret = tag1.str() + "," + tag2.str();
1457     return to_ret;
1458 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!