Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2024-09-27 07:03:40 

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 
0034 #include <iostream>
0035 #include <fstream>
0036 #include <cassert>
0037 #include <cmath>
0038 
0039 #include "gammaavm.h"
0040 //#include "photonNucleusCrossSection.h"
0041 #include "wideResonanceCrossSection.h"
0042 #include "narrowResonanceCrossSection.h"
0043 #include "incoherentVMCrossSection.h"
0044 //
0045 #include "e_narrowResonanceCrossSection.h"
0046 #include "e_wideResonanceCrossSection.h"
0047 
0048 using namespace std;
0049 
0050 
0051 //______________________________________________________________________________
0052 Gammaavectormeson::Gammaavectormeson(const inputParameters& inputParametersInstance, beamBeamSystem& bbsystem):eventChannel(inputParametersInstance, bbsystem), _phaseSpaceGen(0)
0053 {
0054     _VMNPT=inputParametersInstance.nmbPtBinsInterference();
0055     _VMWmax=inputParametersInstance.maxW();
0056     _VMWmin=inputParametersInstance.minW();
0057     //_VMYmax=inputParametersInstance.maxRapidity();
0058     //_VMYmin=-1.*_VMYmax;
0059     _VMnumw=inputParametersInstance.nmbWBins();
0060     //_VMnumy=inputParametersInstance.nmbRapidityBins();
0061     _VMnumega=inputParametersInstance.nmbEnergyBins();
0062     _VMnumQ2=inputParametersInstance.nmbGammaQ2Bins(); 
0063     _VMgamma_em=inputParametersInstance.beamLorentzGamma();
0064     _VMinterferencemode=inputParametersInstance.interferenceEnabled();
0065     _VMbslope=0.;//Will define in wide/narrow constructor
0066         _bslopeDef=inputParametersInstance.bslopeDefinition();
0067     _bslopeVal=inputParametersInstance.bslopeValue();
0068     _pEnergy= inputParametersInstance.protonEnergy();
0069     _beamNucleus = inputParametersInstance.targetBeamA();
0070     // electron energy in CMS frame
0071     _eEnergy= inputParametersInstance.electronEnergy();
0072     _VMpidtest=inputParametersInstance.prodParticleType();
0073     _VMptmax=inputParametersInstance.maxPtInterference();
0074     _VMdpt=inputParametersInstance.ptBinWidthInterference();
0075         _ProductionMode=inputParametersInstance.productionMode();
0076     _targetMaxPhotonEnergy=inputParametersInstance.targetMaxPhotonEnergy();
0077     _targetMinPhotonEnergy=inputParametersInstance.targetMinPhotonEnergy();
0078     // Now saving the photon energy limits
0079     _cmsMaxPhotonEnergy=inputParametersInstance.cmsMaxPhotonEnergy();
0080     _cmsMinPhotonEnergy=inputParametersInstance.cmsMinPhotonEnergy();
0081     _beamLorentzGamma = inputParametersInstance.beamLorentzGamma();
0082     _targetBeamLorentzGamma = inputParametersInstance.targetBeamLorentzGamma();
0083     _rap_CM=inputParametersInstance.rap_CM();
0084     _targetRadius = inputParametersInstance.targetRadius();
0085     //Turn on/off backwards production
0086     _backwardsProduction = inputParametersInstance.backwardsProduction();
0087     
0088         N0 = 0; N1 = 0; N2 = 0; 
0089     if (_VMpidtest == starlightConstants::FOURPRONG){
0090       // create n-body phase-spage generator
0091       _phaseSpaceGen = new nBodyPhaseSpaceGen(_randy);
0092     }
0093     if(_ProductionMode == 12 || _ProductionMode == 13) // Narrow and wide coherent photon-pommeron in eSTARlight
0094       _dummy_pncs = new photonNucleusCrossSection(inputParametersInstance, bbsystem);
0095     //
0096     const double r_04_00 = 0.674;
0097     const double cos_delta = 0.925;
0098     for( int ii = 0; ii < 100; ++ii){ //epsilon 0-1
0099       double epsilon = 0.01*ii;
0100       const double R = (1./epsilon)*r_04_00/(1.-r_04_00);
0101       for(int jj = 0; jj < 200; ++jj){ //psi 0 - 2pi
0102         double psi = jj*starlightConstants::pi/100.;
0103         double max_bin;
0104         for( int kk = 0; kk < 200; ++kk){ //temp
0105           double theta = kk*starlightConstants::pi/100.;
0106           // Fin max
0107           double this_test = std::pow(std::sin(theta),2.)*(1+epsilon*cos(2.*psi)) + 2.*epsilon*R*std::pow(std::cos(theta),2.)
0108         +std::sqrt(2.*epsilon*(1+epsilon))*cos_delta*std::sin(2.*theta)*std::cos(psi);
0109           if(this_test >  max_bin)
0110         max_bin = this_test;
0111         }
0112         _angular_max[ii][jj] = max_bin;
0113       }
0114     }
0115 }
0116 
0117 
0118 //______________________________________________________________________________
0119 Gammaavectormeson::~Gammaavectormeson()
0120 {
0121     if (_phaseSpaceGen)
0122         delete _phaseSpaceGen;
0123     if (_dummy_pncs)
0124       delete _dummy_pncs;
0125 }
0126 
0127 
0128 //______________________________________________________________________________
0129 void Gammaavectormeson::pickwy(double &W, double &Y)
0130 {
0131         double dW, dY, xw,xy,xtest;
0132     int  IW,IY;
0133   
0134     dW = (_VMWmax-_VMWmin)/double(_VMnumw);
0135     dY = (_VMYmax-_VMYmin)/double(_VMnumy);
0136   
0137  L201pwy:
0138 
0139     xw = _randy.Rndom();
0140     W = _VMWmin + xw*(_VMWmax-_VMWmin);
0141 
0142     if (W < 2 * starlightConstants::pionChargedMass)
0143         goto L201pwy;
0144   
0145     IW = int((W-_VMWmin)/dW);
0146     xy = _randy.Rndom();
0147     Y = _VMYmin + xy*(_VMYmax-_VMYmin);
0148     IY = int((Y-_VMYmin)/dY); 
0149     xtest = _randy.Rndom();
0150 
0151     if( xtest > _Farray[IW][IY] )
0152         goto L201pwy;
0153 
0154         N0++; 
0155     // Determine the target nucleus 
0156     // For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2 
0157     _TargetBeam=2; // Always true for eX
0158 }         
0159 
0160 
0161 
0162 //______________________________________________________________________________                                               
0163 void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
0164                                      double  W,
0165                                      double  px0, double  py0, double  pz0,
0166                      double  spin_element,
0167                                      double& px1, double& py1, double& pz1,
0168                                      double& px2, double& py2, double& pz2,
0169                                      int&    iFbadevent)
0170 {
0171     // This routine decays a particle into two particles of mass mdec,
0172     // taking spin into account
0173     double pmag;
0174     double phi,theta,Ecm;
0175     double betax,betay,betaz;
0176     double mdec1=0.0,mdec2=0.0;
0177     double E1=0.0,E2=0.0;
0178 
0179     //    set the mass of the daughter particles
0180     mdec1=getDaughterMass(ipid);
0181 
0182     if(_VMpidtest == starlightConstants::OMEGA_pi0gamma){
0183         mdec2=starlightConstants::pionNeutralMass;
0184         if(W < mdec2){
0185             cout<<" ERROR: W="<<W<<endl;
0186             iFbadevent = 1;
0187             return;
0188         }
0189         //  calculate the magnitude of the momenta
0190         pmag = (W*W - mdec2*mdec2)/(2.0*W);
0191     }else{
0192         mdec2=mdec1;
0193         if(W < 2*mdec1){
0194             cout<<" ERROR: W="<<W<<endl;
0195             iFbadevent = 1;
0196             return;
0197         }
0198         //  calculate the magnitude of the momenta
0199         pmag = sqrt(W*W/4. - mdec2*mdec2);
0200     }
0201       
0202     //  pick an orientation, based on the spin
0203     //  phi has a flat distribution in 2*pi
0204     phi = _randy.Rndom()*2.*starlightConstants::pi;
0205                                                                                                                 
0206     //  find theta, the angle between one of the outgoing particles and
0207     //  the beamline, in the frame of the two photons
0208     //cout<<"spin element: "<<spin_element<<endl;
0209     theta=getTheta(ipid, spin_element);
0210  
0211     //  compute unboosted momenta
0212     px1 = sin(theta)*cos(phi)*pmag;
0213     py1 = sin(theta)*sin(phi)*pmag;
0214     pz1 = cos(theta)*pmag;
0215     px2 = -px1;
0216     py2 = -py1;
0217     pz2 = -pz1;
0218 
0219     Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
0220     E1 = sqrt(mdec1*mdec1+px1*px1+py1*py1+pz1*pz1);
0221     E2 = sqrt(mdec2*mdec2+px2*px2+py2*py2+pz2*pz2);
0222 
0223     //cout<<"Ecm: "<<Ecm<<endl;
0224     //cout<<"W: "<<W<<endl;
0225 
0226     betax = -(px0/Ecm);
0227     betay = -(py0/Ecm);
0228     betaz = -(pz0/Ecm);
0229 
0230     transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
0231     transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
0232 
0233     if(iFbadevent == 1)
0234        return;
0235 
0236 }
0237 
0238 
0239 //______________________________________________________________________________                                               
0240 // decays a particle into four particles with isotropic angular distribution
0241 bool Gammaavectormeson::fourBodyDecay
0242 (starlightConstants::particleTypeEnum& ipid,
0243  const double                  ,           // E (unused)
0244  const double                  W,          // mass of produced particle
0245  const double*                 p,          // momentum of produced particle; expected to have size 3
0246  lorentzVector*                decayVecs,  // array of Lorentz vectors of daughter particles; expected to have size 4
0247  int&                          iFbadevent)
0248 {
0249     const double parentMass = W;
0250 
0251     // set the mass of the daughter particles
0252     const double daughterMass = getDaughterMass(ipid);
0253     if (parentMass < 4 * daughterMass){
0254         cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
0255         iFbadevent = 1;
0256         return false;
0257     }
0258 
0259     // construct parent four-vector
0260     const double        parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
0261                                             + parentMass * parentMass);
0262     const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
0263 
0264     // setup n-body phase-space generator
0265     assert(_phaseSpaceGen);
0266     static bool firstCall = true;
0267     if (firstCall) {
0268         const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
0269         _phaseSpaceGen->setDecay(4, m);
0270         // estimate maximum phase-space weight
0271         _phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
0272         firstCall = false;
0273     }
0274 
0275     // generate phase-space event
0276     if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
0277         return false;
0278 
0279     // set Lorentzvectors of decay daughters
0280     for (unsigned int i = 0; i < 4; ++i)
0281         decayVecs[i] = _phaseSpaceGen->daughter(i);
0282     return true;
0283 }
0284 
0285 //______________________________________________________________________________                                               
0286 void Gammaavectormeson::pi0Decay(double& px_pi0, double& py_pi0, double& pz_pi0,
0287                                  double& e_g1, double& px_g1, double& py_g1, double& pz_g1,
0288                                  double& e_g2, double& px_g2, double& py_g2, double& pz_g2,
0289                                  int&    iFbadevent)
0290 {
0291     double pmag;
0292     double phi,theta,Ecm;
0293     double betax,betay,betaz;
0294     double E1=0.0,E2=0.0;
0295 
0296     // This routine decays a pi0 into two isotropically produced photons
0297     double m_pi0=starlightConstants::pionNeutralMass;
0298     
0299     //  calculate the magnitude of the momenta
0300     pmag = sqrt(m_pi0*m_pi0/4.);
0301       
0302     //  pick an orientation, based on the spin
0303     //  phi has a flat distribution in 2*pi
0304     phi = _randy.Rndom()*2.*starlightConstants::pi;
0305                                                                                                                 
0306     //  find theta, the angle between one of the outgoing photons and
0307     //  the beamline, in the frame of the pi0
0308     theta=getTheta(starlightConstants::PION, 1.0/3.0);
0309  
0310     //  compute unboosted momenta
0311     px_g1 = sin(theta)*cos(phi)*pmag;
0312     py_g1 = sin(theta)*sin(phi)*pmag;
0313     pz_g1 = cos(theta)*pmag;
0314     px_g2 = -px_g1;
0315     py_g2 = -py_g1;
0316     pz_g2 = -pz_g1;
0317 
0318     Ecm = sqrt(m_pi0*m_pi0+px_pi0*px_pi0+py_pi0*py_pi0+pz_pi0*pz_pi0);
0319     E1 = sqrt(px_g1*px_g1+py_g1*py_g1+pz_g1*pz_g1);
0320     E2 = sqrt(px_g2*px_g2+py_g2*py_g2+pz_g2*pz_g2);
0321 
0322     betax = -(px_pi0/Ecm);
0323     betay = -(py_pi0/Ecm);
0324     betaz = -(pz_pi0/Ecm);
0325 
0326     transform (betax,betay,betaz,E1,px_g1,py_g1,pz_g1,iFbadevent);
0327     transform (betax,betay,betaz,E2,px_g2,py_g2,pz_g2,iFbadevent);
0328 
0329     e_g1=E1;
0330     e_g2=E2;
0331 
0332     if(iFbadevent == 1)
0333        return;
0334 }
0335 
0336 
0337 //______________________________________________________________________________
0338 double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
0339 {
0340     //This will return the daughter particles mass, and the final particles outputed id...
0341     double mdec=0.;
0342   
0343     switch(_VMpidtest){
0344     case starlightConstants::RHO:
0345     case starlightConstants::RHOZEUS:
0346     case starlightConstants::FOURPRONG:
0347     case starlightConstants::OMEGA:
0348         mdec = starlightConstants::pionChargedMass;
0349         ipid = starlightConstants::PION;
0350         break;
0351     case starlightConstants::OMEGA_pi0gamma:
0352         mdec = 0.0;
0353         ipid = starlightConstants::PHOTON;
0354         break;
0355     case starlightConstants::PHI:
0356         mdec = starlightConstants::kaonChargedMass;
0357         ipid = starlightConstants::KAONCHARGE;
0358         break;
0359     case starlightConstants::JPSI:
0360         mdec = starlightConstants::mel;
0361         ipid = starlightConstants::ELECTRON;
0362         break; 
0363     case starlightConstants::JPSI_ee:
0364         mdec = starlightConstants::mel;
0365         ipid = starlightConstants::ELECTRON;
0366         break; 
0367     case starlightConstants::JPSI_mumu:
0368         mdec = starlightConstants::muonMass;
0369         ipid = starlightConstants::MUON;
0370         break; 
0371     case starlightConstants::JPSI2S_ee:
0372         mdec = starlightConstants::mel;
0373         ipid = starlightConstants::ELECTRON;
0374         break; 
0375     case starlightConstants::JPSI2S_mumu:
0376         mdec = starlightConstants::muonMass;
0377         ipid = starlightConstants::MUON;
0378         break; 
0379 
0380     case starlightConstants::JPSI2S:
0381     case starlightConstants::UPSILON:
0382     case starlightConstants::UPSILON2S:
0383     case starlightConstants::UPSILON3S:
0384         mdec = starlightConstants::muonMass;
0385         ipid = starlightConstants::MUON;
0386         break;
0387     case starlightConstants::UPSILON_ee:
0388     case starlightConstants::UPSILON2S_ee:
0389     case starlightConstants::UPSILON3S_ee:
0390         mdec = starlightConstants::mel;
0391         ipid = starlightConstants::ELECTRON;
0392         break;
0393     case starlightConstants::UPSILON_mumu:
0394     case starlightConstants::UPSILON2S_mumu:
0395     case starlightConstants::UPSILON3S_mumu:
0396         mdec = starlightConstants::muonMass;
0397         ipid = starlightConstants::MUON;   
0398         break;
0399     default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
0400     }
0401   
0402     return mdec;
0403 }
0404 
0405 
0406 
0407 //______________________________________________________________________________
0408 double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid, double r_04_00)
0409 {
0410     //This depends on the decay angular distribution
0411     //Valid for rho, phi, omega.
0412     double theta=0.;
0413     double xtest=1.;
0414     double dndtheta=0.;
0415 
0416     while(xtest > dndtheta){
0417     
0418       theta = starlightConstants::pi*_randy.Rndom();
0419       xtest = _randy.Rndom();
0420       //  Follow distribution for helicity +/-1 with finite Q2
0421       //  Physics Letters B 449, 328; The European Physical Journal C - Particles and Fields 13, 37; 
0422       //  The European Physical Journal C - Particles and Fields 6, 603
0423   
0424 
0425       switch(ipid){
0426         
0427       case starlightConstants::MUON:
0428       case starlightConstants::ELECTRON:
0429         //primarily for upsilon/j/psi.  VM->ee/mumu
0430         dndtheta = sin(theta)*(1 + r_04_00+( 1-3.*r_04_00 )*cos(theta)*cos(theta));
0431         break;
0432         
0433       case starlightConstants::PION:
0434       case starlightConstants::KAONCHARGE:
0435         //rhos etc
0436         dndtheta=  sin(theta)*(1 - r_04_00+( 3.*r_04_00-1 )*cos(theta)*cos(theta));
0437         break;
0438         
0439       default: if(!_backwardsProduction) cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
0440       }//end of switch
0441 
0442       // Assume unpolarized vector-mesons for now
0443       if(_backwardsProduction){
0444         dndtheta = sin(theta);
0445       }
0446     }
0447     return theta;
0448 
0449 }
0450 
0451 
0452 //______________________________________________________________________________
0453 double Gammaavectormeson::getSpinMatrixElement(double W, double Q2, double epsilon)
0454 {
0455   double m2 = 0.6*(W*W);
0456   double R = starlightConstants::pi/2.*m2/Q2 - std::pow(m2,3./2.)/sqrt(Q2)/(Q2+m2) - m2/Q2*atan(sqrt(m2/Q2));
0457   //////////needed for problem with double precision
0458   double R1 = starlightConstants::pi/2.*m2/Q2;
0459   if(R1>1.0E10)R=0.0;
0460   ////////// 
0461   R = (Q2+m2)*(Q2+m2)*R*R/m2;
0462   return epsilon*R/(1+epsilon*R);
0463 }
0464 
0465 
0466 //______________________________________________________________________________
0467 double Gammaavectormeson::getWidth()
0468 {
0469     return _width;
0470 }
0471 
0472 
0473 //______________________________________________________________________________
0474 double Gammaavectormeson::getMass()
0475 {
0476     return _mass;
0477 }
0478 
0479 
0480 //______________________________________________________________________________
0481 double Gammaavectormeson::getSpin()
0482 {
0483     return 1.0; //VM spins are the same
0484 }
0485 
0486 //______________________________________________________________________________
0487 void Gammaavectormeson::momenta(double W,double Egam,double Q2, double gamma_pz, double gamma_pt, //input conditions
0488                 double &Y,double &E,double &px,double &py,double &pz,  //return vm
0489                 double &t_px, double &t_py, double &t_pz, double &t_E, //return pomeron
0490                 double &e_phi,int &tcheck) //return electron (angle already known by Q2)
0491 {
0492     //     This subroutine calculates momentum and energy of vector meson for electroproduction (eSTARlight)
0493     //     given W and photon 4-vector,   without interference.  No intereference in asymetric eX collisions
0494  
0495     double Epom,Pom_pz,tmin,pt2,phi1,phi2;
0496     double px1,py1;
0497     double xt,xtest,ytest;
0498     double t2;
0499 
0500     double target_px, target_py, target_pz, target_E;
0501 
0502     target_E = _beamNucleus*_pEnergy;
0503     target_px = 0.0;
0504     target_py = 0.0;
0505     target_pz = -_beamNucleus*sqrt(_pEnergy*_pEnergy - pow(starlightConstants::protonMass,2.) );
0506     phi1 = 2.*starlightConstants::pi*_randy.Rndom();
0507     px1 = gamma_pt*cos(phi1);
0508     py1 = gamma_pt*sin(phi1);
0509     int isbadevent = 0;
0510     double betax_cm = ((px1+target_px)/(Egam+target_E));
0511     double betay_cm = ((py1+target_py)/(Egam+target_E));
0512     double betaz_cm = ((gamma_pz+target_pz)/(Egam+target_E));
0513     transform (betax_cm,betay_cm,betaz_cm,target_E,target_px,target_py,target_pz,isbadevent);
0514     transform (betax_cm,betay_cm,betaz_cm,Egam, px1, py1, gamma_pz, isbadevent);
0515       
0516     e_phi = starlightConstants::pi+phi1;
0517     // Pomeron pz is != than its energy in eSTARlight, in order to conserve energy/momentum of scattered
0518     // target
0519     //Pom_pz = 0.5*(W*W-Q2)/(Egam + gamma_pz);
0520     Epom = 0.5*(W*W+Q2)/(Egam + target_E);
0521 
0522     L522vm:
0523     while( e_phi > 2.*starlightConstants::pi ) e_phi-= 2.*starlightConstants::pi;
0524     //
0525     if ( ( _bbs.targetBeam().A()==1 ) 
0526           || (_ProductionMode == 4) ) {
0527         if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
0528           // Use dipole form factor for light VM
0529         L613vm:
0530           xtest = 2.0*_randy.Rndom();
0531               double ttest = xtest*xtest; 
0532               ytest = _randy.Rndom();
0533               double t0 = 1./2.23; 
0534               double yprob = xtest*_bbs.electronBeam().dipoleFormFactor(ttest,t0)*_bbs.electronBeam().dipoleFormFactor(ttest,t0); 
0535               if( ytest > yprob ) goto L613vm; 
0536               t2 = ttest; 
0537               pt2 = xtest;              
0538         }else{
0539         //Use dsig/dt= exp(-_VMbslope*t) for heavy VM
0540                 double bslope_tdist = _VMbslope; 
0541         double Wgammap = 0.0; 
0542                 switch(_bslopeDef){
0543           case 0:
0544             //This is the default, as before
0545             bslope_tdist = _VMbslope;
0546             break;
0547           case 1:
0548             //User defined value of bslope. BSLOPE_VALUE default is 4.0 if not set. 
0549                     bslope_tdist = _bslopeVal;
0550             if( N0 <= 1 )cout<<" ATTENTION: Using user defined value of bslope = "<<_bslopeVal<<endl;
0551                     break; 
0552           case 2:
0553                     //This is Wgammap dependence of b from H1 (Eur. Phys. J. C 46 (2006) 585)
0554             Wgammap = sqrt(4.*Egam*_pEnergy); 
0555             bslope_tdist = 4.63 + 4.*0.164*log(Wgammap/90.0);
0556             if( N0 <= 1 )cout<<" ATTENTION: Using energy dependent value of bslope!"<<endl; 
0557             break;
0558           default:
0559             cout<<" Undefined setting for BSLOPE_DEFINITION "<<endl;
0560         }
0561 
0562             xtest = _randy.Rndom(); 
0563         // t2 = (-1./_VMbslope)*log(xtest);
0564         t2 = (-1./bslope_tdist)*log(xtest);
0565         pt2 = sqrt(1.*t2);
0566         }
0567     } else {
0568         // >> Check tmin
0569         tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
0570 
0571         if(tmin > 0.5){
0572         cout<<" WARNING: tmin= "<<tmin<<endl;
0573                 cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl; 
0574         cout<<" Will pick a new W,Y "<<endl;
0575         tcheck = 1;
0576         return;
0577         }
0578  L663vm:
0579         xt = _randy.Rndom(); 
0580             if( _bbs.targetBeam().A() != 1){ 
0581           pt2 = 8.*xt*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius();
0582         } 
0583         else{
0584           std::cout<<"Can't find the electron for eX"<<std::endl;
0585         }  
0586 
0587         xtest = _randy.Rndom();
0588         t2 = tmin + pt2*pt2;
0589 
0590         double comp=0.0; 
0591             if( _bbs.targetBeam().A() != 1){ 
0592           comp = _bbs.targetBeam().formFactor(t2)*_bbs.targetBeam().formFactor(t2)*pt2;
0593         }
0594         else 
0595           std::cout<<"Can't find the electron for eX"<<std::endl;
0596             if( xtest > comp ) goto L663vm;
0597             
0598     }
0599     phi2 = 2.*starlightConstants::pi*_randy.Rndom();
0600 
0601     //
0602     t_px = pt2*cos(phi2);
0603     t_py = pt2*sin(phi2);
0604     // Used to return the pomeron pz to generator
0605     // Compute scattered target kinematics p_(initial ion) = p_(pomeron) + p_(final ion)
0606     double newion_E = target_E-Epom;
0607     double newion_px = target_px - t_px;
0608     double newion_py = target_py - t_py;
0609     double newion_pz_squared = newion_E*newion_E - newion_px*newion_px - newion_py*newion_py - pow(_beamNucleus*starlightConstants::protonMass,2.);
0610     if(newion_pz_squared<0)goto L522vm;
0611     double newion_pz = -sqrt( newion_pz_squared );
0612     Pom_pz = target_pz - newion_pz;
0613     t_pz = Pom_pz;
0614     // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
0615     px = px1 + t_px;
0616     py = py1 + t_py;
0617 
0618     t_E = Epom;
0619     // Finally V.M. energy, pz and rapidity from photon + pommeron.
0620     E = Egam + t_E;
0621     pz = gamma_pz + t_pz;
0622     //transform back to e-ion CM frame
0623     transform (-betax_cm,-betay_cm,-betaz_cm,target_E,target_px,target_py,target_pz,isbadevent);
0624     transform (-betax_cm,-betay_cm,-betaz_cm,newion_E,newion_px,newion_py,newion_pz,isbadevent);
0625     transform (-betax_cm,-betay_cm,-betaz_cm,Egam,    px1,      py1,      gamma_pz, isbadevent);
0626     transform (-betax_cm,-betay_cm,-betaz_cm,E,       px,       py,       pz,       isbadevent);
0627     t_px = target_px-newion_px;
0628     t_py = target_py-newion_py;
0629     t_pz = target_pz-newion_pz;
0630     t_E  = target_E-newion_E;
0631     Y = 0.5*std::log( (E+fabs(pz))/(E-fabs(pz)) );
0632     //  cout << "pz_final = " << pz << endl;
0633 
0634 
0635     // Backwards Production... updated by Zachary Sweger on 9/21/2021
0636     if (_backwardsProduction) {
0637       double is_proton_E  = target_E;
0638       double is_proton_px = target_px;
0639       double is_proton_py = target_py;
0640       double is_proton_pz = target_pz;
0641       double is_gamma_E  = Egam;
0642       double is_gamma_px = px1;
0643       double is_gamma_py = py1;
0644       double is_gamma_pz = gamma_pz;
0645       double is_tot_E  = is_proton_E  + is_gamma_E;
0646       double is_tot_px = is_proton_px + is_gamma_px;
0647       double is_tot_py = is_proton_py + is_gamma_py;
0648       double is_tot_pz = is_proton_pz + is_gamma_pz;
0649 
0650       double fs_proton_E  = newion_E;
0651       double fs_proton_px = newion_px;
0652       double fs_proton_py = newion_py;
0653       double fs_proton_pz = newion_pz;
0654       
0655       double fs_vm_E  = E;
0656       double fs_vm_px = px;
0657       double fs_vm_py = py;
0658       double fs_vm_pz = pz;
0659 
0660       int isbadevent = 0;
0661       double betax_cm = (is_tot_px/is_tot_E);
0662       double betay_cm = (is_tot_py/is_tot_E);
0663       double betaz_cm = (is_tot_pz/is_tot_E);
0664       transform (betax_cm,betay_cm,betaz_cm,is_proton_E,is_proton_px,is_proton_py,is_proton_pz,isbadevent);
0665       transform (betax_cm,betay_cm,betaz_cm,is_gamma_E, is_gamma_px, is_gamma_py, is_gamma_pz, isbadevent);
0666       transform (betax_cm,betay_cm,betaz_cm,fs_proton_E,fs_proton_px,fs_proton_py,fs_proton_pz,isbadevent);
0667       transform (betax_cm,betay_cm,betaz_cm,fs_vm_E,    fs_vm_px,    fs_vm_py,    fs_vm_pz,    isbadevent);
0668 
0669       double u_fs_proton_px = -1.0*fs_proton_px;
0670       double u_fs_proton_py = -1.0*fs_proton_py;
0671       double u_fs_proton_pz = -1.0*fs_proton_pz;
0672       double u_fs_proton_E  = fs_proton_E;
0673       
0674       double u_fs_vm_px = -1.0*fs_vm_px;
0675       double u_fs_vm_py = -1.0*fs_vm_py;
0676       double u_fs_vm_pz = -1.0*fs_vm_pz;
0677       double u_fs_vm_E  = fs_vm_E;
0678 
0679       betax_cm = -1.0*betax_cm;
0680       betay_cm = -1.0*betay_cm;
0681       betaz_cm = -1.0*betaz_cm;
0682       transform (betax_cm,betay_cm,betaz_cm,is_proton_E,  is_proton_px,  is_proton_py,  is_proton_pz,  isbadevent);
0683       transform (betax_cm,betay_cm,betaz_cm,is_gamma_E,   is_gamma_px,   is_gamma_py,   is_gamma_pz,   isbadevent);
0684       transform (betax_cm,betay_cm,betaz_cm,u_fs_proton_E,u_fs_proton_px,u_fs_proton_py,u_fs_proton_pz,isbadevent);
0685       transform (betax_cm,betay_cm,betaz_cm,u_fs_vm_E,    u_fs_vm_px,    u_fs_vm_py,    u_fs_vm_pz,    isbadevent);
0686 
0687       px=u_fs_vm_px;
0688       py=u_fs_vm_py;
0689       pz=u_fs_vm_pz;
0690       E =u_fs_vm_E;
0691 
0692       t_px=-u_fs_proton_px;
0693       t_py=-u_fs_proton_py;
0694       t_pz=is_proton_pz-u_fs_proton_pz;
0695       t_E =is_proton_E-u_fs_proton_E;
0696       //cout<<"Energy Violation: "<<(E+u_fs_proton_E)-(is_tot_E)<<endl;
0697       //cout<<"Pz Violation: "<<(pz+u_fs_proton_pz)-(is_tot_pz)<<endl;
0698       //cout<<"Px Violation: "<<(px+u_fs_proton_px)-(is_tot_px)<<endl;
0699       //cout<<"Py Violation: "<<(py+u_fs_proton_py)-(is_tot_py)<<endl;
0700       //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;
0701       //cout << "VM MASS = " << sqrt(E*E-px*px - py*py-pz*pz) << endl;
0702 
0703       //Calculate the rapidity
0704       Y = 0.5*std::log( (E+fabs(pz))/(E-fabs(pz)) );
0705       if(Y!=Y){
0706         //FIXME the following should be upated if not using omegas
0707         if(_VMpidtest==starlightConstants::OMEGA_pi0gamma || _VMpidtest==starlightConstants::OMEGA) E=sqrt(starlightConstants::OmegaMass*starlightConstants::OmegaMass+px*px+py*py+pz*pz); 
0708         else if(_VMpidtest==starlightConstants::RHO) E=sqrt(starlightConstants::rho0Mass*starlightConstants::rho0Mass+px*px+py*py+pz*pz);
0709         //cout<<"Energy Violation: "<<(E+u_fs_proton_E)-(is_tot_E)<<endl;
0710       }
0711     }
0712     
0713 }
0714 
0715 
0716 //______________________________________________________________________________
0717 double Gammaavectormeson::pTgamma(double E)
0718 {
0719     // returns on random draw from pp(E) distribution
0720     double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
0721     double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
0722     int satisfy =0;
0723         
0724     ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
0725     //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
0726     Cm = sqrt(3.)*E/_VMgamma_em;
0727     // If E is very small, the drawing of a pT below is extremely slow. 
0728     // ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E. 
0729     // Should have no observable consequences (JN, SRK 11-Sep-2014)
0730     if( E < 0.0005 )return Cm; 
0731  
0732     //the amplitude of the p_t spectrum at the maximum
0733 
0734     if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A() != 1){ 
0735       if( _ProductionMode == 2 || _ProductionMode ==3 ){
0736          singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
0737       }else{
0738          singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
0739       }  
0740     } else if( _bbs.targetBeam().A()==1 && _bbs.electronBeam().A() != 1 ){
0741       if( _ProductionMode == 2 || _ProductionMode ==3){
0742          singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
0743       }else{
0744          singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
0745       }  
0746     } else if (_TargetBeam == 1) {
0747       singleformfactorCm=_bbs.targetBeam().formFactor(Cm*Cm+ereds);
0748     } else {
0749       singleformfactorCm=_bbs.electronBeam().formFactor(Cm*Cm+ereds);
0750     }
0751 
0752     Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
0753         
0754     //pick a test value pp, and find the amplitude there
0755     x = _randy.Rndom();
0756 
0757     if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A() != 1){ 
0758       if( _ProductionMode == 2 || _ProductionMode ==3){
0759         pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); 
0760         singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
0761       }else{
0762         pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); 
0763         singleformfactorpp1=_bbs.targetBeam().formFactor(pp*pp+ereds);
0764       }  
0765     } else if( _bbs.targetBeam().A()==1 && _bbs.electronBeam().A() != 1 ){
0766       if( _ProductionMode == 2 || _ProductionMode ==3){
0767         pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); 
0768         singleformfactorpp1=_bbs.targetBeam().formFactor(pp*pp+ereds);
0769       }else{
0770         pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); 
0771         singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
0772       }  
0773     } else if (_TargetBeam == 1) {
0774         pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); 
0775         singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
0776     } else {
0777         pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); 
0778         singleformfactorpp1=_bbs.electronBeam().formFactor(pp*pp+ereds);
0779     }
0780 
0781     test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
0782 
0783     while(satisfy==0){
0784     u = _randy.Rndom();
0785     if(u*Coef <= test)
0786     {
0787         satisfy =1;
0788     }
0789     else{
0790         x =_randy.Rndom();
0791             if( _bbs.electronBeam().A()==1 && _bbs.targetBeam().A() != 1){ 
0792               if( _ProductionMode == 2 || _ProductionMode ==3){
0793                 pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); 
0794                 singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
0795               }else{
0796                 pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); 
0797                 singleformfactorpp2=_bbs.targetBeam().formFactor(pp*pp+ereds);
0798               }  
0799             } else if( _bbs.targetBeam().A()==1 && _bbs.electronBeam().A() != 1 ){
0800               if( _ProductionMode == 2 || _ProductionMode ==3){
0801                 pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); 
0802                 singleformfactorpp2=_bbs.targetBeam().formFactor(pp*pp+ereds);
0803               }else{
0804                 pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); 
0805                 singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
0806               }  
0807             } else if (_TargetBeam == 1) {
0808               pp = x*5.*starlightConstants::hbarc/_bbs.targetBeam().nuclearRadius(); 
0809               singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
0810             } else {
0811               pp = x*5.*starlightConstants::hbarc/_bbs.electronBeam().nuclearRadius(); 
0812               singleformfactorpp2=_bbs.electronBeam().formFactor(pp*pp+ereds);
0813             }
0814         test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
0815     }
0816     }
0817 
0818     return pp;
0819 }
0820 
0821 
0822 //______________________________________________________________________________
0823 void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
0824                              int&) // tcheck (unused)
0825 {
0826     //    This function calculates momentum and energy of vector meson
0827     //    given W and Y, including interference.
0828     //    It gets the pt distribution from a lookup table.
0829     double dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
0830     int IY=0,j=0;
0831   
0832     dY  = (_VMYmax-_VMYmin)/double(_VMnumy);
0833   
0834     //  Y is already fixed; choose a pt
0835     //  Follow the approach in pickwy
0836     //  in  _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
0837     // Changed,  now works -y to +y.
0838     IY=int((Y-_VMYmin)/dY);
0839     if (IY > (_VMnumy)-1){
0840             IY=(_VMnumy)-1;
0841     }
0842 
0843     yleft=(Y-_VMYmin)-(IY)*dY;
0844 
0845     yfract=yleft*dY;
0846   
0847     xpt=_randy.Rndom();
0848     for(j=0;j<_VMNPT+1;j++){
0849         if (xpt < _fptarray[IY][j]) goto L60;
0850     }
0851  L60:
0852   
0853     //  now do linear interpolation - start with extremes
0854     if (j == 0){
0855         pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
0856         goto L80;
0857     }
0858     if (j == _VMNPT){
0859         pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
0860         goto L80;
0861     }
0862   
0863     //  we're in the middle
0864     ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
0865     pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
0866   
0867     //  at an extreme in y?
0868     if (IY == (_VMnumy/2)-1){
0869         pt=pt1;
0870         goto L120;
0871     }
0872  L80:
0873 
0874     //  interpolate in y repeat for next fractional y bin      
0875     for(j=0;j<_VMNPT+1;j++){
0876         if (xpt < _fptarray[IY+1][j]) goto L90;
0877     }
0878  L90:
0879   
0880     //  now do linear interpolation - start with extremes
0881     if (j == 0){
0882         pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
0883         goto L100;
0884     }
0885     if (j == _VMNPT){
0886         pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
0887         goto L100;
0888     }
0889   
0890     //  we're in the middle
0891     ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
0892     pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
0893  L100:
0894 
0895     //  now interpolate in y  
0896     pt=yfract*pt2+(1-yfract)*pt1;
0897  L120:
0898 
0899     //  we have a pt 
0900     theta=2.*starlightConstants::pi*_randy.Rndom();
0901     px=pt*cos(theta);
0902     py=pt*sin(theta);
0903 
0904     E  = sqrt(W*W+pt*pt)*cosh(Y);
0905     pz = sqrt(W*W+pt*pt)*sinh(Y);
0906     //      randomly choose to make pz negative 50% of the time
0907     if(_randy.Rndom()>=0.5) pz = -pz;
0908 }
0909 
0910 
0911 
0912 //______________________________________________________________________________
0913 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
0914 {
0915     double pT = sqrt(px*px + py*py);
0916     double p = sqrt(pz*pz + pT*pT);
0917     double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
0918     return eta;
0919 }
0920 
0921 //______________________________________________________________________________
0922 Gammaanarrowvm::Gammaanarrowvm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
0923 {
0924     cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
0925     read();
0926     cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
0927     narrowResonanceCrossSection sigma(input, bbsystem);
0928     sigma.crossSectionCalculation(_bwnormsave);
0929     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
0930     _VMbslope=sigma.slopeParameter(); 
0931 }
0932 
0933 
0934 //______________________________________________________________________________
0935 Gammaanarrowvm::~Gammaanarrowvm()
0936 { }
0937 
0938 
0939 //______________________________________________________________________________
0940 Gammaaincoherentvm::Gammaaincoherentvm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
0941 {
0942         cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
0943         read();
0944         cout<<"Creating and calculating crosssection. Gammaaincoherentvm()"<<endl;
0945         incoherentVMCrossSection sigma(input, bbsystem);
0946     sigma.crossSectionCalculation(_bwnormsave);
0947     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
0948         _VMbslope=sigma.slopeParameter(); 
0949 }
0950 
0951 
0952 //______________________________________________________________________________
0953 Gammaaincoherentvm::~Gammaaincoherentvm()
0954 { }
0955 
0956 
0957 //______________________________________________________________________________
0958 Gammaawidevm::Gammaawidevm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
0959 {
0960     cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
0961     read();
0962     cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
0963     wideResonanceCrossSection sigma(input, bbsystem);
0964     sigma.crossSectionCalculation(_bwnormsave);
0965     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
0966     _VMbslope=sigma.slopeParameter();
0967 }
0968 
0969 
0970 //______________________________________________________________________________
0971 Gammaawidevm::~Gammaawidevm()
0972 { }
0973 
0974 
0975 //______________________________________________________________________________
0976 e_Gammaanarrowvm::e_Gammaanarrowvm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
0977 {
0978     cout<<"Reading in luminosity tables. e_Gammaanarrowvm()"<<endl;
0979     e_read();
0980     cout<<"Creating and calculating crosssection. e_Gammaanarrowvm()"<<endl;
0981     e_narrowResonanceCrossSection sigma(input, bbsystem);
0982     sigma.crossSectionCalculation(_bwnormsave);
0983     setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
0984     _VMbslope=sigma.slopeParameter(); 
0985 }
0986 
0987 
0988 //______________________________________________________________________________
0989 e_Gammaanarrowvm::~e_Gammaanarrowvm()
0990 { }
0991 
0992 
0993 //______________________________________________________________________________
0994 e_Gammaawidevm::e_Gammaawidevm(const inputParameters& input, beamBeamSystem& bbsystem):Gammaavectormeson(input, bbsystem)
0995 {
0996     cout<<"Reading in luminosity tables. e_Gammaawidevm()"<<endl;
0997     e_read();
0998     cout<<"Creating and calculating crosssection. e_Gammaawidevm()"<<endl;
0999     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 }
[ 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!