+++ /dev/null
-///////////////////////////////////////////////////////////////////////////
-//
-// 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(beamBeamSystem& bbsystem):eventChannel(bbsystem), _phaseSpaceGen(0) //:readLuminosity(input),_bbs(bbsystem)
-{
- _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
- _VMpidtest=inputParametersInstance.prodParticleType();
- _VMptmax=inputParametersInstance.maxPtInterference();
- _VMdpt=inputParametersInstance.ptBinWidthInterference();
- _VMCoherence=inputParametersInstance.coherentProduction();
- _VMCoherenceFactor=inputParametersInstance.coherentProduction();
- _ProductionMode=inputParametersInstance.productionMode();
-
- switch(_VMpidtest){
- case starlightConstants::RHO:
- case starlightConstants::RHOZEUS:
- _width=0.1507;
- _mass=0.7685;
- break;
- case starlightConstants::FOURPRONG:
- // create n-body phase-space generator instance
- _phaseSpaceGen = new nBodyPhaseSpaceGen();
- _width = 0.360;
- _mass = 1.350;
- break;
- case starlightConstants::OMEGA:
- _width=0.00843;
- _mass=0.78194;
- break;
- case starlightConstants::PHI:
- _width=0.00443;
- _mass=1.019413;
- break;
- case starlightConstants::JPSI:
- case starlightConstants::JPSI_ee:
- case starlightConstants::JPSI_mumu:
- _width=0.000091;
- _mass=3.09692;
- break;
- case starlightConstants::JPSI2S:
- case starlightConstants::JPSI2S_ee:
- case starlightConstants::JPSI2S_mumu:
- _width=0.000337;
- _mass=3.686093;
- break;
- case starlightConstants::UPSILON:
- case starlightConstants::UPSILON_ee:
- case starlightConstants::UPSILON_mumu:
- _width=0.00005402;
- _mass=9.46030;
- break;
- case starlightConstants::UPSILON2S:
- case starlightConstants::UPSILON2S_ee:
- case starlightConstants::UPSILON2S_mumu:
- _width=0.00003198;
- _mass=10.02326;
- break;
- case starlightConstants::UPSILON3S:
- case starlightConstants::UPSILON3S_ee:
- case starlightConstants::UPSILON3S_mumu:
- _width=0.00002032;
- _mass=10.3552;
- break;
- default: cout<<"No proper vector meson defined, gammaavectormeson::gammaavectormeson"<<endl;
- }
-}
-
-
-//______________________________________________________________________________
-Gammaavectormeson::~Gammaavectormeson()
-{
- if (_phaseSpaceGen)
- delete _phaseSpaceGen;
-}
-
-
-//______________________________________________________________________________
-void Gammaavectormeson::pickwy(double &W, double &Y)
-{
- double dW, dY, xw,xy,xtest;
- int IW,IY;
-
- dW = (_VMWmax-_VMWmin)/double(_VMnumw);
- dY = (_VMYmax-_VMYmin)/double(_VMnumy);
-
- L201pwy:
-
- xw = randyInstance.Rndom();// random()/(RAND_MAX+1.0);
- W = _VMWmin + xw*(_VMWmax-_VMWmin);
-
- if (W < 2 * starlightConstants::pionChargedMass)
- goto L201pwy;
-
- IW = int((W-_VMWmin)/dW); //+ 1;
- xy = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
- Y = _VMYmin + xy*(_VMYmax-_VMYmin);
- IY = int((Y-_VMYmin)/dY); //+ 1;
- xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-
- if( xtest > _Farray[IW][IY] )
- goto L201pwy;
-
-}
-
-
-//______________________________________________________________________________
-void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
- double, // E (unused)
- 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 anglelep[20001],xtest,ytest=0.,dndtheta;
- 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 = randyInstance.Rndom()*2.*starlightConstants::pi;//(random()/(RAND_MAX+1.0))* 2.*pi;
-
- // find theta, the angle between one of the outgoing particles and
- // the beamline, in the frame of the two photons
-
- 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 ytest=0.,mdec=0.;
-
- switch(_VMpidtest){
- case starlightConstants::RHO:
- case starlightConstants::RHOZEUS:
- case starlightConstants::FOURPRONG:
- case starlightConstants::OMEGA:
- mdec = starlightConstants::pionChargedMass;
- ipid = starlightConstants::PION;
- break;
- case starlightConstants::PHI:
- mdec = starlightConstants::kaonChargedMass;
- ipid = starlightConstants::KAONCHARGE;
- break;
- case starlightConstants::JPSI:
- mdec = starlightConstants::mel;
- ipid = starlightConstants::ELECTRON;
- break;
- case starlightConstants::JPSI_ee:
- mdec = starlightConstants::mel;
- ipid = starlightConstants::ELECTRON;
- break;
- case starlightConstants::JPSI_mumu:
- mdec = starlightConstants::muonMass;
- ipid = starlightConstants::MUON;
- break;
- case starlightConstants::JPSI2S_ee:
- mdec = starlightConstants::mel;
- ipid = starlightConstants::ELECTRON;
- break;
- case starlightConstants::JPSI2S_mumu:
- mdec = starlightConstants::muonMass;
- ipid = starlightConstants::MUON;
- break;
-
- case starlightConstants::JPSI2S:
- case starlightConstants::UPSILON:
- case starlightConstants::UPSILON2S:
- case starlightConstants::UPSILON3S:
- // decays 50% to e+/e-, 50% to mu+/mu-
- ytest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-
- mdec = starlightConstants::muonMass;
- ipid = starlightConstants::MUON;
- break;
- case starlightConstants::UPSILON_ee:
- case starlightConstants::UPSILON2S_ee:
- case starlightConstants::UPSILON3S_ee:
- mdec = starlightConstants::mel;
- ipid = starlightConstants::ELECTRON;
- break;
- case starlightConstants::UPSILON_mumu:
- case starlightConstants::UPSILON2S_mumu:
- case starlightConstants::UPSILON3S_mumu:
- mdec = starlightConstants::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*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
- xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
- // 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:
- //primarily for upsilon/j/psi. VM->ee/mumu
- dndtheta = sin(theta)*(1.+((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.f handles
- // production with interference
-
- double dW,dY;
- double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
- double px1,py1,px2,py2;
- double pt,xt,xtest,ytest;
- // double photon_spectrum;
- double t2;
-
- dW = (_VMWmax-_VMWmin)/double(_VMnumw);
- dY = (_VMYmax-_VMYmin)/double(_VMnumy);
-
- //Find Egam,Epom in CM frame
- if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
- if( _ProductionMode == 2 ){
- 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 ){
- if( _ProductionMode == 2 ){
- 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 {
- Egam = 0.5*W*exp(Y);
- Epom = 0.5*W*exp(-Y);
- }
-
- // cout<<" Y: "<<Y<<" W: "<<W<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl;
- pt1 = pTgamma(Egam);
- phi1 = 2.*starlightConstants::pi*randyInstance.Rndom();
-
- // cout<<" pt1: "<<pt1<<endl;
-
- 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*randyInstance.Rndom();
- double ttest = xtest*xtest;
- ytest = randyInstance.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
- xtest = randyInstance.Rndom();
- t2 = (-1./_VMbslope)*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 = randyInstance.Rndom();
- if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
- if( _ProductionMode == 2 ){
- pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
- }else{
- pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
- }
- } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
- if( _ProductionMode == 2 ){
- pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
- }else{
- pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
- }
- } else {
- pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
- }
-
- xtest = randyInstance.Rndom();
- t2 = tmin + pt2*pt2;
-
- if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
- if(1.0 < _bbs.beam2().formFactor(t2)*pt2) cout <<"POMERON"<<endl;
- if( xtest > _bbs.beam2().formFactor(t2)*pt2) goto L203vm;
- }
- else{
-
- double comp=0.0;
- if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
- if( _ProductionMode == 2 ){
- 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 ){
- comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
- }else{
- comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
- }
- } else {
- comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
- }
-
- if( xtest > comp ) goto L203vm;
- }
-
- /*
- if(_VMCoherence==0 && (!(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2))){
- //dsig/dt= exp(-_VMbslope*t)
- xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
- t2 = (-1./_VMbslope)*log(xtest);
- pt2 = sqrt(1.*t2);
- }
- */
-
- }//else end from pp
- phi2 = 2.*starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
-
- // cout<<" pt2: "<<pt2<<endl;
-
- 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);
-
- // cout<< " Y = "<<Y<<" W = "<<W<<" Egam = "<<Egam<<" gamma = "<<_VMgamma_em<<endl;
-
- // Randomly choose to make pz negative 50% of the time
- if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
- pz = -pz;
- }else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
- // Don't switch
- }
- else{
- if (randyInstance.Rndom() >= 0.5) pz = -pz;
- }
-
-}
-
-//______________________________________________________________________________
-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;
-
- //the amplitude of the p_t spectrum at the maximum
-
- if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
- if( _ProductionMode == 2 ){
- 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 ){
- singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
- }else{
- singleformfactorCm=_bbs.beam1().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 = randyInstance.Rndom();
-
- if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
- if( _ProductionMode == 2 ){
- 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 ){
- 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 {
- 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 = randyInstance.Rndom();
- if(u*Coef <= test)
- {
- satisfy =1;
- }
- else{
- x =randyInstance.Rndom();
- if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
- if( _ProductionMode == 2 ){
- 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 ){
- 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 {
- 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 dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
- int IY=0,j=0;
-
- dW = (_VMWmax-_VMWmin)/double(_VMnumw);
- dY = (_VMYmax-_VMYmin)/double(_VMnumy);
-
- // Y is already fixed; choose a pt
- // Follow the approavh in pickwy.f
- // in _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
-
- IY=int(fabs(Y)/dY);//+1;
- if (IY > (_VMnumy/2)-1){
- IY=(_VMnumy/2)-1;
- }
-
- yleft=fabs(Y)-(IY)*dY;
- yfract=yleft*dY;
-
- xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
-
- for(j=0;j<_VMNPT+1;j++){
- if (xpt < _fptarray[IY][j]) goto L60;
- }
- L60:
-
- // now do linear interpolation - start with extremes
-
- if (j == 0){
- pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
- goto L80;
- }
- if (j == _VMNPT){
- 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/2)-1){
- pt=pt1;
- goto L120;
- }
- L80:
- // interpolate in y repeat for next fractional y bin
-
- for(j=0;j<_VMNPT+1;j++){
- if (xpt < _fptarray[IY+1][j]) goto L90;
- }
- L90:
-
- // now do linear interpolation - start with extremes
-
- if (j == 0){
- pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
- goto L100;
- }
- if (j == _VMNPT){
- 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*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
- px=pt*cos(theta);
- py=pt*sin(theta);
-
- // I guess W is the mass of the vector meson (not necessarily
- // on-mass-shell), and E is the energy
-
- 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(randyInstance.Rndom()>=0.5) pz = -pz;
-}
-
-
-//______________________________________________________________________________
-starlightConstants::event Gammaavectormeson::produceEvent(int&)
-{
- // Not 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);
- 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;
- // if(_accCut){
- 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);
- }
-
- // cout << "_ptCutMin: " << _ptCutMin << " _ptCutMax: " << _ptCutMax << " _etaCutMin: " << _etaCutMin << " _etaCutMax: " << _etaCutMax << endl;
- _nmbAttempts++;
- //cout << "n tries: " << _nmbAttempts<< endl;
- vmpid = ipid;
- twoBodyDecay(ipid,E,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);
-
- //cout << "pt1: " << pt1chk << " pt2: " << pt2chk << endl;
- double eta1 = pseudoRapidity(px1, py1, pz1);
- double eta2 = pseudoRapidity(px2, py2, pz2);
- //cout << "eta1: " << eta1 << " eta2: " << eta2 << endl;
- 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);
- /* }else{
- twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
- }*/
- if (iFbadevent==0&&tcheck==0) {
- int q1=0,q2=0;
- int ipid1,ipid2=0;
-
- double xtest = randyInstance.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;
- }
- // The new stuff
- 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);
- // End of the new stuff
-
- }
- }
-
- 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(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
-{
- cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
- read();
- cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
- narrowResonanceCrossSection sigma(bbsystem);
- sigma.crossSectionCalculation(_bwnormsave);
- _VMbslope=sigma.slopeParameter();
-}
-
-
-//______________________________________________________________________________
-Gammaanarrowvm::~Gammaanarrowvm()
-{ }
-
-
-//______________________________________________________________________________
-Gammaaincoherentvm::Gammaaincoherentvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
-{
- cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
- read();
- cout<<"Creating and calculating crosssection. Gammaainkoherentvm()"<<endl;
- incoherentVMCrossSection sigma(bbsystem); sigma.crossSectionCalculation(_bwnormsave);
- _VMbslope=sigma.slopeParameter();
-}
-
-
-//______________________________________________________________________________
-Gammaaincoherentvm::~Gammaaincoherentvm()
-{ }
-
-
-//______________________________________________________________________________
-Gammaawidevm::Gammaawidevm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
-{
- cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
- read();
- cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
- wideResonanceCrossSection sigma(bbsystem);
- sigma.crossSectionCalculation(_bwnormsave);
- _VMbslope=sigma.slopeParameter();
-}
-
-
-//______________________________________________________________________________
-Gammaawidevm::~Gammaawidevm()
-{ }
-
-
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff