Index: trunk/src/photonNucleusCrossSection.cpp
===================================================================
--- trunk/src/photonNucleusCrossSection.cpp	(revision 305)
+++ trunk/src/photonNucleusCrossSection.cpp	(revision 306)
@@ -1,925 +1,925 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 //    Copyright 2010
 //
 //    This file is part of starlight.
 //
 //    starlight is free software: you can redistribute it and/or modify
 //    it under the terms of the GNU General Public License as published by
 //    the Free Software Foundation, either version 3 of the License, or
 //    (at your option) any later version.
 //
 //    starlight is distributed in the hope that it will be useful,
 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 //    GNU General Public License for more details.
 //
 //    You should have received a copy of the GNU General Public License
 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
 //
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
 // $Rev::                             $: revision of last commit
 // $Author::                          $: author of last commit
 // $Date::                            $: date of last commit
 //
 // Description:
 //
 //
 ///////////////////////////////////////////////////////////////////////////
 
 
 #include <iostream>
 #include <fstream>
 #include <cmath>
 
 #include "reportingUtils.h"
 #include "starlightconstants.h"
 #include "bessel.h"
 #include "photonNucleusCrossSection.h"
 
 
 using namespace std;
 using namespace starlightConstants;
 
 
 //______________________________________________________________________________
 photonNucleusCrossSection::photonNucleusCrossSection(const inputParameters& inputParametersInstance, const beamBeamSystem& bbsystem)
         : _ip(&inputParametersInstance),
           _nWbins            (inputParametersInstance.nmbWBins()          ),
 	  _nYbins            (inputParametersInstance.nmbRapidityBins()   ),
 	  _wMin              (inputParametersInstance.minW()              ),
 	  _wMax              (inputParametersInstance.maxW()              ),
 	  _yMax              (inputParametersInstance.maxRapidity()       ),
 	  _beamLorentzGamma  (inputParametersInstance.beamLorentzGamma()  ),
 	  _bbs               (bbsystem                                    ),
 	  _protonEnergy      (inputParametersInstance.protonEnergy()      ),
 	  _particleType      (inputParametersInstance.prodParticleType()  ),
 	  _beamBreakupMode   (inputParametersInstance.beamBreakupMode()   ),
           _productionMode    (inputParametersInstance.productionMode()    ),
 	  _sigmaNucleus      (_bbs.beam2().A()          )
 {
   // new options - impulse aproximation (per Joakim) and Quantum Glauber (per SK) SKQG
          _impulseSelected = inputParametersInstance.impulseVM();
 	 _quantumGlauber = inputParametersInstance.quantumGlauber();
 	switch(_particleType) {
 	case RHO:
 		_slopeParameter = 11.0;  // [(GeV/c)^{-2}]
 		_vmPhotonCoupling = 2.02;
 		_ANORM       = -2.75;
 		_BNORM       = 0.0;
 		_defaultC    = 1.0;
 		_channelMass = _ip->rho0Mass(); 
 		_width       = _ip->rho0Width(); 
 		break;
 	case RHOZEUS:
 		_slopeParameter =11.0;
 		_vmPhotonCoupling=2.02;
 		_ANORM=-2.75;
 		_BNORM=1.84;
 		_defaultC=1.0;
 		_channelMass  = _ip->rho0Mass();
 		_width        = _ip->rho0Width();
 		break;
 	case FOURPRONG:
 		_slopeParameter      = 11.0;
 		_vmPhotonCoupling      = 2.02;
 		_ANORM       = -2.75;
 		_BNORM       = 0;  
-		_defaultC    = 11.0;
+		_defaultC    = 1.0;
 		_channelMass  = _ip->rho0PrimeMass();
 		_width        = _ip->rho0PrimeWidth();
 		break;
 	case OMEGA:
 		_slopeParameter=10.0;
 		_vmPhotonCoupling=23.13;
 		_ANORM=-2.75;
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->OmegaMass();
 		_width        = _ip->OmegaWidth();
 		break;
 	case PHI:
 		_slopeParameter=7.0;
 		_vmPhotonCoupling=13.71;
 		_ANORM=-2.75;
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->PhiMass();
 		_width        = _ip->PhiWidth();
 		break;
 	case JPSI:
 	case JPSI_ee:
 	case JPSI_mumu:
 	case JPSI_ppbar:
 		_slopeParameter=4.0;
 		_vmPhotonCoupling=10.45;
 		_ANORM=-2.75; 
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->JpsiMass(); 
 		_width        = _ip->JpsiWidth(); 
 		break;
 	case JPSI2S:
 	case JPSI2S_ee:
 	case JPSI2S_mumu:
 		_slopeParameter=4.3;
 		_vmPhotonCoupling=26.39;
 		_ANORM=-2.75; 
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->Psi2SMass();
 		_width        = _ip->Psi2SWidth();
 		break;
 	case UPSILON:
 	case UPSILON_ee:
 	case UPSILON_mumu:
 		_slopeParameter=4.0;
 		_vmPhotonCoupling=125.37;
 		_ANORM=-2.75; 
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->Upsilon1SMass();
 		_width        = _ip->Upsilon1SWidth();
 		break;
 	case UPSILON2S:
 	case UPSILON2S_ee:
 	case UPSILON2S_mumu:
 		_slopeParameter=4.0;
 		_vmPhotonCoupling=290.84;
 		_ANORM=-2.75;
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->Upsilon2SMass();
 		_width        = _ip->Upsilon2SWidth();		
 		break;
 	case UPSILON3S:
 	case UPSILON3S_ee:
 	case UPSILON3S_mumu:
 		_slopeParameter=4.0;
 		_vmPhotonCoupling=415.10;
 		_ANORM=-2.75;
 		_BNORM=0.0;
 		_defaultC=1.0;
 		_channelMass  = _ip->Upsilon3SMass();
 		_width        = _ip->Upsilon3SWidth();
 		break;
 	default:
 		cout <<"No sigma constants parameterized for pid: "<<_particleType
 		     <<" GammaAcrosssection"<<endl;
 	}
 
 	_maxPhotonEnergy = 12. * _beamLorentzGamma * hbarc/(_bbs.beam1().nuclearRadius()+_bbs.beam2().nuclearRadius()); 
 }
 
 
 //______________________________________________________________________________
 photonNucleusCrossSection::~photonNucleusCrossSection()
 { }
 
 
 //______________________________________________________________________________
 void
 photonNucleusCrossSection::crossSectionCalculation(const double)
 {
 	cout << "Neither narrow/wide resonance cross-section calculation.--Derived" << endl;
 }
 
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::getcsgA(const double Egamma,
                                    const double W, 
                                    const int beam)
 {
 	//This function returns the cross-section for photon-nucleus interaction 
 	//producing vectormesons
   
 	double Av,Wgp,cs,cvma;
 	double t,tmin,tmax;
 	double csgA,ax,bx;
 	int NGAUSS; 
   
 	//     DATA FOR GAUSS INTEGRATION
 	double xg[6] = {0, 0.1488743390, 0.4333953941, 0.6794095683, 0.8650633667, 0.9739065285};
 	double ag[6] = {0, 0.2955242247, 0.2692667193, 0.2190863625, 0.1494513492, 0.0666713443};
 	NGAUSS = 6;
   
 	//       Find gamma-proton CM energy
 	Wgp = sqrt(2. * Egamma * (_protonEnergy
 	                          + sqrt(_protonEnergy * _protonEnergy - _ip->protonMass() * _ip->protonMass()))
 	           + _ip->protonMass() * _ip->protonMass());
 	
 	//Used for A-A
 	tmin = (W * W / (4. * Egamma * _beamLorentzGamma)) * (W * W / (4. * Egamma * _beamLorentzGamma));
   
 	if ((_bbs.beam1().A() == 1) && (_bbs.beam2().A() == 1)){
 	   // proton-proton, no scaling needed
 	   csgA = sigmagp(Wgp);
 	} else {
 	   // coherent AA interactions
            int A_1 = _bbs.beam1().A(); 
            int A_2 = _bbs.beam2().A(); 
 
 	   // Calculate V.M.+proton cross section
            // cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha); 
            cs = sigma_N(Wgp); //Use member function instead 
     
 	   // Calculate V.M.+nucleus cross section
 	   cvma = sigma_A(cs,beam); 
 
            // Do impulse approximation here
            if( _impulseSelected == 1){
              if( beam == 1 ){
 	       cvma = A_1*cs;
 	     } else if ( beam == 2 ){
                cvma = A_2*cs;
 	     }   
            }	   
 
 	   // Calculate Av = dsigma/dt(t=0) Note Units: fm**s/Gev**2
 	   Av = (alpha * cvma * cvma) / (16. * pi * _vmPhotonCoupling * hbarc * hbarc);
 
 	   tmax   = tmin + 0.25;
 	   ax     = 0.5 * (tmax - tmin);
 	   bx     = 0.5 * (tmax + tmin);
 	   csgA   = 0.;
 	   for (int k = 1; k < NGAUSS; ++k) { 
 
 	       t    = ax * xg[k] + bx;
                if( A_1 == 1 && A_2 != 1){ 
 		  csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
                }else if(A_2 ==1 && A_1 != 1){
 		  csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
                }else{     
                   if( beam==1 ){
  		     csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
                   }else if(beam==2){
  		     csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);	
   		  }else{
 		     cout<<"Something went wrong here, beam= "<<beam<<endl; 
                   }
                }
 
 	       t    = ax * (-xg[k]) + bx;
                if( A_1 == 1 && A_2 != 1){ 
 			  csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);
                }else if(A_2 ==1 && A_1 != 1){
 			  csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
                }else{     
                   if( beam==1 ){
  			    csgA = csgA + ag[k] * _bbs.beam1().formFactor(t) * _bbs.beam1().formFactor(t);
                   }else if(beam==2){
  			    csgA = csgA + ag[k] * _bbs.beam2().formFactor(t) * _bbs.beam2().formFactor(t);	
   		  }else{
 			    cout<<"Something went wrong here, beam= "<<beam<<endl; 
                   }
 	       }
 	   }
 	   csgA = 0.5 * (tmax - tmin) * csgA;
 	   csgA = Av * csgA;
 	}
 	return csgA;	
 }
 
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::photonFlux(const double Egamma, const int beam)
 {
 	// This routine gives the photon flux as a function of energy Egamma
 	// It works for arbitrary nuclei and gamma; the first time it is
 	// called, it calculates a lookup table which is used on
 	// subsequent calls.
 	// It returns dN_gamma/dE (dimensions 1/E), not dI/dE
 	// energies are in GeV, in the lab frame
 	// rewritten 4/25/2001 by SRK
 
         // NOTE: beam (=1,2) defines the photon TARGET
 
 	double lEgamma,Emin,Emax;
 	static double lnEmax, lnEmin, dlnE;
 	double stepmult,energy,rZ;
 	int nbstep,nrstep,nphistep,nstep;
 	double bmin,bmax,bmult,biter,bold,integratedflux;
 	double fluxelement,deltar,riter;
 	double deltaphi,phiiter,dist;
 	static double dide[401];
 	double lnElt;
 	double flux_r; 
 	double Xvar;
 	int Ilt;
 	double RNuc=0.,RSum=0.;
 
 	RSum=_bbs.beam1().nuclearRadius()+_bbs.beam2().nuclearRadius();
         if( beam == 1){
           rZ=double(_bbs.beam2().Z());
           RNuc = _bbs.beam1().nuclearRadius();
         } else { 
 	  rZ=double(_bbs.beam1().Z());
           RNuc = _bbs.beam2().nuclearRadius();
         }
 
 	static int  Icheck = 0;
 	static int  Ibeam  = 0; 
   
 	//Check first to see if pp 
 	if( _bbs.beam1().A()==1 && _bbs.beam2().A()==1 ){
 		int nbsteps = 400;
 		double bmin = 0.5;
 		double bmax = 5.0 + (5.0*_beamLorentzGamma*hbarc/Egamma);
 		double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
  
 		double local_sum=0.0;
 
 		// Impact parameter loop 
 		for (int i = 0; i<=nbsteps;i++){
 
 			double bnn0 = bmin*exp(i*dlnb);
 			double bnn1 = bmin*exp((i+1)*dlnb);
 			double db   = bnn1-bnn0;
         
 			double ppslope = 19.0; 
 			double GammaProfile = exp(-bnn0*bnn0/(2.*hbarc*hbarc*ppslope));  
 			double PofB0 = 1. - (1. - GammaProfile)*(1. - GammaProfile);   
 			GammaProfile = exp(-bnn1*bnn1/(2.*hbarc*hbarc*ppslope));  
 			double PofB1 = 1. - (1. - GammaProfile)*(1. - GammaProfile);   
 
                         double loc_nofe0 = _bbs.beam1().photonDensity(bnn0,Egamma);
 			double loc_nofe1 = _bbs.beam2().photonDensity(bnn1,Egamma);
 
 			local_sum += 0.5*loc_nofe0*(1. - PofB0)*2.*starlightConstants::pi*bnn0*db; 
 			local_sum += 0.5*loc_nofe1*(1. - PofB1)*2.*starlightConstants::pi*bnn1*db; 
 
 		}
 		// End Impact parameter loop 
 		return local_sum;
 	}
 
 	//   first call or new beam?  - initialize - calculate photon flux
 	Icheck=Icheck+1;
 
 	// Do the numerical integration only once for symmetric systems. 
         if( Icheck > 1 && _bbs.beam1().A() == _bbs.beam2().A() && _bbs.beam1().Z() == _bbs.beam2().Z() ) goto L1000f;
         // For asymmetric systems check if we have another beam 
 	if( Icheck > 1 && beam == Ibeam ) goto L1000f; 
         Ibeam = beam; 
   
   	//  Nuclear breakup is done by PofB
 	//  collect number of integration steps here, in one place
   
 	nbstep=1200;
 	nrstep=60;
 	nphistep=40;
   
 	//  this last one is the number of energy steps
 	nstep=100;
  
 	// following previous choices, take Emin=10 keV at LHC, Emin = 1 MeV at RHIC
   	Emin=1.E-5;
 	if (_beamLorentzGamma < 500) 
 		Emin=1.E-3;
   
         Emax=_maxPhotonEnergy; 
 	// Emax=12.*hbarc*_beamLorentzGamma/RSum;
  
 	//     >> lnEmin <-> ln(Egamma) for the 0th bin
 	//     >> lnEmax <-> ln(Egamma) for the last bin
   	lnEmin=log(Emin);
 	lnEmax=log(Emax);
 	dlnE=(lnEmax-lnEmin)/nstep; 
 
         printf("Calculating photon flux from Emin = %e GeV to Emax = %e GeV (CM frame) for source with Z = %3.0f \n", Emin, Emax, rZ);
 	
 	stepmult= exp(log(Emax/Emin)/double(nstep));
 	energy=Emin;
   
 	for (int j = 1; j<=nstep;j++){
 		energy=energy*stepmult;
     
 		//  integrate flux over 2R_A < b < 2R_A+ 6* gamma hbar/energy
 		//  use exponential steps
     
 		bmin=0.8*RSum; //Start slightly below 2*Radius 
 		bmax=bmin + 6.*hbarc*_beamLorentzGamma/energy;
     
 		bmult=exp(log(bmax/bmin)/double(nbstep));
 		biter=bmin;
 		integratedflux=0.;
 
 		if( (_bbs.beam1().A() == 1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A() == 1 && _bbs.beam1().A() != 1) ){
 		    // This is pA 
 
 		  if( _productionMode == PHOTONPOMERONINCOHERENT ){
 		      // This pA incoherent, proton is the target
 
 		      int nbsteps = 400;
 		      double bmin = 0.7*RSum;
 		      double bmax = 2.0*RSum + (8.0*_beamLorentzGamma*hbarc/energy);
 		      double dlnb = (log(bmax)-log(bmin))/(1.*nbsteps);
 
 		      double local_sum=0.0;
 
 		      // Impact parameter loop 
 		      for (int i = 0; i<=nbsteps; i++){
 
         		  double bnn0 = bmin*exp(i*dlnb);
 			  double bnn1 = bmin*exp((i+1)*dlnb);
 			  double db   = bnn1-bnn0;
 
                           double PofB0 = _bbs.probabilityOfBreakup(bnn0); 
                           double PofB1 = _bbs.probabilityOfBreakup(bnn1); 
       
 			  double loc_nofe0 = 0.0;
 			  double loc_nofe1 = 0.0; 
                           if( _bbs.beam1().A() == 1 ){
 			    loc_nofe0 = _bbs.beam2().photonDensity(bnn0,energy);
 			    loc_nofe1 = _bbs.beam2().photonDensity(bnn1,energy);
                           }
 		          else if( _bbs.beam2().A() == 1 ){
 			    loc_nofe0 = _bbs.beam1().photonDensity(bnn0,energy);
 			    loc_nofe1 = _bbs.beam1().photonDensity(bnn1,energy);			    
                           }
 
                           // cout<<" i: "<<i<<" bnn0: "<<bnn0<<" PofB0: "<<PofB0<<" loc_nofe0: "<<loc_nofe0<<endl; 
 
 			  local_sum += 0.5*loc_nofe0*PofB0*2.*starlightConstants::pi*bnn0*db; 
 			  local_sum += 0.5*loc_nofe1*PofB1*2.*starlightConstants::pi*bnn1*db;  
 		      }  // End Impact parameter loop 
 	              integratedflux = local_sum; 
                     } else if ( _productionMode == PHOTONPOMERONNARROW ||  _productionMode == PHOTONPOMERONWIDE ){
                       // This is pA coherent, nucleus is the target 
                       double localbmin = 0.0;   
                       if( _bbs.beam1().A() == 1 ){
 			localbmin = _bbs.beam2().nuclearRadius() + 0.7; 
 		      }
                       if( _bbs.beam2().A() == 1 ){ 
 			localbmin = _bbs.beam1().nuclearRadius() + 0.7; 
                       }
                       integratedflux = nepoint(energy,localbmin); 
 		    }
 		}else{ 
 		// This is AA
 		for (int jb = 1; jb<=nbstep;jb++){
 		    bold=biter;
 		    biter=biter*bmult;
 		    // When we get to b>20R_A change methods - just take the photon flux
 		    //  at the center of the nucleus.
 		    if (biter > (10.*RNuc)){
 		       // if there is no nuclear breakup or only hadronic breakup, which only
 		       // occurs at smaller b, we can analytically integrate the flux from b~20R_A
 		       // to infinity, following Jackson (2nd edition), Eq. 15.54
 		       Xvar=energy*biter/(hbarc*_beamLorentzGamma);
 		       // Here, there is nuclear breakup.  So, we can't use the integrated flux
 		       // However, we can do a single flux calculation, at the center of the nucleus
 		       // Eq. 41 of Vidovic, Greiner and Soff, Phys.Rev.C47,2308(1993), among other places
 		       // this is the flux per unit area
 		       fluxelement  = (rZ*rZ*alpha*energy)*
 				       (bessel::dbesk1(Xvar))*(bessel::dbesk1(Xvar))/
 				       ((pi*_beamLorentzGamma*hbarc)*
 				       (pi*_beamLorentzGamma*hbarc));
 	    
                    } else {
 		       // integrate over nuclear surface. n.b. this assumes total shadowing -
 		       // treat photons hitting the nucleus the same no matter where they strike
 		       fluxelement=0.;
 		       deltar=RNuc/double(nrstep);
 		       riter=-deltar/2.;
           
 		       for (int jr =1; jr<=nrstep;jr++){
 			   riter=riter+deltar;
 			   // use symmetry;  only integrate from 0 to pi (half circle)
 			   deltaphi=pi/double(nphistep);
 			   phiiter=0.;
             
 			   for( int jphi=1;jphi<= nphistep;jphi++){
 			       phiiter=(double(jphi)-0.5)*deltaphi;
 			       // dist is the distance from the center of the emitting nucleus 
 			       // to the point in question
 			       dist=sqrt((biter+riter*cos(phiiter))*(biter+riter*
 				     cos(phiiter))+(riter*sin(phiiter))*(riter*sin(phiiter)));
 			       Xvar=energy*dist/(hbarc*_beamLorentzGamma);  
 			       flux_r = (rZ*rZ*alpha*energy)*
 				     (bessel::dbesk1(Xvar)*bessel::dbesk1(Xvar))/
 				     ((pi*_beamLorentzGamma*hbarc)*
 				     (pi*_beamLorentzGamma*hbarc));
 	      
 				     //  The surface  element is 2.* delta phi* r * delta r
 				     //  The '2' is because the phi integral only goes from 0 to pi
 				     fluxelement=fluxelement+flux_r*2.*deltaphi*riter*deltar;
 				     //  end phi and r integrations
 			   }//for(jphi)
 		       }//for(jr)
 			  //  average fluxelement over the nuclear surface
 			  fluxelement=fluxelement/(pi*RNuc*RNuc);
 		   }//else
 			  //  multiply by volume element to get total flux in the volume element
 			  fluxelement=fluxelement*2.*pi*biter*(biter-bold);
 			  //  modulate by the probability of nuclear breakup as f(biter)
                           // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl; 
 			  if (_beamBreakupMode > 1){
 			      fluxelement=fluxelement*_bbs.probabilityOfBreakup(biter);
 			  }
                           // cout<<" jb: "<<jb<<" biter: "<<biter<<" fluxelement: "<<fluxelement<<endl; 
 			  integratedflux=integratedflux+fluxelement;
       
 		} //end loop over impact parameter 
 	    }  //end of else (pp, pA, AA) 
     
 	    //  In lookup table, store k*dN/dk because it changes less
 	    //  so the interpolation should be better    
 	    dide[j]=integratedflux*energy;                                     
 	}//end loop over photon energy 
        
 	//  for 2nd and subsequent calls, use lookup table immediately
   
  L1000f:
   
 	lEgamma=log(Egamma);
 	if (lEgamma < (lnEmin+dlnE) ||  lEgamma  > lnEmax){
 		flux_r=0.0;
 		// cout<<"  WARNING: Egamma outside defined range. Egamma= "<<Egamma
 		//    <<"   "<<lnEmax<<" "<<(lnEmin+dlnE)<<endl;
 	}
 	else{
 		//       >> Egamma between Ilt and Ilt+1
 		Ilt = int((lEgamma-lnEmin)/dlnE);
 		//       >> ln(Egamma) for first point 
 		lnElt = lnEmin + Ilt*dlnE; 
 		//       >> Interpolate
 		flux_r = dide[Ilt] + ((lEgamma-lnElt)/dlnE)*(dide[Ilt+1]- dide[Ilt]);
 		flux_r = flux_r/Egamma;
 	}
   
 	return flux_r;
 }
 
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::nepoint(const double Egamma,
                                    const double bmin)
 {
 	// Function for the spectrum of virtual photons,
 	// dn/dEgamma, for a point charge q=Ze sweeping
 	// past the origin with velocity gamma
 	// (=1/SQRT(1-(V/c)**2)) integrated over impact
 	// parameter from bmin to infinity
 	// See Jackson eq15.54 Classical Electrodynamics
 	// Declare Local Variables
 	double beta,X,C1,bracket,nepoint_r;
   
 	beta = sqrt(1.-(1./(_beamLorentzGamma*_beamLorentzGamma)));
 	X = (bmin*Egamma)/(beta*_beamLorentzGamma*hbarc);
   
 	bracket = -0.5*beta*beta*X*X*(bessel::dbesk1(X)*bessel::dbesk1(X)
 	                              -bessel::dbesk0(X)*bessel::dbesk0(X));
 
 	bracket = bracket+X*bessel::dbesk0(X)*bessel::dbesk1(X);
   
 	// Note: NO  Z*Z!!
 	C1=(2.*alpha)/pi;
   
 	nepoint_r = C1*(1./beta)*(1./beta)*(1./Egamma)*bracket;
   
 	return nepoint_r;
   
 }
 
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::sigmagp(const double Wgp)
 {
 	// Function for the gamma-proton --> VectorMeson
 	// cross section. Wgp is the gamma-proton CM energy.
 	// Unit for cross section: fm**2
   
 	double sigmagp_r=0.;
   
 	switch(_particleType)
 		{ 
 		case RHO:
 		case RHOZEUS:
 		case FOURPRONG:
 			sigmagp_r=1.E-4*(5.0*exp(0.22*log(Wgp))+26.0*exp(-1.23*log(Wgp)));
 			break;
 		case OMEGA:
 			sigmagp_r=1.E-4*(0.55*exp(0.22*log(Wgp))+18.0*exp(-1.92*log(Wgp)));
 			break;                                                      
 		case PHI:
 			sigmagp_r=1.E-4*0.34*exp(0.22*log(Wgp));
 			break;
 		case JPSI:
 		case JPSI_ee:
 		case JPSI_mumu:
 		case JPSI_ppbar:
 			sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
 			sigmagp_r*=sigmagp_r;
 			sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
 			// sigmagp_r=1.E-4*0.0015*exp(0.80*log(Wgp));
 			break;
 		case JPSI2S:
 		case JPSI2S_ee:
 		case JPSI2S_mumu:
 			sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
 			sigmagp_r*=sigmagp_r;
 			sigmagp_r*=1.E-4*0.00406*exp(0.65*log(Wgp));
 			sigmagp_r*=0.166;  
 			//      sigmagp_r=0.166*(1.E-4*0.0015*exp(0.80*log(Wgp)));
 			break;
 		case UPSILON:
 		case UPSILON_ee:
 		case UPSILON_mumu:
 			//       >> This is W**1.7 dependence from QCD calculations
 			//  sigmagp_r=1.E-10*(0.060)*exp(1.70*log(Wgp));
 			sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
 			sigmagp_r*=sigmagp_r;
 			sigmagp_r*=1.E-10*6.4*exp(0.74*log(Wgp));
 			break;
 		case UPSILON2S:
 		case UPSILON2S_ee:
 		case UPSILON2S_mumu:
 		        // sigmagp_r=1.E-10*(0.0259)*exp(1.70*log(Wgp));
 		        sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
 			sigmagp_r*=sigmagp_r;
 			sigmagp_r*=1.E-10*2.9*exp(0.74*log(Wgp)); 
 			break;
 		case UPSILON3S:
 		case UPSILON3S_ee:
 		case UPSILON3S_mumu:
 		        // sigmagp_r=1.E-10*(0.0181)*exp(1.70*log(Wgp));
 		        sigmagp_r=(1.0-((_channelMass+_ip->protonMass())*(_channelMass+_ip->protonMass()))/(Wgp*Wgp));
 			sigmagp_r*=sigmagp_r;
 			sigmagp_r*=1.E-10*2.1*exp(0.74*log(Wgp)); 
 			break;
 		default: cout<< "!!!  ERROR: Unidentified Vector Meson: "<< _particleType <<endl;
 		}                                                                  
 	return sigmagp_r;
 }
 
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::sigma_A(const double sig_N, const int beam)
 {                                                         
 	// Nuclear Cross Section
 	// sig_N,sigma_A in (fm**2) 
 
 	double sum;
 	double b,bmax,Pint,arg,sigma_A_r;
   
 	int NGAUSS;
   
 	double xg[17]=
 		{.0,
 		 .0483076656877383162,.144471961582796493,
 		 .239287362252137075, .331868602282127650,
 		 .421351276130635345, .506899908932229390,
 		 .587715757240762329, .663044266930215201,
 		 .732182118740289680, .794483795967942407,
 		 .849367613732569970, .896321155766052124,
 		 .934906075937739689, .964762255587506430,
 		 .985611511545268335, .997263861849481564
 		};
   
 	double ag[17]=
 		{.0,
 		 .0965400885147278006, .0956387200792748594,
 		 .0938443990808045654, .0911738786957638847,
 		 .0876520930044038111, .0833119242269467552,
 		 .0781938957870703065, .0723457941088485062,
 		 .0658222227763618468, .0586840934785355471,
 		 .0509980592623761762, .0428358980222266807,
 		 .0342738629130214331, .0253920653092620595,
 		 .0162743947309056706, .00701861000947009660
 		};
   
 	NGAUSS=16;
  
 	// Check if one or both beams are nuclei 
         int A_1 = _bbs.beam1().A(); 
         int A_2 = _bbs.beam2().A(); 
         if( A_1 == 1 && A_2 == 1)cout<<" This is pp, you should not be here..."<<endl;  
 
 	// CALCULATE P(int) FOR b=0.0 - bmax (fm)
 	bmax = 25.0;
 	sum  = 0.;
 	for(int IB=1;IB<=NGAUSS;IB++){
     
 		b = 0.5*bmax*xg[IB]+0.5*bmax;
 
                 if( A_1 == 1 && A_2 != 1){  
                   arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
                 }else if(A_2 == 1 && A_1 != 1){
                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
                 }else{
 		  // Check which beam is target 
                   if( beam == 1 ){ 
                     arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
                   }else if( beam==2 ){
                     arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
                   }else{
                     cout<<" Something went wrong here, beam= "<<beam<<endl; 
                   } 
                 }
     
 		Pint=1.0-exp(arg);
 		// If this is a quantum Glauber calculation, use the quantum Glauber formula
 		if (_quantumGlauber == 1){Pint=2.0*(1.0-exp(arg/2.0));}
 						    
 		sum=sum+2.*pi*b*Pint*ag[IB];
     
 		b = 0.5*bmax*(-xg[IB])+0.5*bmax;
 
                 if( A_1 == 1 && A_2 != 1){  
                   arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
                 }else if(A_2 == 1 && A_1 != 1){
                   arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
                 }else{ 
 		  // Check which beam is target 
                   if( beam == 1 ){ 
                     arg=-sig_N*_bbs.beam1().rho0()*_bbs.beam1().thickness(b);
                   }else if(beam==2){
                     arg=-sig_N*_bbs.beam2().rho0()*_bbs.beam2().thickness(b);
                   }else{
                     cout<<" Something went wrong here, beam= "<<beam<<endl; 
                   } 
                 }
 
 		Pint=1.0-exp(arg);
 		// If this is a quantum Glauber calculation, use the quantum Glauber formula
 		if (_quantumGlauber == 1){Pint=2.0*(1.0-exp(arg/2.0));}
 		sum=sum+2.*pi*b*Pint*ag[IB];
 
 	}
 
 	sum=0.5*bmax*sum;
   
 	sigma_A_r=sum;
  
 	return sigma_A_r;
 }
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::sigma_N(const double Wgp)
 {                                                         
         // Nucleon Cross Section in (fm**2) 
         double cs = sqrt(16. * pi * _vmPhotonCoupling * _slopeParameter * hbarc * hbarc * sigmagp(Wgp) / alpha);
         return cs;
 }
 
 
 //______________________________________________________________________________
 double
 photonNucleusCrossSection::breitWigner(const double W,
                                        const double C)
 {
-	// use simple fixed-width s-wave Breit-Wigner without coherent backgorund for rho'
+	// use simple fixed-width s-wave Breit-Wigner without coherent background for rho'
 	// (PDG '08 eq. 38.56)
 	if(_particleType==FOURPRONG) {
 		if (W < 4.01 * _ip->pionChargedMass())
 			return 0;
 		const double termA  = _channelMass * _width;
 		const double termA2 = termA * termA;
 		const double termB  = W * W - _channelMass * _channelMass;
 		return C * _ANORM * _ANORM * termA2 / (termB * termB + termA2);
 	}
 
 	// Relativistic Breit-Wigner according to J.D. Jackson,
 	// Nuovo Cimento 34, 6692 (1964), with nonresonant term. A is the strength
 	// of the resonant term and b the strength of the non-resonant
 	// term. C is an overall normalization.
 
 	double ppi=0.,ppi0=0.,GammaPrim,rat;
 	double aa,bb,cc;
   
 	double nrbw_r;
 
 	// width depends on energy - Jackson Eq. A.2
 	// if below threshold, then return 0.  Added 5/3/2001 SRK
 	// 0.5% extra added for safety margin
         // omega added here 10/26/2014 SRK
 	if( _particleType==RHO ||_particleType==RHOZEUS || _particleType==OMEGA){  
 		if (W < 2.01*_ip->pionChargedMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt( ((W/2.)*(W/2.)) - _ip->pionChargedMass() * _ip->pionChargedMass());
 		ppi0=0.358;
 	}
   
 	// handle phi-->K+K- properly
 	if (_particleType  ==  PHI){
 		if (W < 2.*_ip->kaonChargedMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt( ((W/2.)*(W/2.))- _ip->kaonChargedMass()*_ip->kaonChargedMass());
 		ppi0=sqrt( ((_channelMass/2.)*(_channelMass/2.))-_ip->kaonChargedMass()*_ip->kaonChargedMass());
 	}
 
 	//handle J/Psi-->e+e- properly
 	if (_particleType==JPSI || _particleType==JPSI2S){
 		if(W<2.*_ip->mel()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
 	}
 	if (_particleType==JPSI_ee){
 		if(W<2.*_ip->mel()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());   
 	}
 	if (_particleType==JPSI_mumu){
 		if(W<2.*_ip->muonMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
 	}
 	if (_particleType==JPSI_ppbar){
 		if(W<2.*_ip->protonMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->protonMass()*_ip->protonMass());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->protonMass()*_ip->protonMass());
 	}
 	if (_particleType==JPSI2S_ee){
 		if(W<2.*_ip->mel()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());   
 	}
 	if (_particleType==JPSI2S_mumu){
 		if(W<2.*_ip->muonMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
 	}
 
 	if(_particleType==UPSILON || _particleType==UPSILON2S ||_particleType==UPSILON3S ){ 
 		if (W<2.*_ip->muonMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
 	}
   
 	if(_particleType==UPSILON_mumu || _particleType==UPSILON2S_mumu ||_particleType==UPSILON3S_mumu ){ 
 		if (W<2.*_ip->muonMass()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->muonMass()*_ip->muonMass());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->muonMass()*_ip->muonMass());
 	}
   
 	if(_particleType==UPSILON_ee || _particleType==UPSILON2S_ee ||_particleType==UPSILON3S_ee ){ 
 		if (W<2.*_ip->mel()){
 			nrbw_r=0.;
 			return nrbw_r;
 		}
 		ppi=sqrt(((W/2.)*(W/2.))-_ip->mel()*_ip->mel());
 		ppi0=sqrt(((_channelMass/2.)*(_channelMass/2.))-_ip->mel()*_ip->mel());
 	}
   
 	if(ppi==0.&&ppi0==0.) 
 		cout<<"Improper Gammaacrosssection::breitwigner, ppi&ppi0=0."<<endl;
   
 	rat=ppi/ppi0;
 	GammaPrim=_width*(_channelMass/W)*rat*rat*rat;
   
 	aa=_ANORM*sqrt(GammaPrim*_channelMass*W);
 	bb=W*W-_channelMass*_channelMass;
 	cc=_channelMass*GammaPrim;
   
 	// First real part squared 
 	nrbw_r = (( (aa*bb)/(bb*bb+cc*cc) + _BNORM)*( (aa*bb)/(bb*bb+cc*cc) + _BNORM));
   
 	// Then imaginary part squared 
 	nrbw_r = nrbw_r + (( (aa*cc)/(bb*bb+cc*cc) )*( (aa*cc)/(bb*bb+cc*cc) ));
 
 	//  Alternative, a simple, no-background BW, following J. Breitweg et al.
 	//  Eq. 15 of Eur. Phys. J. C2, 247 (1998).  SRK 11/10/2000
 	//      nrbw_r = (_ANORM*_mass*GammaPrim/(bb*bb+cc*cc))**2
 
   
 	nrbw_r = C*nrbw_r;
   
 	return nrbw_r;    
 }
Index: trunk/src/inputParameters.cpp
===================================================================
--- trunk/src/inputParameters.cpp	(revision 305)
+++ trunk/src/inputParameters.cpp	(revision 306)
@@ -1,841 +1,842 @@
 /////////////////////////////////////////////////////////////////////////// 
 // 
 // Copyright 2010 
 // 
 // This file is part of starlight. 
 // 
 // starlight is free software: you can redistribute it and/or modify 
 // it under the terms of the GNU General Public License as published by 
 // the Free Software Foundation, either version 3 of the License, or 
 // (at your option) any later version. 
 // 
 // starlight is distributed in the hope that it will be useful, 
 // but WITHOUT ANY WARRANTY; without even the implied warranty of 
 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the 
 // GNU General Public License for more details. 
 // 
 // You should have received a copy of the GNU General Public License 
 // along with starlight. If not, see <http://www.gnu.org/licenses/>. 
 // 
 /////////////////////////////////////////////////////////////////////////// 
 // 
 // File and Version Information: 
 // $Rev::     $: revision of last commit // $Author::         $: author of last commit 
 // $Date::                            $: date of last commit // 
 // Description: 
 // 
 ///////////////////////////////////////////////////////////////////////////
 
 
 #include <iostream>
 #include <fstream>
 
 #include "reportingUtils.h"
 #include "starlightconstants.h"
 #include "inputParameters.h"
 #include "inputParser.h"
 #include "starlightconfig.h"
 #include <cmath>
 #include <cstring>
 #include "randomgenerator.h"
 
 using namespace std;
 using namespace starlightConstants;
 
 parameterlist parameterbase::_parameters;
 
 #define REQUIRED true
 #define NOT_REQUIRED false
 
 //______________________________________________________________________________
 inputParameters::inputParameters()
         : _baseFileName          ("baseFileName","slight"),
  	  _beam1Z                ("BEAM_1_Z",0),
 	  _beam1A                ("BEAM_1_A",0),
 	  _beam2Z                ("BEAM_2_Z",0),
 	  _beam2A                ("BEAM_2_A",0),
 	  _beam1LorentzGamma     ("BEAM_1_GAMMA",0),
 	  _beam2LorentzGamma     ("BEAM_2_GAMMA",0),
 	  _maxW                  ("W_MAX",0),
 	  _minW                  ("W_MIN",0),
 	  _nmbWBins              ("W_N_BINS",0),
 	  _maxRapidity           ("RAP_MAX",0),
 	  _nmbRapidityBins       ("RAP_N_BINS",0),
 	  _ptCutEnabled          ("CUT_PT",false, NOT_REQUIRED),
 	  _ptCutMin              ("PT_MIN",0, NOT_REQUIRED),
 	  _ptCutMax              ("PT_MAX",0, NOT_REQUIRED),
 	  _etaCutEnabled         ("CUT_ETA",false, NOT_REQUIRED),
 	  _etaCutMin             ("ETA_MIN",0, NOT_REQUIRED),
 	  _etaCutMax             ("ETA_MAX",0, NOT_REQUIRED),
 	  _productionMode        ("PROD_MODE",0),
 	  _nmbEventsTot          ("N_EVENTS",0),
 	  _prodParticleId        ("PROD_PID",0),
 	  _randomSeed            ("RND_SEED",0),
 	  _beamBreakupMode       ("BREAKUP_MODE",0),
 	  _interferenceEnabled   ("INTERFERENCE",false),
 	  _interferenceStrength  ("IF_STRENGTH",0),
 	  _maxPtInterference     ("INT_PT_MAX",0),
 	  _nmbPtBinsInterference ("INT_PT_N_BINS",0),
 	  _ptBinWidthInterference("INT_PT_WIDTH",0),
 	  _protonEnergy          ("PROTON_ENERGY",0, NOT_REQUIRED),
 	  _minGammaEnergy	 ("MIN_GAMMA_ENERGY",6.0, NOT_REQUIRED),
 	  _maxGammaEnergy	 ("MAX_GAMMA_ENERGY",600000.0, NOT_REQUIRED),
 	  _pythiaParams          ("PYTHIA_PARAMS","", NOT_REQUIRED),
 	  _pythiaFullEventRecord ("PYTHIA_FULL_EVENTRECORD",false, NOT_REQUIRED),
 	  _xsecCalcMethod	 ("XSEC_METHOD",0, NOT_REQUIRED),
           _axionMass             ("AXION_MASS",50, NOT_REQUIRED),  // AXION HACK
           _bslopeDefinition      ("BSLOPE_DEFINITION",0, NOT_REQUIRED),
 	  _bslopeValue           ("BSLOPE_VALUE",4.0,NOT_REQUIRED),
 	  _printVM               ("PRINT_VM",0,NOT_REQUIRED),
 	  _impulseVM             ("SELECT_IMPULSE_VM",0,NOT_REQUIRED),
 	  _quantumGlauber        ("QUANTUM_GLAUBER",0,NOT_REQUIRED),
 	  _bmin                  ("BMIN",0,NOT_REQUIRED),
           _bmax                  ("BMAX",0,NOT_REQUIRED),   
 
           _deuteronSlopePar      ("deuteronSlopePar"      , 9.5           , NOT_REQUIRED),
           _protonMass            ("protonMass"            , 0.938272046   , NOT_REQUIRED),
           _pionChargedMass       ("pionChargedMass"       , 0.13957018    , NOT_REQUIRED),
           _pionNeutralMass       ("pionNeutralMass"       , 0.1349766     , NOT_REQUIRED),
           _kaonChargedMass       ("kaonChargedMass"       , 0.493677      , NOT_REQUIRED),
           _mel                   ("mel"                   , 0.000510998928, NOT_REQUIRED),
           _muonMass              ("muonMass"              , 0.1056583715  , NOT_REQUIRED),
           _tauMass               ("tauMass"               , 1.77682       , NOT_REQUIRED),
           _f0Mass                ("f0Mass"                , 0.990         , NOT_REQUIRED),
           _f0Width               ("f0Width"               , 0.100         , NOT_REQUIRED),
           _f0BrPiPi              ("f0BrPiPi"              , 1.0           , NOT_REQUIRED),
           _etaMass               ("etaMass"               , 0.547862      , NOT_REQUIRED),
           _etaWidth              ("etaWidth"              , 0.00000131    , NOT_REQUIRED),
           _etaPrimeMass          ("etaPrimeMass"          , 0.95778       , NOT_REQUIRED),
           _etaPrimeWidth         ("etaPrimeWidth"         , 0.000198      , NOT_REQUIRED),
           _etaCMass              ("etaCMass"              , 2.9836        , NOT_REQUIRED),
           _etaCWidth             ("etaCWidth"             , 0.0322        , NOT_REQUIRED),
           _f2Mass                ("f2Mass"                , 1.2751        , NOT_REQUIRED),
           _f2Width               ("f2Width"               , 0.1851        , NOT_REQUIRED),
           _f2BrPiPi              ("f2BrPiPi"              , 0.561         , NOT_REQUIRED),
           _a2Mass                ("a2Mass"                , 1.3183        , NOT_REQUIRED),
           _a2Width               ("a2Width"               , 0.105         , NOT_REQUIRED),
           _f2PrimeMass           ("f2PrimeMass"           , 1.525         , NOT_REQUIRED),
           _f2PrimeWidth          ("f2PrimeWidth"          , 0.073         , NOT_REQUIRED),
           _f2PrimeBrKK           ("f2PrimeBrKK"           , 0.887         , NOT_REQUIRED),
           _zoverz03Mass          ("zoverz03Mass"          , 1.540         , NOT_REQUIRED),
           _f0PartialggWidth      ("f0PartialggWidth"      , 0.29E-6       , NOT_REQUIRED),
           _etaPartialggWidth     ("etaPartialggWidth"     , 0.516E-6      , NOT_REQUIRED),
           _etaPrimePartialggWidth("etaPrimePartialggWidth", 4.35E-6       , NOT_REQUIRED),
           _etaCPartialggWidth    ("etaCPartialggWidth"    , 5.0E-6        , NOT_REQUIRED),
           _f2PartialggWidth      ("f2PartialggWidth"      , 3.03E-6       , NOT_REQUIRED),
           _a2PartialggWidth      ("a2PartialggWidth"      , 1.0E-6        , NOT_REQUIRED),
           _f2PrimePartialggWidth ("f2PrimePartialggWidth" , 0.081E-6      , NOT_REQUIRED),
           _zoverz03PartialggWidth("zoverz03PartialggWidth", 0.1E-6        , NOT_REQUIRED),
           _f0Spin                ("f0Spin"                , 0.0           , NOT_REQUIRED),
           _etaSpin               ("etaSpin"               , 0.0           , NOT_REQUIRED),
           _etaPrimeSpin          ("etaPrimeSpin"          , 0.0           , NOT_REQUIRED),
           _etaCSpin              ("etaCSpin"              , 0.0           , NOT_REQUIRED),
           _f2Spin                ("f2Spin"                , 2.0           , NOT_REQUIRED),
           _a2Spin                ("a2Spin"                , 2.0           , NOT_REQUIRED),
           _f2PrimeSpin           ("f2PrimeSpin"           , 2.0           , NOT_REQUIRED),
           _zoverz03Spin          ("zoverz03Spin"          , 2.0           , NOT_REQUIRED),
           _axionSpin             ("axionSpin"             , 0.0           , NOT_REQUIRED),
           _rho0Mass              ("rho0Mass"              , 0.769         , NOT_REQUIRED),
           _rho0Width             ("rho0Width"             , 0.1517        , NOT_REQUIRED),
           _rho0BrPiPi            ("rho0BrPiPi"            , 1.0           , NOT_REQUIRED),
           _rho0PrimeMass         ("rho0PrimeMass"         , 1.540         , NOT_REQUIRED),
           _rho0PrimeWidth        ("rho0PrimeWidth"        , 0.570         , NOT_REQUIRED),
           _rho0PrimeBrPiPi       ("rho0PrimeBrPiPi"       , 1.0           , NOT_REQUIRED),
           _OmegaMass             ("OmegaMass"             , 0.78265       , NOT_REQUIRED),
           _OmegaWidth            ("OmegaWidth"            , 0.00849       , NOT_REQUIRED),
           _OmegaBrPiPi           ("OmegaBrPiPi"           , 0.0153        , NOT_REQUIRED),
           _PhiMass               ("PhiMass"               , 1.019461      , NOT_REQUIRED),
           _PhiWidth              ("PhiWidth"              , 0.004266      , NOT_REQUIRED),
           _PhiBrKK               ("PhiBrKK"               , 0.489         , NOT_REQUIRED),
           _JpsiMass              ("JpsiMass"              , 3.096916      , NOT_REQUIRED),
           _JpsiWidth             ("JpsiWidth"             , 0.0000929     , NOT_REQUIRED),
           _JpsiBree              ("JpsiBree"              , 0.05971       , NOT_REQUIRED),
           _JpsiBrmumu            ("JpsiBrmumu"            , 0.05961       , NOT_REQUIRED),
           _JpsiBrppbar           ("JpsiBrppbar"           , 0.002120      , NOT_REQUIRED),
           _Psi2SMass             ("Psi2SMass"             , 3.686109      , NOT_REQUIRED),
           _Psi2SWidth            ("Psi2SWidth"            , 0.000299      , NOT_REQUIRED),
           _Psi2SBree             ("Psi2SBree"             , 0.00789       , NOT_REQUIRED),
           _Psi2SBrmumu           ("Psi2SBrmumu"           , 0.0079        , NOT_REQUIRED),
           _Upsilon1SMass         ("Upsilon1SMass"         , 9.46030       , NOT_REQUIRED),
           _Upsilon1SWidth        ("Upsilon1SWidth"        , 0.00005402    , NOT_REQUIRED),
           _Upsilon1SBree         ("Upsilon1SBree"         , 0.0238        , NOT_REQUIRED),
           _Upsilon1SBrmumu       ("Upsilon1SBrmumu"       , 0.0248        , NOT_REQUIRED),
           _Upsilon2SMass         ("Upsilon2SMass"         , 10.02326      , NOT_REQUIRED),
           _Upsilon2SWidth        ("Upsilon2SWidth"        , 0.00003198    , NOT_REQUIRED),
           _Upsilon2SBree         ("Upsilon2SBree"         , 0.0191        , NOT_REQUIRED),
           _Upsilon2SBrmumu       ("Upsilon2SBrmumu"       , 0.0193        , NOT_REQUIRED),
           _Upsilon3SMass         ("Upsilon3SMass"         , 10.3552       , NOT_REQUIRED),
           _Upsilon3SWidth        ("Upsilon3SWidth"        , 0.00002032    , NOT_REQUIRED),
           _Upsilon3SBree         ("Upsilon3SBree"         , 0.0218        , NOT_REQUIRED),
           _Upsilon3SBrmumu       ("Upsilon3SBrmumu"       , 0.0218        , NOT_REQUIRED)
 {
   // All parameters must be initialised in initialisation list! 
   // If not: error: 'parameter<T, validate>::parameter() [with T = unsigned int, bool validate = true]' is private
   // or similar
 	
         _ip.addParameter(_baseFileName);
 
 	_ip.addParameter(_beam1Z);
 	_ip.addParameter(_beam2Z);
 	_ip.addParameter(_beam1A);
 	_ip.addParameter(_beam2A);
 
 	_ip.addParameter(_beam1LorentzGamma);
 	_ip.addParameter(_beam2LorentzGamma);
 	
 	_ip.addParameter(_maxW);
 	_ip.addParameter(_minW);
 	
 	_ip.addParameter(_nmbWBins);
 
 	_ip.addParameter(_maxRapidity);
 	_ip.addParameter(_nmbRapidityBins);
 	
 	_ip.addParameter(_ptCutEnabled);
 	_ip.addParameter(_ptCutMin);
 	_ip.addParameter(_ptCutMax);
 	
 	_ip.addParameter(_etaCutEnabled);
 	_ip.addParameter(_etaCutMax);
 	_ip.addParameter(_etaCutMin);
 	
 	_ip.addParameter(_productionMode);
 	_ip.addParameter(_nmbEventsTot);
 	_ip.addParameter(_prodParticleId);
 	_ip.addParameter(_randomSeed);
 	_ip.addParameter(_beamBreakupMode);
 	_ip.addParameter(_interferenceEnabled);
 	_ip.addParameter(_interferenceStrength);
 	_ip.addParameter(_maxPtInterference);
 	_ip.addParameter(_nmbPtBinsInterference);
 	_ip.addParameter(_minGammaEnergy);
 	_ip.addParameter(_maxGammaEnergy);
 	_ip.addParameter(_pythiaParams);
 	_ip.addParameter(_pythiaFullEventRecord);
 	_ip.addParameter(_xsecCalcMethod);
         _ip.addParameter(_axionMass);     // AXION HACK
         _ip.addParameter(_bslopeDefinition); 
         _ip.addParameter(_bslopeValue); 
         _ip.addParameter(_printVM); 
         _ip.addParameter(_impulseVM);
 	_ip.addParameter(_quantumGlauber);
 	_ip.addParameter(_bmin);
 	_ip.addParameter(_bmax);
 
         _ip.addParameter(_deuteronSlopePar      );
         _ip.addParameter(_protonMass            );
         _ip.addParameter(_pionChargedMass       );
         _ip.addParameter(_pionNeutralMass       );
         _ip.addParameter(_kaonChargedMass       );
         _ip.addParameter(_mel                   );
         _ip.addParameter(_muonMass              );
         _ip.addParameter(_tauMass               );
         _ip.addParameter(_f0Mass                );
         _ip.addParameter(_f0Width               );
         _ip.addParameter(_f0BrPiPi              );
         _ip.addParameter(_etaMass               );
         _ip.addParameter(_etaWidth              );
         _ip.addParameter(_etaPrimeMass          );
         _ip.addParameter(_etaPrimeWidth         );
         _ip.addParameter(_etaCMass              );
         _ip.addParameter(_etaCWidth             );
         _ip.addParameter(_f2Mass                );
         _ip.addParameter(_f2Width               );
         _ip.addParameter(_f2BrPiPi              );
         _ip.addParameter(_a2Mass                );
         _ip.addParameter(_a2Width               );
         _ip.addParameter(_f2PrimeMass           );
         _ip.addParameter(_f2PrimeWidth          );
         _ip.addParameter(_f2PrimeBrKK           );
         _ip.addParameter(_zoverz03Mass          );
         _ip.addParameter(_f0PartialggWidth      );
         _ip.addParameter(_etaPartialggWidth     );
         _ip.addParameter(_etaPrimePartialggWidth);
         _ip.addParameter(_etaCPartialggWidth    );
         _ip.addParameter(_f2PartialggWidth      );
         _ip.addParameter(_a2PartialggWidth      );
         _ip.addParameter(_f2PrimePartialggWidth );
         _ip.addParameter(_zoverz03PartialggWidth);
         _ip.addParameter(_f0Spin                );
         _ip.addParameter(_etaSpin               );
         _ip.addParameter(_etaPrimeSpin          );
         _ip.addParameter(_etaCSpin              );
         _ip.addParameter(_f2Spin                );
         _ip.addParameter(_a2Spin                );
         _ip.addParameter(_f2PrimeSpin           );
         _ip.addParameter(_zoverz03Spin          );
         _ip.addParameter(_axionSpin             );
         _ip.addParameter(_rho0Mass              );
         _ip.addParameter(_rho0Width             );
         _ip.addParameter(_rho0BrPiPi            );
         _ip.addParameter(_rho0PrimeMass         );
         _ip.addParameter(_rho0PrimeWidth        );
         _ip.addParameter(_rho0PrimeBrPiPi       );
         _ip.addParameter(_OmegaMass             );
         _ip.addParameter(_OmegaWidth            );
         _ip.addParameter(_OmegaBrPiPi           );
         _ip.addParameter(_PhiMass               );
         _ip.addParameter(_PhiWidth              );
         _ip.addParameter(_PhiBrKK               );
         _ip.addParameter(_JpsiMass              );
         _ip.addParameter(_JpsiWidth             );
         _ip.addParameter(_JpsiBree              );
         _ip.addParameter(_JpsiBrmumu            );
         _ip.addParameter(_JpsiBrppbar           );
         _ip.addParameter(_Psi2SMass             );
         _ip.addParameter(_Psi2SWidth            );
         _ip.addParameter(_Psi2SBree             );
         _ip.addParameter(_Psi2SBrmumu           );
         _ip.addParameter(_Upsilon1SMass         );
         _ip.addParameter(_Upsilon1SWidth        );
         _ip.addParameter(_Upsilon1SBree         );
         _ip.addParameter(_Upsilon1SBrmumu       );
         _ip.addParameter(_Upsilon2SMass         );
         _ip.addParameter(_Upsilon2SWidth        );
         _ip.addParameter(_Upsilon2SBree         );
         _ip.addParameter(_Upsilon2SBrmumu       );
         _ip.addParameter(_Upsilon3SMass         );
         _ip.addParameter(_Upsilon3SWidth        );
         _ip.addParameter(_Upsilon3SBree         );
         _ip.addParameter(_Upsilon3SBrmumu       );
 }
 
 
 //______________________________________________________________________________
 inputParameters::~inputParameters()
 { }
 
 
 //______________________________________________________________________________
 bool
 inputParameters::configureFromFile(const std::string &_configFileName)
 {
 
 	int nParameters = _ip.parseFile(_configFileName);
 	
 	if(nParameters == -1) 
 	{
 		printWarn << "could not open file '" << _configFileName << "'" << endl;
 	  return false;
 	}
 	
 	
 	if(_ip.validateParameters(cerr))
 		printInfo << "successfully read input parameters from '" << _configFileName << "'" << endl;
 	else {
  		printWarn << "problems reading input parameters from '" << _configFileName << "'" << endl
  		          << *this;
  		return false;
  	}
 
 	/// temporary output test
 	printWarn<< * this<<endl;
 	return true;
 }
  bool inputParameters::init()
  {
 	
  	// Calculate beam gamma in CMS frame
  	double rap1 = acosh(beam1LorentzGamma());
 	double rap2 = -acosh(beam2LorentzGamma());
 	_beamLorentzGamma = cosh((rap1-rap2)/2);
 	
 	std::cout << "Rapidity beam 1: " << rap1 << ", rapidity beam 2: " << rap2 << ", rapidity CMS system: " << (rap1+rap2)/2 << ", beam gamma in CMS: " << _beamLorentzGamma<< std::endl;
 	_ptBinWidthInterference = maxPtInterference() / nmbPtBinsInterference();
 	_protonEnergy           = _beamLorentzGamma * protonMass();
 
 	// check for deuteron or tritium - these must be the second beam
 	if((beam1Z()==1) && (beam1A()==2)){
 		if((beam2Z()==1) && (beam2A()==2)){
 		printWarn << "deuteron-deuteron collisions are not supported" << endl;
 		return false;}
 		printWarn << "deuterium must be listed as the second nucleus" << endl;
 		return false;}
 
 	if( ((beam1Z()==1) && (beam1A()==3)) || ((beam2Z()==1) && (beam2A()==3)) ){
 		printWarn << "tritium is not currently supported" << endl;
 		return false;}
 	// check that rho production uses wide resonance option
 	if(_prodParticleId.value()==113 && _productionMode.value()==2){
 		printWarn << endl<< "For rho meson production, you should choose the wide resonance option (production mode = 3)" << endl;
 		return false;}
 
 	// define interaction type
 	switch (productionMode()) {
 	case 1:
 		_interactionType = PHOTONPHOTON;
 		break;
 	case 2:
 		_interactionType = PHOTONPOMERONNARROW;
 		break;
 	case 3:
 		_interactionType = PHOTONPOMERONWIDE;
 		break;
         case 4:
                 _interactionType = PHOTONPOMERONINCOHERENT;
                 break;
 	case 5:
 		_interactionType = PHOTONUCLEARSINGLE;
 		break;
 	case 6:
 		_interactionType = PHOTONUCLEARDOUBLE;
 		break;
 	case 7:
 		_interactionType = PHOTONUCLEARSINGLEPA;
 		break;
 	case 8:
 		_interactionType = PHOTONUCLEARSINGLEPAPY;
 		break;
 // 	case 9:
 // 		_interactionType = PHOTONPHOTONKNIEHL;
 // 		break;
 //         case 10:
 //                 _interactionType = PHOTONPHOTONKNIEHLMODIFIED;
 //                 break;
 	default:
 		printWarn << "unknown production mode '" << _productionMode << "'" << endl;
 		return false;
 	}
 	if( (_productionMode.value() ==4) && (_interferenceEnabled.value())) {
 		printWarn << " cannot enable interference for incoherent production " << endl;
 		return false;
 	}
   
 	double mass        = 0;
 	double width       = 0;
 	double defaultMinW = 0;  // default for _minW, unless it is defined later [GeV/c^2]
 	double defaultMaxW = 0;  // default for _maxW, unless it is defined later [GeV/c^2]
 	switch (prodParticleId()) {
 	case 11:  // e+e- pair
 		_particleType = ELECTRON;
 		_decayType    = LEPTONPAIR;
 		mass          = mel();
 		defaultMinW   = 2*mass; // default is 0.01; up to 0.15 is safe for Summer 2000 triggering for e+e- pairs
                 defaultMaxW     = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 13:  // mu+mu- pair
 		_particleType = MUON;
 		_decayType    = LEPTONPAIR;
 		defaultMinW   = 2 * muonMass();
                 defaultMaxW     = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 15:  // tau+tau- pair
 		_particleType = TAUON;
 		_decayType    = LEPTONPAIR;
 		defaultMinW   = 2 * tauMass();
                 defaultMaxW     = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 10015:  // tau+tau- pair
 		_particleType = TAUONDECAY;
 		_decayType    = LEPTONPAIR;
 		defaultMinW   = 2 * tauMass();
                 defaultMaxW   = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 
 // 	case 24:  // W+W- pair
 // 		_particleType = W;
 // 		_decayType    = WW;
 // 		defaultMinW   = 2 * muonMass();
 // 		break;
 	case 115:  // a_2(1320)
 		_particleType = A2;
 		_decayType    = SINGLEMESON;
 		mass          = a2Mass();
 		width         = a2Width();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW   = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 221:  // eta
 		_particleType = ETA;
 		_decayType    = SINGLEMESON;
 		mass          = etaMass();
 		width         = etaWidth();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW   = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 225:  // f_2(1270)
 		_particleType = F2;
 		_decayType    = SINGLEMESON;
 		mass          = f2Mass();
 		width         = f2Width();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW         = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = f2BrPiPi(); 
 		break;
 	case 331:  // eta'(958)
 		_particleType = ETAPRIME;
 		_decayType    = SINGLEMESON;
 		mass          = etaPrimeMass();
 		width         = etaPrimeWidth();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW         = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 335:  // f_2'(1525)
 		_particleType = F2PRIME;
 		_decayType    = SINGLEMESON;
 		mass          = f2PrimeMass();
 		width         = f2PrimeWidth();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW         = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = f2PrimeBrKK(); 
 		break;
 	case 441:  // eta_c(1s)
 		_particleType = ETAC;
 		_decayType    = SINGLEMESON;
 		mass          = etaCMass();
 		width         = etaCWidth();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW         = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
 	case 9010221:  // f_0(980), was orginally called 10221? updated to standard number
 		_particleType = F0;
 		_decayType    = SINGLEMESON;
 		mass          = f0Mass();
 		width         = f0Width();
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW         = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                _inputBranchingRatio = f0BrPiPi(); 
 		break;
 	case 33:  // Z"/Z0  This is the rho^0 rho^0 final state SRK
 		_particleType = ZOVERZ03;
 		_decayType    = SINGLEMESON;
                 defaultMinW   = 4*pionChargedMass();
                 defaultMaxW         = 1.6; // JES 6.17.2015 to avoid problems with no default
                 _inputBranchingRatio = 1.0; 
 		break;
         case 88:  // axion// AXION HACK, till break statement
                 _particleType = AXION;
                 _decayType    = SINGLEMESON;
                 mass          = _axionMass.value();
                 width         = 1/(64*starlightConstants::pi)*mass*mass*mass/(1000*1000);//Fix Lambda=1000 GeV, rescaling is trivial.
                 defaultMinW   = mass - 5*width; // JES 6.17.2015 to avoid problems with default of 0
                 defaultMaxW         = mass + 5*width; // JES 6.17.2015 to avoid problems with no default
                 break;    // AXION HACK, end
 // 	case 25: // Higgs
 // 		_particleType = HIGGS;
 // 		_decayType    = SINGLEMESON;
 // 		break;
 	case 113:  // rho(770)
 		_particleType = RHO;
 		_decayType    = WIDEVMDEFAULT;
 		mass          = rho0Mass();
 		width         = rho0Width();
 		defaultMinW   = 2 * pionChargedMass();
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = rho0BrPiPi(); 
 		break;
 	case 913:  // rho(770) with direct pi+pi- decay, interference given by ZEUS data
 		_particleType = RHOZEUS;
 		_decayType    = WIDEVMDEFAULT;
 		mass          = rho0Mass();
 		width         = rho0Width();
 		defaultMinW   = 2 * pionChargedMass();
 		defaultMaxW         = mass + 5 * width;  // use the same 1.5GeV max mass as ZEUS
 		_inputBranchingRatio = rho0BrPiPi(); 
 		break;
 	case 999:  // pi+pi-pi+pi- phase space decay
 		_particleType = FOURPRONG;
 		_decayType    = WIDEVMDEFAULT;
 		mass          = rho0PrimeMass();
 		width         = rho0PrimeWidth();		
 		defaultMinW   = 4 * pionChargedMass();
                 defaultMaxW     = sqrt(beam1LorentzGamma()*beam2LorentzGamma())*2*(starlightConstants::hbarc)/(1.2*pow(float(beam1A()),1./6.)*pow(float(beam2A()),1./6.)); // JES 6.17.2015 to avoid problems with no default
+		if (defaultMaxW>10.){defaultMaxW=10.;}  //SRK May 29, 2019, to avoid an unduly high defaultMaxW
 		_inputBranchingRatio = 1.0; 
 		break;
 	case 223:  // omega(782)
 		_particleType = OMEGA;
 		_decayType    = NARROWVMDEFAULT;  
 		mass          = OmegaMass();
 		width         = OmegaWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = OmegaBrPiPi(); 
 		break;
 	case 333:  // phi(1020)
 		_particleType = PHI;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = PhiMass();
 		width         = PhiWidth();
 		defaultMinW   = 2 * kaonChargedMass();
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = PhiBrKK(); 
 		break;
 	case 443:  // J/psi
 		_particleType = JPSI;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = JpsiMass();
 		width         = JpsiWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = (JpsiBree() + JpsiBrmumu())/2.; 
 		break;
    	case 443011:  // J/psi
 		_particleType = JPSI_ee;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = JpsiMass();
 		width         = JpsiWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = JpsiBree(); 
 		break;
 	case 443013:  // J/psi
 	        cout<<"In inputParameters setting J/psi mass!"<<endl; 
 		_particleType = JPSI_mumu;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = JpsiMass();
 		width         = JpsiWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = JpsiBrmumu(); 
 		break;
 	case 4432212:  // J/psi
 	        cout<<"In inputParameters setting J/psi mass!"<<endl; 
 		_particleType = JPSI_ppbar;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = JpsiMass();
 		width         = JpsiWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = JpsiBrppbar(); 
 		break;
 	case 444:  // psi(2S) 
 		_particleType = JPSI2S;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Psi2SMass();
 		width         = Psi2SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW         = mass + 5 * width;
 		_inputBranchingRatio = (Psi2SBree() + Psi2SBrmumu())/2.; 
 		break;
 	case 444011: // psi(2S)
 		_particleType = JPSI2S_ee;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Psi2SMass();
 		width         = Psi2SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Psi2SBree(); 
 		break;
 	case 444013:  // psi(2S)
 		_particleType = JPSI2S_mumu;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Psi2SMass();
 		width         = Psi2SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Psi2SBrmumu(); 
 		break;
 	case 553:  // Upsilon(1S) 
 		_particleType = UPSILON;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon1SMass();
 		width         = Upsilon1SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = (Upsilon1SBree() + Upsilon1SBrmumu())/2.; 
 		break;
 	case 553011:  // Upsilon
 		_particleType = UPSILON_ee;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon1SMass();
 		width         = Upsilon1SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Upsilon1SBree(); 
 		break;
 	case 553013:  // Upsilon
 		_particleType = UPSILON_mumu;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon1SMass();
 		width         = Upsilon1SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Upsilon1SBrmumu(); 
 		break;
 	case 554:  // Upsilon(2S)
 		_particleType = UPSILON2S;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon2SMass();
 		width         = Upsilon2SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = (Upsilon2SBree() + Upsilon2SBrmumu())/2.; 
 		break;
 	case 554011:  // Upsilon(2S)
 		_particleType = UPSILON2S_ee;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon2SMass();
 		width         = Upsilon2SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Upsilon2SBree(); 
 		break;
 	case 554013:  // Upsilon(2S)
 		_particleType = UPSILON2S_mumu;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon2SMass();
 		width         = Upsilon2SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Upsilon2SBrmumu(); 
 		break;
 	case 555:  // Upsilon(3S)
 		_particleType = UPSILON3S;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon3SMass();
 		width         = Upsilon3SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = (Upsilon3SBree() + Upsilon3SBrmumu())/2.; 
 		break;
 	case 555011:  // Upsilon(3S)
 		_particleType = UPSILON3S_ee;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon3SMass();
 		width         = Upsilon3SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Upsilon3SBree(); 
 		break;
 	case 555013:  // Upsilon(3S)
 		_particleType = UPSILON3S_mumu;
 		_decayType    = NARROWVMDEFAULT;
 		mass          = Upsilon3SMass();
 		width         = Upsilon3SWidth();
 		defaultMinW   = mass - 5 * width;
 		defaultMaxW   = mass + 5 * width;
 		_inputBranchingRatio = Upsilon3SBrmumu(); 
 		break;
 	default:
 		printWarn << "unknown particle ID " << _prodParticleId << endl;
 		return false;
 	}  // _prodParticleId
 
 	if (_minW.value() == -1)
 		_minW = defaultMinW;
 	if (_maxW.value() == -1)
 		_maxW = defaultMaxW;
 	if ( _maxW.value() <= _minW.value() ) {
 		printWarn << "maxW must be greater than minW" << endl;
 		return false;
 	}
 
 	printInfo << "using the following " << *this;
 	
 	return true;
 }
 
 
 //______________________________________________________________________________
 ostream&
 inputParameters::print(ostream& out) const
 {
 	out << "starlight parameters:" << endl
 	    << "    base file name  ...................... '"  << _baseFileName.value() << "'" << endl
 	    << "    beam 1 atomic number ................... " << _beam1Z.value() << endl
 	    << "    beam 1 atomic mass number .............. " << _beam1A.value() << endl
 	    << "    beam 2 atomic number ................... " << _beam2Z.value() << endl
 	    << "    beam 2 atomic mass number .............. " << _beam2A.value() << endl
 	    << "    Lorentz gamma of beams in CM frame ..... " << _beamLorentzGamma << endl
 	    << "    mass W of produced hadronic system ..... " << _minW.value() << " < W < " << _maxW.value() << " GeV/c^2" << endl
 	    << "    # of W bins ............................ " << _nmbWBins.value() << endl
 	    << "    maximum absolute value for rapidity .... " << _maxRapidity.value() << endl
 	    << "    # of rapidity bins ..................... " << _nmbRapidityBins.value() << endl
 	    << "    cut in pT............................... " << yesNo(_ptCutEnabled.value()) << endl;
     if (_ptCutEnabled.value()) {
 	out << "        minumum pT.......................... " << _ptCutMin.value() << " GeV/c" << endl
 	    << "        maximum pT.......................... " << _ptCutMax.value() << " GeV/c" << endl;}
 	out << "    cut in eta.............................. " << yesNo(_etaCutEnabled.value()) << endl;
     if (_etaCutEnabled.value()) {
 	out << "        minumum eta......................... " << _etaCutMin.value() << endl
 	    << "        maximum eta......................... " << _etaCutMax.value() << endl;}
         out << "    production mode ........................ " << _productionMode.value() << endl
 	    << "    number of events to generate ........... " << _nmbEventsTot.value() << endl
 	    << "    PDG ID of produced particle ............ " << _prodParticleId.value() << endl
 	    << "    seed for random generator .............. " << _randomSeed.value() << endl
 	    << "    breakup mode for beam particles ........ " << _beamBreakupMode.value() << endl
 	    << "    interference enabled ................... " << yesNo(_interferenceEnabled.value()) << endl;
     if (_interferenceEnabled.value()) {
 	out << "    interference strength .................. " << _interferenceStrength.value() << endl
 	    << "    maximum p_T for interference calc. ..... " << _maxPtInterference.value() << " GeV/c" << endl
 	    << "    # of p_T bins for interference calc. ... " << _nmbPtBinsInterference.value() << endl;}
     if (_productionMode.value()!=1) {
 		if (_productionMode.value()==4) {
 	 	  out  << "    coherent scattering off nucleus ........ no" << endl;}
 		else {
 	 	  out  << "    coherent scattering off nucleus ........ yes" << endl;
 		}
 	}
     out     <<"    Quantum Glauber parameter...............  "<<_quantumGlauber.value()<<endl;
     out     <<"    Impulse VM parameter....................  "<<_impulseVM.value()<<endl;
     //   if (_quantumGlauber.value()==1) {out << "    Quantum Glauber calculation being used"<< endl;}
     //if (_quantumGlauber.value()!=1) {out << "    Classical Glauber calculation being used"<< endl;}
     if (_beamBreakupMode.value()==8) {
       out <<"    Minimum impact parameter.................."<<_bmin.value()<<" fm"<<endl;
       out <<"    Maximum impact parameter.................."<<_bmax.value()<<" fm"<<endl;
     }
 
     // Add some checks here  SRK September, 2017
     if (_beamBreakupMode.value()==8 && _bmin.value() > _bmax.value()) {
       out <<"****Error bmin > bmax"<<endl;
       exit(-1);
     }
     if (_beamBreakupMode.value() == 8 && _productionMode.value() !=1) {
 	out <<"****Error.  Cannot use beam breakup mode 8 (with bmin and bmax) with photonuclear interactions"<<endl;
 	exit(-1);
     }
 	return out;
 }
 
 
 //______________________________________________________________________________
 ostream&
 inputParameters::write(ostream& out) const
 {
   
         out << "baseFileName"  << baseFileName         () <<endl
 	    << "BEAM_1_Z"      << beam1Z               () <<endl
 	    << "BEAM_2_Z"      << beam1A               () <<endl
 	    << "BEAM_1_A"      << beam2Z               () <<endl
 	    << "BEAM_2_A"      << beam2A               () <<endl
 	    << "BEAM_GAMMA"    << beamLorentzGamma     () <<endl
 	    << "W_MAX"         << maxW                 () <<endl
 	    << "W_MIN"         << minW                 () <<endl
 	    << "W_N_BINS"      << nmbWBins             () <<endl
 	    << "RAP_MAX"       << maxRapidity          () <<endl
 	    << "RAP_N_BINS"    << nmbRapidityBins      () <<endl
 	    << "CUT_PT"        << ptCutEnabled         () <<endl
 	    << "PT_MIN"        << ptCutMin             () <<endl
 	    << "PT_MAX"        << ptCutMax             () <<endl
 	    << "CUT_ETA"       << etaCutEnabled        () <<endl
 	    << "ETA_MIN"       << etaCutMin            () <<endl
 	    << "ETA_MAX"       << etaCutMax            () <<endl
 	    << "PROD_MODE"     << productionMode       () <<endl
 	    << "N_EVENTS"      << nmbEvents            () <<endl
 	    << "PROD_PID"      << prodParticleId       () <<endl
 	    << "RND_SEED"      << randomSeed           () <<endl
 	    << "BREAKUP_MODE"  << beamBreakupMode      () <<endl
 	    << "INTERFERENCE"  << interferenceEnabled  () <<endl
 	    << "IF_STRENGTH"   << interferenceStrength () <<endl
 	    << "INT_PT_MAX"    << maxPtInterference    () <<endl
 	    << "INT_PT_N_BINS" << nmbPtBinsInterference() <<endl
 	    << "QUANTUM_GLAUBER"<<quantumGlauber       () <<endl
 	    << "BMIN"          <<bmin                  ()<<endl
 	    << "BMAX"          <<bmax                  ()<<endl;
 
 	return out;
 }
 
 std::string
 inputParameters::parameterValueKey() const 
 {
   return parameterbase::_parameters.validationKey();
 }
Index: trunk/src/gammaavm.cpp
===================================================================
--- trunk/src/gammaavm.cpp	(revision 305)
+++ trunk/src/gammaavm.cpp	(revision 306)
@@ -1,949 +1,949 @@
 ///////////////////////////////////////////////////////////////////////////
 //
 //    Copyright 2010
 //
 //    This file is part of starlight.
 //
 //    starlight is free software: you can redistribute it and/or modify
 //    it under the terms of the GNU General Public License as published by
 //    the Free Software Foundation, either version 3 of the License, or
 //    (at your option) any later version.
 //
 //    starlight is distributed in the hope that it will be useful,
 //    but WITHOUT ANY WARRANTY; without even the implied warranty of
 //    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
 //    GNU General Public License for more details.
 //
 //    You should have received a copy of the GNU General Public License
 //    along with starlight. If not, see <http://www.gnu.org/licenses/>.
 //
 ///////////////////////////////////////////////////////////////////////////
 //
 // File and Version Information:
 // $Rev::                             $: revision of last commit
 // $Author::                          $: author of last commit
 // $Date::                            $: date of last commit
 //
 // Description:
 //    Added incoherent t2-> pt2 selection.  Following pp selection scheme
 //
 //
 ///////////////////////////////////////////////////////////////////////////
 
 
 #include <iostream>
 #include <fstream>
 #include <cassert>
 #include <cmath>
 
 #include "gammaavm.h"
 #include "photonNucleusCrossSection.h"
 #include "wideResonanceCrossSection.h"
 #include "narrowResonanceCrossSection.h"
 #include "incoherentVMCrossSection.h"
 
 using namespace std;
 
 
 //______________________________________________________________________________
 Gammaavectormeson::Gammaavectormeson(const inputParameters& inputParametersInstance, randomGenerator* randy, beamBeamSystem& bbsystem):eventChannel(inputParametersInstance, randy, bbsystem), _phaseSpaceGen(0)
 {
 	_VMNPT=inputParametersInstance.nmbPtBinsInterference();
 	_VMWmax=inputParametersInstance.maxW();
 	_VMWmin=inputParametersInstance.minW();
 	_VMYmax=inputParametersInstance.maxRapidity();
 	_VMYmin=-1.*_VMYmax;
 	_VMnumw=inputParametersInstance.nmbWBins();
 	_VMnumy=inputParametersInstance.nmbRapidityBins();
 	_VMgamma_em=inputParametersInstance.beamLorentzGamma();
 	_VMinterferencemode=inputParametersInstance.interferenceEnabled();
 	_VMbslope=0.;//Will define in wide/narrow constructor
         _bslopeDef=inputParametersInstance.bslopeDefinition();
 	_bslopeVal=inputParametersInstance.bslopeValue();
 	_pEnergy= inputParametersInstance.protonEnergy();
 	_VMpidtest=inputParametersInstance.prodParticleType();
 	_VMptmax=inputParametersInstance.maxPtInterference();
 	_VMdpt=inputParametersInstance.ptBinWidthInterference();
         _ProductionMode=inputParametersInstance.productionMode();
 
         N0 = 0; N1 = 0; N2 = 0; 
 	  if (_VMpidtest == starlightConstants::FOURPRONG){
 		// create n-body phase-spage generator
 		_phaseSpaceGen = new nBodyPhaseSpaceGen(_randy);
          }
 }
 
 
 //______________________________________________________________________________
 Gammaavectormeson::~Gammaavectormeson()
 {
 	if (_phaseSpaceGen)
 		delete _phaseSpaceGen;
 }
 
 
 //______________________________________________________________________________
 void Gammaavectormeson::pickwy(double &W, double &Y)
 {
         double dW, dY, xw,xy,xtest,btest;
 	int  IW,IY;
   
 	dW = (_VMWmax-_VMWmin)/double(_VMnumw);
 	dY = (_VMYmax-_VMYmin)/double(_VMnumy);
   
  L201pwy:
 
 	xw = _randy->Rndom();
 	W = _VMWmin + xw*(_VMWmax-_VMWmin);
 
 	if (W < 2 * _ip->pionChargedMass())
 		goto L201pwy;
   
 	IW = int((W-_VMWmin)/dW);
 	xy = _randy->Rndom();
 	Y = _VMYmin + xy*(_VMYmax-_VMYmin);
 	IY = int((Y-_VMYmin)/dY); 
 	xtest = _randy->Rndom();
 
 	if( xtest > _Farray[IW][IY] )
 		goto L201pwy;
 
         N0++; 
 	// Determine the target nucleus 
 	// For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2 
         if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
            if( _ProductionMode == 2 || _ProductionMode ==3){
 	     _TargetBeam = 2;
 	   } else {
              _TargetBeam = 1;
            }
         } else if(  _bbs.beam1().A() != 1 && _bbs.beam2().A()==1 ){
            if( _ProductionMode == 2 || _ProductionMode ==3){
 	     _TargetBeam = 1;
 	   } else {
              _TargetBeam = 2;
            }
         } else {
           btest = _randy->Rndom();
 	  if ( btest < _Farray1[IW][IY]/_Farray[IW][IY] ){
             _TargetBeam = 2;
             N2++;
           }  else {
             _TargetBeam = 1;
             N1++; 
           }
         }
 }         
 
 
 //______________________________________________________________________________                                               
 void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
                                      double  W,
                                      double  px0, double  py0, double  pz0,
                                      double& px1, double& py1, double& pz1,
                                      double& px2, double& py2, double& pz2,
                                      int&    iFbadevent)
 {
 	// This routine decays a particle into two particles of mass mdec,
 	// taking spin into account
 
 	double pmag;
 	double phi,theta,Ecm;
 	double betax,betay,betaz;
 	double mdec=0.0;
 	double E1=0.0,E2=0.0;
 
 	//    set the mass of the daughter particles
 	mdec=getDaughterMass(ipid);
 
 	//  calculate the magnitude of the momenta
 	if(W < 2*mdec){
 		cout<<" ERROR: W="<<W<<endl;
 		iFbadevent = 1;
 		return;
 	}
 	pmag = sqrt(W*W/4. - mdec*mdec);
   
 	//  pick an orientation, based on the spin
 	//  phi has a flat distribution in 2*pi
 	phi = _randy->Rndom()*2.*starlightConstants::pi;
                                                                                                                 
 	//  find theta, the angle between one of the outgoing particles and
 	//  the beamline, in the outgoing particles' rest frame
 
 	theta=getTheta(ipid);
  
 	//  compute unboosted momenta
 	px1 = sin(theta)*cos(phi)*pmag;
 	py1 = sin(theta)*sin(phi)*pmag;
 	pz1 = cos(theta)*pmag;
 	px2 = -px1;
 	py2 = -py1;
 	pz2 = -pz1;
 
 	Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
 	E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
 	E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
 
 	betax = -(px0/Ecm);
 	betay = -(py0/Ecm);
 	betaz = -(pz0/Ecm);
 
 	transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
 	transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
 
 	if(iFbadevent == 1)
 	   return;
 
 }
 
 
 //______________________________________________________________________________                                               
 // decays a particle into four particles with isotropic angular distribution
 bool Gammaavectormeson::fourBodyDecay
 (starlightConstants::particleTypeEnum& ipid,
  const double                  ,           // E (unused)
  const double                  W,          // mass of produced particle
  const double*                 p,          // momentum of produced particle; expected to have size 3
  lorentzVector*                decayVecs,  // array of Lorentz vectors of daughter particles; expected to have size 4
  int&                          iFbadevent)
 {
 	const double parentMass = W;
 
 	// set the mass of the daughter particles
 	const double daughterMass = getDaughterMass(ipid);
 	if (parentMass < 4 * daughterMass){
 		cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
 		iFbadevent = 1;
 		return false;
 	}
 
 	// construct parent four-vector
 	const double        parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
 	                                        + parentMass * parentMass);
 	const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
 
 	// setup n-body phase-space generator
 	assert(_phaseSpaceGen);
 	static bool firstCall = true;
 	if (firstCall) {
 		const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
 		_phaseSpaceGen->setDecay(4, m);
 		// estimate maximum phase-space weight
 		_phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
 		firstCall = false;
 	}
 
 	// generate phase-space event
 	if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
 		return false;
 
 	// set Lorentzvectors of decay daughters
 	for (unsigned int i = 0; i < 4; ++i)
 		decayVecs[i] = _phaseSpaceGen->daughter(i);
 	return true;
 }
 
 
 //______________________________________________________________________________
 double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
 {
 	//This will return the daughter particles mass, and the final particles outputed id...
 	double mdec=0.;
   
 	switch(_VMpidtest){
 	case starlightConstants::RHO:
 	case starlightConstants::RHOZEUS:
 	case starlightConstants::FOURPRONG:
 	case starlightConstants::OMEGA:
 		mdec = _ip->pionChargedMass();
 		ipid = starlightConstants::PION;
 		break;
 	case starlightConstants::PHI:
 		mdec = _ip->kaonChargedMass();
 		ipid = starlightConstants::KAONCHARGE;
 		break;
 	case starlightConstants::JPSI:
 		mdec = _ip->mel();
 		ipid = starlightConstants::ELECTRON;
 		break; 
 	case starlightConstants::JPSI_ee:
 		mdec = _ip->mel();
 		ipid = starlightConstants::ELECTRON;
 		break; 
 	case starlightConstants::JPSI_mumu:
 		mdec = _ip->muonMass();
 		ipid = starlightConstants::MUON;
 		break; 
 	case starlightConstants::JPSI_ppbar:
 		mdec = _ip->protonMass();
 		ipid = starlightConstants::PROTON;
 		break; 
 	case starlightConstants::JPSI2S_ee:
 		mdec = _ip->mel();
 		ipid = starlightConstants::ELECTRON;
 		break; 
 	case starlightConstants::JPSI2S_mumu:
 		mdec = _ip->muonMass();
 		ipid = starlightConstants::MUON;
 		break; 
 
 	case starlightConstants::JPSI2S:
 	case starlightConstants::UPSILON:
 	case starlightConstants::UPSILON2S:
 	case starlightConstants::UPSILON3S:
 		mdec = _ip->muonMass();
 		ipid = starlightConstants::MUON;
 		break;
 	case starlightConstants::UPSILON_ee:
 	case starlightConstants::UPSILON2S_ee:
 	case starlightConstants::UPSILON3S_ee:
 		mdec = _ip->mel();
 		ipid = starlightConstants::ELECTRON;
 		break;
 	case starlightConstants::UPSILON_mumu:
 	case starlightConstants::UPSILON2S_mumu:
 	case starlightConstants::UPSILON3S_mumu:
 		mdec = _ip->muonMass();
 		ipid = starlightConstants::MUON;   
 		break;
 	default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
 	}
   
 	return mdec;
 }
 
 
 //______________________________________________________________________________
 double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid)
 {
 	//This depends on the decay angular distribution
 	//Valid for rho, phi, omega.
 	double theta=0.;
 	double xtest=0.;
 	double dndtheta=0.;
 
  L200td:
     
 	theta = starlightConstants::pi*_randy->Rndom();
 	xtest = _randy->Rndom();
 	//  Follow distribution for helicity +/-1
 	//  Eq. 19 of J. Breitweg et al., Eur. Phys. J. C2, 247 (1998)
 	//  SRK 11/14/2000
   
 	switch(ipid){
 	  
 	case starlightConstants::MUON:
 	case starlightConstants::ELECTRON:
 		//VM->ee/mumu
  	        dndtheta = sin(theta)*(1.+(cos(theta)*cos(theta)));
                 break; 
 		
 	case starlightConstants::PROTON:
 		//Pick this angular distribution for J/psi --> ppbar 
 	        dndtheta = sin(theta)*(1.+(0.605*cos(theta)*cos(theta)));
 		break;
     
 	case starlightConstants::PION:
 	case starlightConstants::KAONCHARGE:
 		//rhos etc
 		dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
 		break;
     
 	default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
 	}//end of switch
   
 	if(xtest > dndtheta)
 		goto L200td;
   
 	return theta;
   
 }
 
 
 //______________________________________________________________________________
 double Gammaavectormeson::getWidth()
 {
 	return _width;
 }
 
 
 //______________________________________________________________________________
 double Gammaavectormeson::getMass()
 {
 	return _mass;
 }
 
 
 //______________________________________________________________________________
 double Gammaavectormeson::getSpin()
 {
 	return 1.0; //VM spins are the same
 }
 
 
 //______________________________________________________________________________
 void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)
 {
 	//     This subroutine calculates momentum and energy of vector meson
 	//     given W and Y,   without interference.  Subroutine vmpt handles
 	//     production with interference
  
 	double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
 	double px1,py1,px2,py2;
 	double pt,xt,xtest,ytest;
 	double t2;
 
   
 	//Find Egam,Epom in CM frame
         if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
           // This is pA
           if( _ProductionMode == 2 || _ProductionMode ==3 ){
     	    Egam = 0.5*W*exp(Y);
   	    Epom = 0.5*W*exp(-Y);
           }else{
     	    Egam = 0.5*W*exp(-Y);
   	    Epom = 0.5*W*exp(Y);
           }  
         } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
           // This is Ap
           if( _ProductionMode == 2 || _ProductionMode == 3 ){
   	    Egam = 0.5*W*exp(-Y);
   	    Epom = 0.5*W*exp(Y);
           }else{
     	    Egam = 0.5*W*exp(Y);
   	    Epom = 0.5*W*exp(-Y);
           }
 	} else {
           // This is pp or AA 
           if( _TargetBeam == 1 ){
             Egam = 0.5*W*exp(-Y);
 	    Epom = 0.5*W*exp(Y);
 	  }
           else {
             Egam = 0.5*W*exp(Y);
 	    Epom = 0.5*W*exp(-Y);
 	  }
 	}
 
 	//        } else if( _ProductionMode == 2 || _ProductionMode==3){
 	//	  Egam = 0.5*W*exp(-Y);
 	//	  Epom = 0.5*W*exp(Y);
 	//        } else { 
 	//          Egam = 0.5*W*exp(Y);
 	//	  Epom = 0.5*W*exp(-Y);
 	//	 }
 
         pt1 = pTgamma(Egam);  
 	phi1 = 2.*starlightConstants::pi*_randy->Rndom();
 
 	if( (_bbs.beam1().A()==1 && _bbs.beam2().A()==1) || 
             (_ProductionMode == 4) ) {
 	    if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
 	      // Use dipole form factor for light VM
 	      L555vm:
 	      xtest = 2.0*_randy->Rndom();
               double ttest = xtest*xtest; 
               ytest = _randy->Rndom();
               double t0 = 1./2.23; 
               double yprob = xtest*_bbs.beam1().dipoleFormFactor(ttest,t0)*_bbs.beam1().dipoleFormFactor(ttest,t0); 
               if( ytest > yprob ) goto L555vm; 
               t2 = ttest; 
               pt2 = xtest;              
 	    }else{
 		//Use dsig/dt= exp(-_VMbslope*t) for heavy VM
                 double bslope_tdist = _VMbslope; 
 		double Wgammap = 0.0; 
                 switch(_bslopeDef){
 		  case 0:
 		    //This is the default, as before
 		    bslope_tdist = _VMbslope;
 		    break;
 		  case 1:
 		    //User defined value of bslope. BSLOPE_VALUE default is 4.0 if not set. 
                     bslope_tdist = _bslopeVal;
 		    if( N0 <= 1 )cout<<" ATTENTION: Using user defined value of bslope = "<<_bslopeVal<<endl;
                     break; 
 		  case 2:
                     //This is Wgammap dependence of b from H1 (Eur. Phys. J. C 46 (2006) 585)
 		    Wgammap = sqrt(4.*Egam*_pEnergy); 
 		    bslope_tdist = 4.63 + 4.*0.164*log(Wgammap/90.0);
 		    if( N0 <= 1 )cout<<" ATTENTION: Using energy dependent value of bslope!"<<endl; 
 		    break;
 		  default:
 		    cout<<" Undefined setting for BSLOPE_DEFINITION "<<endl;
 		}
 
 	        xtest = _randy->Rndom(); 
 		// t2 = (-1./_VMbslope)*log(xtest);
 		t2 = (-1./bslope_tdist)*log(xtest);
 		pt2 = sqrt(1.*t2);
 	    }
 	} else {
 	    // >> Check tmin
 	    tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
 	
 	    if(tmin > 0.5){
 		cout<<" WARNING: tmin= "<<tmin<<endl;
                 cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl; 
 		cout<<" Will pick a new W,Y "<<endl;
 		tcheck = 1;
 		return;
 	    }
  L203vm:
 	    xt = _randy->Rndom(); 
             if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
               if( _ProductionMode == 2 || _ProductionMode ==3){
 		
  // Changed '8' to '32' 6 times below to extend the region of the p_T calculation up to 1 GeV.c  SRK May 28, 2019
 //  'pt2' is the maximum vector meson momentum.  For heavy nuclei, the '32'coefficient corresonds to about 1 GeV/c
 //  The downside of the larger coefficient is that the sampling runs more slowly.  This could be made into a parameter		
 		
 
 		pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();  
               }else{
                 pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();  
               }   
             } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
               if( _ProductionMode == 2 || _ProductionMode ==3){
                 pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();  
               }else{
                 pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();  
               }  
             } else if (_TargetBeam==1) {
                 pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();  
             } else {
                 pt2 = 32.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();  
             }
 
 	    xtest = _randy->Rndom();
 	    t2 = tmin + pt2*pt2;
 
 	    double comp=0.0; 
             if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
               if( _ProductionMode == 2 || _ProductionMode ==3){
                 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
               }else{
                 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
               }   
             } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
               if( _ProductionMode == 2 || _ProductionMode ==3){
                 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
               }else{
                 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
               }  
             } else if (_TargetBeam==1) {
               comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
             } else {
               comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2; 
             }	      
             if( xtest > comp ) goto L203vm;
        		
 	}//else end from pp
 
 	phi2 = 2.*starlightConstants::pi*_randy->Rndom();
 
 	px1 = pt1*cos(phi1);
 	py1 = pt1*sin(phi1);
 	px2 = pt2*cos(phi2);
 	py2 = pt2*sin(phi2);
         
 	// Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
 	px = px1 + px2;
 	py = py1 + py2;
 	pt = sqrt( px*px + py*py );
        
 	E  = sqrt(W*W+pt*pt)*cosh(Y);
 	pz = sqrt(W*W+pt*pt)*sinh(Y);
 
 }
 
 //______________________________________________________________________________
 double Gammaavectormeson::pTgamma(double E)
 {
     // returns on random draw from pp(E) distribution
     double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
     double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
     int satisfy =0;
         
     ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
     //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
     Cm = sqrt(3.)*E/_VMgamma_em;
     // If E is very small, the drawing of a pT below is extre_ip->mel()y slow. 
     // ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E. 
     // Should have no observable consequences (JN, SRK 11-Sep-2014)
     if( E < 0.0005 )return Cm; 
  
     //the amplitude of the p_t spectrum at the maximum
 
     if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
       if( _ProductionMode == 2 || _ProductionMode ==3 ){
          singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
       }else{
          singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
       }  
     } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
       if( _ProductionMode == 2 || _ProductionMode ==3){
          singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
       }else{
          singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
       }  
     } else if (_TargetBeam == 1) {
       singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
     } else {
       singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
     }
 
     Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
         
     //pick a test value pp, and find the amplitude there
     x = _randy->Rndom();
 
     if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
       if( _ProductionMode == 2 || _ProductionMode ==3){
         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
       }else{
         pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
       }  
     } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
       if( _ProductionMode == 2 || _ProductionMode ==3){
         pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
       }else{
         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
       }  
     } else if (_TargetBeam == 1) {
         pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
     } else {
         pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
         singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
     }
 
     test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
 
     while(satisfy==0){
 	u = _randy->Rndom();
 	if(u*Coef <= test)
 	{
 	    satisfy =1;
 	}
 	else{
 	    x =_randy->Rndom();
             if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){ 
               if( _ProductionMode == 2 || _ProductionMode ==3){
                 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
                 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
               }else{
                 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
                 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
               }  
             } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
               if( _ProductionMode == 2 || _ProductionMode ==3){
                 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
                 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
               }else{
                 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
                 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
               }  
             } else if (_TargetBeam == 1) {
               pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); 
               singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
             } else {
               pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); 
               singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
             }
 	    test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
 	}
     }
 
     return pp;
 }
 
 
 //______________________________________________________________________________
 void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
                              int&) // tcheck (unused)
 {
 	//    This function calculates momentum and energy of vector meson
 	//    given W and Y, including interference.
 	//    It gets the pt distribution from a lookup table.
 	double dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
 	int IY=0,j=0;
   
 	dY  = (_VMYmax-_VMYmin)/double(_VMnumy);
   
 	//  Y is already fixed; choose a pt
 	//  Follow the approach in pickwy
 	//  in  _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
  	//  Changed,  now works -y to +y.
 	IY=int((Y-_VMYmin)/dY);
 	if (IY > (_VMnumy)-1){
         	IY=(_VMnumy)-1;
 	}
 
 	yleft=(Y-_VMYmin)-(IY)*dY;
 
 	yfract=yleft*dY;
   
 	xpt=_randy->Rndom();
 	for(j=0;j<_VMNPT;j++){
 		if (xpt < _fptarray[IY][j]) goto L60;
 	}
 	if(j == _VMNPT) j = _VMNPT-1;
  L60:
   
 	//  now do linear interpolation - start with extremes
   	if (j == 0){
 		pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
 		goto L80;
 	}
 	if (j == _VMNPT-1){
 		pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
 		goto L80;
 	}
   
 	//  we're in the middle
   	ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
 	pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
   
 	//  at an extreme in y?
 	if (IY == (_VMnumy)-1){
 		pt=pt1;
 		goto L120;
 	}
  L80:
 
 	//  interpolate in y repeat for next fractional y bin      
 	for(j=0;j<_VMNPT;j++){
 		if (xpt < _fptarray[IY+1][j]) goto L90;
 	}
         if(j == _VMNPT) j = _VMNPT-1;
  L90:
   
 	//  now do linear interpolation - start with extremes
 	if (j == 0){
 		pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
 		goto L100;
 	}
 	if (j == _VMNPT-1){
 		pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
 		goto L100;
 	}
   
 	//  we're in the middle
 	ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
 	pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
  L100:
 
 	//  now interpolate in y  
 	pt=yfract*pt2+(1-yfract)*pt1;
  L120:
 
 	//  we have a pt 
 	theta=2.*starlightConstants::pi*_randy->Rndom();
 	px=pt*cos(theta);
 	py=pt*sin(theta);
 
 	E  = sqrt(W*W+pt*pt)*cosh(Y);
 	pz = sqrt(W*W+pt*pt)*sinh(Y);
 	//      randomly choose to make pz negative 50% of the time
 	if(_randy->Rndom()>=0.5) pz = -pz;
 }
 
 
 //______________________________________________________________________________
 starlightConstants::event Gammaavectormeson::produceEvent(int&)
 {
 	//Note used; return default event
 	return starlightConstants::event();
 }
 
 
 //______________________________________________________________________________
 upcEvent Gammaavectormeson::produceEvent()
 {
 	// The new event type
 	upcEvent event;
 
 	int iFbadevent=0;
 	int tcheck=0;
 	starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
         starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN; 
 
 	if (_VMpidtest == starlightConstants::FOURPRONG) {
 		double        comenergy = 0;
 		double        mom[3]    = {0, 0, 0};
 		double        E         = 0;
 		lorentzVector decayVecs[4];
 		do {
 			double rapidity = 0;
 			pickwy(comenergy, rapidity);
 			if (_VMinterferencemode == 0)
 				momenta(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
 			else if (_VMinterferencemode==1)
 				vmpt(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
 		} while (!fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent));
 		if ((iFbadevent == 0) and (tcheck == 0))
 			for (unsigned int i = 0; i < 4; ++i) {
 				starlightParticle daughter(decayVecs[i].GetPx(),
 				                           decayVecs[i].GetPy(),
 				                           decayVecs[i].GetPz(),
-				                           starlightConstants::UNKNOWN,  // energy 
-				                           starlightConstants::UNKNOWN,  // _mass
-				                           ipid,
-				                           (i < 2) ? -1 : +1);
+							   sqrt(decayVecs[i].GetPx()*decayVecs[i].GetPx()+decayVecs[i].GetPy()*decayVecs[i].GetPy()+decayVecs[i].GetPz()*decayVecs[i].GetPz()+0.139*0.139),//energy 
+							   starlightConstants::UNKNOWN,  // _mass
+							   ipid*(2*(i/2)-1),   // make half of the particles pi^+, half pi^-
+							   i/2);
 				event.addParticle(daughter);
 			}
 	} else {
 		double comenergy = 0.;
 		double rapidity = 0.;
 		double E = 0.;
 		double momx=0.,momy=0.,momz=0.;
 
 		double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
 		bool accepted = false;
 		do{
 			pickwy(comenergy,rapidity);
 
 			if (_VMinterferencemode==0){
 				momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
 			
 			} else if (_VMinterferencemode==1){
 				vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck);
 			}
 	   
 			_nmbAttempts++;
 
                         vmpid = ipid;
 			twoBodyDecay(ipid,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
 			double pt1chk = sqrt(px1*px1+py1*py1);
 			double pt2chk = sqrt(px2*px2+py2*py2);
 			double eta1 = pseudoRapidity(px1, py1, pz1);
 			double eta2 = pseudoRapidity(px2, py2, pz2);
 
 			if(_ptCutEnabled && !_etaCutEnabled){
 				if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
 					accepted = true;
 					_nmbAccepted++;
 				}
 			}
 			else if(!_ptCutEnabled && _etaCutEnabled){
 				if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
 					accepted = true;
 					_nmbAccepted++;
 				}
 			}
 			else if(_ptCutEnabled && _etaCutEnabled){
 				if(pt1chk > _ptCutMin && pt1chk < _ptCutMax &&  pt2chk > _ptCutMin && pt2chk < _ptCutMax){
 					if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
 						accepted = true;
 						_nmbAccepted++;
 					}
 				}
 			}
 			else if(!_ptCutEnabled && !_etaCutEnabled)
 				_nmbAccepted++;
 		}while((_ptCutEnabled || _etaCutEnabled) && !accepted);
 		if (iFbadevent==0&&tcheck==0) {
 			int q1=0,q2=0;
                         int ipid1,ipid2=0;
 
 			double xtest = _randy->Rndom(); 
 			if (xtest<0.5)
 				{
 					q1=1;
 					q2=-1;
 				}
 			else {
 				q1=-1;
 				q2=1;
 			}
 
                         if ( ipid == 11 || ipid == 13 ){
                           ipid1 = -q1*ipid;
                           ipid2 = -q2*ipid;
                         } else {
                           ipid1 = q1*ipid;
                           ipid2 = q2*ipid;
                         }
 
 			double md = getDaughterMass(vmpid);
                         double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1); 
 			starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
 			event.addParticle(particle1);
 
                         double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2); 
 			starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
 			event.addParticle(particle2);
 
 
 		}
 	}
 
 	return event;
 
 }
 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
 {
 	double pT = sqrt(px*px + py*py);
 	double p = sqrt(pz*pz + pT*pT);
 	double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
 	return eta;
 }
 
 //______________________________________________________________________________
 Gammaanarrowvm::Gammaanarrowvm(const inputParameters& input, randomGenerator* randy, beamBeamSystem& bbsystem):Gammaavectormeson(input, randy, bbsystem)
 {
 	cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
 	read();
 	cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
 	narrowResonanceCrossSection sigma(input, bbsystem);
 	sigma.crossSectionCalculation(_bwnormsave);
 	setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
 	_VMbslope=sigma.slopeParameter(); 
 }
 
 
 //______________________________________________________________________________
 Gammaanarrowvm::~Gammaanarrowvm()
 { }
 
 
 //______________________________________________________________________________
 Gammaaincoherentvm::Gammaaincoherentvm(const inputParameters& input, randomGenerator* randy, beamBeamSystem& bbsystem):Gammaavectormeson(input, randy, bbsystem)
 {
         cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
         read();
         cout<<"Creating and calculating crosssection. Gammaaincoherentvm()"<<endl;
         incoherentVMCrossSection sigma(input, bbsystem);
 	sigma.crossSectionCalculation(_bwnormsave);
 	setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
         _VMbslope=sigma.slopeParameter(); 
 }
 
 
 //______________________________________________________________________________
 Gammaaincoherentvm::~Gammaaincoherentvm()
 { }
 
 
 //______________________________________________________________________________
 Gammaawidevm::Gammaawidevm(const inputParameters& input, randomGenerator* randy, beamBeamSystem& bbsystem):Gammaavectormeson(input, randy, bbsystem)
 {
 	cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
 	read();
 	cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
 	wideResonanceCrossSection sigma(input, bbsystem);
 	sigma.crossSectionCalculation(_bwnormsave);
 	setTotalChannelCrossSection(sigma.getPhotonNucleusSigma());
 	_VMbslope=sigma.slopeParameter();
 }
 
 
 //______________________________________________________________________________
 Gammaawidevm::~Gammaawidevm()
 { }
Index: trunk/Readme.docx
===================================================================
Cannot display: file marked as a binary type.
svn:mime-type = application/octet-stream