Implemented more accurate nuclear density models for some light nuclei and added the ρ′(1450) meson production and decay to eSTARlight.

include/nucleus.h

@@ -66,21 +66,46 @@ class nucleus
 
 	double nuclearRadius() const { return _Radius; }
 	double rho0() const { return _rho0; }
+	double woodSaxonSkinDepth() const { return _woodSaxonSkinDepth;}
+	/** The "wine_bottle" parameter
+	 */
+	double w() const {return _w;}
+
 
 private:
 
-	double woodSaxonSkinDepth() const { return 0.53; } // 0.53 fm skin depth
-	double rws(const double r) const;
-
+	 
+
+    double rws(const double r) const;
 	int    _Z;                      ///< atomic number of nucleus
 	int    _A;                      ///< nucleon number of nucleus
-        int    _productionMode;
+    int    _productionMode;
 
 	double _r0;
 	double _Radius;
 	double _rho0;
-
+	double _w;               ///< 3pF fermi model parameter
+	double _woodSaxonSkinDepth;
+	double _v; ///< Sum of Gaussian(SOG) model parameter 
+	double *_a;  // SOG model parameter
+        double *_b;  // SOG model parameter
+        double *_c;  // SOG model parameter
+
+//SOG model parameters https:doi.org/10.1016/0092-640X(87)90013-1
+
+	double pkCarb[12] = {0.009560250143852663, 0.02161993792818375, 0.023894326875090147, 0.022871610300732596, 0.01787510497352304, 0.013266599310805718, 0.0020798983410400215, 0.0012487276970433168, 0.0001121651485588463, 0.00001856133380347859, 0.00000004241988368721027, 0.0};
+        double okCarb[12] = {0.016690, 0.050325, 0.128621, 0.180515, 0.219097, 0.278416, 0.058779, 0.057817, 0.007739, 0.002001, 0.000007, 0.0};
+        double qkCarb[12] = {0.0, 0.4, 1.0, 1.3, 1.7, 2.3, 2.7, 3.5, 4.3, 5.4, 6.7, 0.0};
+
+	double rkHel[12] = {0.2, 0.6, 0.9, 1.4, 1.9, 2.3, 2.6, 3.1, 3.5, 4.2, 4.9, 5.2};
+        double tkHel[12] = {0.034724, 0.430761, 0.203166, 0.192986, 0.083866, 0.033007, 0.014201, 0.0, 0.006860, 0.0, 0.000438, 0.0};
+	double skHel[12] = {0.010229060308680777, 0.0683272860306785, 0.01954234673886537, 0.009254567928457054, 0.002338942473693466, 0.0006455191775197494, 0.00022017367746850975, 0.0, 0.000059954958779193, 0.0, 0.000001978748821697672, 0.0};
 };
 
 
 #endif  // NUCLEUS_H
+
+
+
+
+

include/starlightconstants.h

@@ -102,8 +102,8 @@ namespace starlightConstants
         static const double rho0Mass  = 0.769;               // [GeV/c^2]
 	static const double rho0Width = 0.1517;              // [GeV/c^2]
 	static const double rho0BrPiPi = 1.0;                // Branching ratio pi+pi-
-        static const double rho0PrimeMass  = 1.540;          // [GeV/c^2]
-	static const double rho0PrimeWidth = 0.570;          // [GeV/c^2]
+        static const double rho0PrimeMass  = 1.465;          // [GeV/c^2]
+	static const double rho0PrimeWidth = 0.4;          // [GeV/c^2]
 	static const double rho0PrimeBrPiPi = 1.0;           // Branching ratio pi+pi- (set to 100%)
         static const double OmegaMass  = 0.78265;            // [GeV/c^2]
 	static const double OmegaWidth = 0.00849;            // [GeV/c^2]

src/e_wideResonanceCrossSection.cpp

@@ -125,6 +125,11 @@ e_wideResonanceCrossSection::crossSectionCalculation(const double bwnormsave)
 		W = _wideWmin + double(iW)*dW + 0.5*dW;
 		int nQ2 = 1000;
 		for(iEgamma = 0 ; iEgamma < nEgamma; ++iEgamma){    // Integral over photon energy
+		  // Displaying the percentage progress
+		  float ratio = float(iW)/float(nW)  + 1/float(nW)*float(iEgamma)/float(nEgamma);
+		 printf("calculating cross section :%3.2f %%\r",float(ratio *100.0));
+
+
 		  // Target frame photon energies
 		  ega[0] = exp(minEgamma + iEgamma*dEgamma );
 		  ega[1] = exp(minEgamma + (iEgamma+1)*dEgamma );

src/gammaavm.cpp

@@ -35,6 +35,7 @@
 #include <fstream>
 #include <cassert>
 #include <cmath>
+#include <sstream>
 
 #include "gammaavm.h"
 //#include "photonNucleusCrossSection.h"
@@ -1124,8 +1125,11 @@ eXEvent Gammaavectormeson::e_produceEvent()
 	int iFbadevent=0;
 	int tcheck=0;
 	starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
-        starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN; 
-	// at present 4 prong decay is not implemented
+    starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN; 
+	// at present 4 prong decay is implemented (rhoprime(1450))
+	double ptCutMin2 = _ptCutMin*_ptCutMin;//used for ptCut comparison without using square roots - to reduce processing time
+	double ptCutMax2 = _ptCutMax*_ptCutMax;//same as above
+
 	double comenergy = 0.;
 	double rapidity = 0.;
 	double Q2 = 0;
@@ -1137,6 +1141,103 @@ eXEvent Gammaavectormeson::e_produceEvent()
 	double e_E=0., e_phi=0;
 	double t_px =0, t_py=0., t_pz=0, t_E;
 	bool accepted = false;
+
+	if (_VMpidtest == starlightConstants::FOURPRONG) {
+	double mom[3] ={0,0,0};
+	lorentzVector decayVecs[4];
+	bool accepted;
+	double rapidity = 0;
+	do {
+	tcheck = 0;//reinitialized after every loop cycle - to avoid infinite loop
+	iFbadevent = 0;//same as above.
+	//pickwy(comenergy, rapidity); 
+	pickwEgamq2(comenergy,cmsEgamma, targetEgamma, 
+		   Q2, gamma_pz, gamma_pt, //photon infor in CMS frame
+		   e_E, e_theta);
+
+	//Vector meson is created and its four momentum is determined below
+	//if (_VMinterferencemode == 0)
+	 momenta(comenergy,cmsEgamma, Q2, gamma_pz, gamma_pt, //input
+		  rapidity, E, mom[0], mom[1], mom[2], //VM
+		  t_px, t_py, t_pz, t_E, //pomeron
+		  e_phi,tcheck);
+
+	 _nmbAttempts++;
+	 accepted = true;//re-initialized after every loop cycle -to avoid infinite loop
+
+	if(tcheck != 0 || !fourBodyDecay(ipid,E,comenergy,mom,decayVecs,iFbadevent))
+	{//if either vector meson creation, or further decay into four pions, is impossible
+	accepted = false;
+	continue;//this skips the etaCut and ptCut checks.
+	}
+
+	if (_ptCutEnabled) {
+	for (int i = 0; i < 4; i++) {
+	double pt_chk2 = 0;
+	pt_chk2 += pow( decayVecs[i].GetPx() , 2);
+	pt_chk2 += pow( decayVecs[i].GetPy() , 2);// compute transverse momentum (squared) for the particle, for ptCut checks.
+
+	if (pt_chk2 < ptCutMin2 || pt_chk2 > ptCutMax2) {//if particle does not fall into ptCut range.
+	accepted = false;
+	break;//skips checking other daughter particles
+	}
+	}
+	}
+	if (_etaCutEnabled) {
+	for (int i = 0; i < 4; i++) {
+	double eta_chk = pseudoRapidity(
+	decayVecs[i].GetPx(),
+	decayVecs[i].GetPy(),
+	decayVecs[i].GetPz()
+	);
+	//computes the pseudorapidity in the laboratory frame.
+	if (eta_chk < _etaCutMin || eta_chk > _etaCutMax) {//if this particle does not fall into range
+	accepted = false;
+	break;//skips checking other daughter particles
+	}
+	}
+	}
+	if (accepted and (tcheck == 0)) {
+	_nmbAccepted++;//maintain counts of accepted events.
+	}
+
+	} while (!accepted || tcheck != 0);//repeats loop if VM creation, decay, ptcut or etaCut criterias are not fulfilled. Important to avoid situations where events produced is less than requested.
+	if (iFbadevent==0&&tcheck==0) {
+	  double e_ptot = sqrt(e_E*e_E - starlightConstants::mel*starlightConstants::mel);
+	  double e_px = e_ptot*sin(e_theta)*cos(e_phi);
+	  double e_py = e_ptot*sin(e_theta)*sin(e_phi);
+	  double e_pz = e_ptot*cos(e_theta);
+	  lorentzVector electron(e_px, e_py, e_pz, e_E);
+	  event.addSourceElectron(electron);
+	  // - Generated photon - target frame
+	  double gamma_x = gamma_pt*cos(e_phi+starlightConstants::pi);
+	  double gamma_y = gamma_pt*sin(e_phi+starlightConstants::pi);
+	  lorentzVector gamma(gamma_x,gamma_y,gamma_pz,cmsEgamma);
+	  vector3 boostVector(0, 0, tanh(_rap_CM));
+	  (gamma).Boost(boostVector);
+	  event.addGamma(gamma, targetEgamma, Q2);
+
+	double md = getDaughterMass(ipid);
+	//adds daughters as particles into the output event.
+	for (unsigned int i = 0; i < 4; ++i) {
+	starlightParticle daughter(decayVecs[i].GetPx(),
+	                          decayVecs[i].GetPy(),
+	                          decayVecs[i].GetPz(),
+	  sqrt(decayVecs[i].GetPx()*decayVecs[i].GetPx()+decayVecs[i].GetPy()*decayVecs[i].GetPy()+decayVecs[i].GetPz()*decayVecs[i].GetPz()+md*md),//energy
+	  md,  // _mass
+	  ipid*(2*(i/2)-1),   // make half of the particles pi^+, half pi^-
+	  (2*(i/2)-1));//charge
+	event.addParticle(daughter);
+	}
+	// - Scattered target and transfered momenta at target vertex
+	  double target_pz =  - _beamNucleus*sqrt(_pEnergy*_pEnergy - pow(starlightConstants::protonMass,2.) ) - t_pz;
+	  //Sign of t_px in following equation changed to fix sign error and conserve p_z.  Change made by Spencer Klein based on a bug report from Ya-Ping Xie.  Nov. 14, 2019
+	  lorentzVector target(-t_px, -t_py, target_pz, _beamNucleus*_pEnergy - t_E);
+	  double t_var = t_E*t_E - t_px*t_px - t_py*t_py - t_pz*t_pz;
+	  event.addScatteredTarget(target, t_var);
+	}
+	}
+else {
 	do{
 	  pickwEgamq2(comenergy,cmsEgamma, targetEgamma, 
 		   Q2, gamma_pz, gamma_pt, //photon infor in CMS frame
@@ -1259,6 +1360,7 @@ eXEvent Gammaavectormeson::e_produceEvent()
 	  double t_var = t_E*t_E - t_px*t_px - t_py*t_py - t_pz*t_pz;
 	  event.addScatteredTarget(target, t_var);
 	}
+}
 	return event;
 
 }

src/inputParameters.cpp

@@ -386,15 +386,16 @@ inputParameters::configureFromFile(const std::string &_configFileName)
 		defaultMaxW         = mass + 5 * width;  // use the same 1.5GeV max mass as ZEUS
 		_inputBranchingRatio = starlightConstants::rho0BrPiPi; 
 		break;
-		//	case 999:  // pi+pi-pi+pi- phase space decay
-		//		_particleType = FOURPRONG;
-		//		_decayType    = WIDEVMDEFAULT;
-		//		mass          = starlightConstants::rho0PrimeMass;
-		//		width         = starlightConstants::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
-		//		_inputBranchingRatio = 1.0; 
-		//		break;
+	case 999:  // rho'(1450)-> pi+pi-pi+pi- phase space decay
+		_particleType = FOURPRONG;
+		_decayType    = WIDEVMDEFAULT;
+		mass          = starlightConstants::rho0PrimeMass;
+		width         = starlightConstants::rho0PrimeWidth;		
+		defaultMinW   = 4 * pionChargedMass;
+                defaultMaxW   = mass + 5 * width;	    
+    	        _inputBranchingRatio = starlightConstants::rho0PrimeBrPiPi; 
+		break;
+
 	case 223:  // omega(782)
 		_particleType = OMEGA;
 		_decayType    = NARROWVMDEFAULT;