1 ///////////////////////////////////////////////////////////////////////////
5 // This file is part of starlight.
7 // starlight is free software: you can redistribute it and/or modify
8 // it under the terms of the GNU General Public License as published by
9 // the Free Software Foundation, either version 3 of the License, or
10 // (at your option) any later version.
12 // starlight is distributed in the hope that it will be useful,
13 // but WITHOUT ANY WARRANTY; without even the implied warranty of
14 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
15 // GNU General Public License for more details.
17 // You should have received a copy of the GNU General Public License
18 // along with starlight. If not, see <http://www.gnu.org/licenses/>.
20 ///////////////////////////////////////////////////////////////////////////
22 // File and Version Information:
23 // $Rev:: 179 $: revision of last commit
24 // $Author:: jnystrand $: author of last commit
25 // $Date:: 2014-09-11 20:58:05 +0200 #$: date of last commit
28 // Added incoherent t2-> pt2 selection. Following pp selection scheme
31 ///////////////////////////////////////////////////////////////////////////
40 #include "photonNucleusCrossSection.h"
41 #include "wideResonanceCrossSection.h"
42 #include "narrowResonanceCrossSection.h"
43 #include "incoherentVMCrossSection.h"
48 //______________________________________________________________________________
49 Gammaavectormeson::Gammaavectormeson(beamBeamSystem& bbsystem):eventChannel(bbsystem), _phaseSpaceGen(0) //:readLuminosity(input),_bbs(bbsystem)
51 _VMNPT=inputParametersInstance.nmbPtBinsInterference();
52 _VMWmax=inputParametersInstance.maxW();
53 _VMWmin=inputParametersInstance.minW();
54 _VMYmax=inputParametersInstance.maxRapidity();
56 _VMnumw=inputParametersInstance.nmbWBins();
57 _VMnumy=inputParametersInstance.nmbRapidityBins();
58 _VMgamma_em=inputParametersInstance.beamLorentzGamma();
59 _VMinterferencemode=inputParametersInstance.interferenceEnabled();
60 _VMbslope=0.;//Will define in wide/narrow constructor
61 _VMpidtest=inputParametersInstance.prodParticleType();
62 _VMptmax=inputParametersInstance.maxPtInterference();
63 _VMdpt=inputParametersInstance.ptBinWidthInterference();
64 _VMCoherence=inputParametersInstance.coherentProduction();
65 _VMCoherenceFactor=inputParametersInstance.coherentProduction();
66 _ProductionMode=inputParametersInstance.productionMode();
69 case starlightConstants::RHO:
70 case starlightConstants::RHOZEUS:
74 case starlightConstants::FOURPRONG:
75 // create n-body phase-space generator instance
76 _phaseSpaceGen = new nBodyPhaseSpaceGen();
80 case starlightConstants::OMEGA:
84 case starlightConstants::PHI:
88 case starlightConstants::JPSI:
89 case starlightConstants::JPSI_ee:
90 case starlightConstants::JPSI_mumu:
94 case starlightConstants::JPSI2S:
95 case starlightConstants::JPSI2S_ee:
96 case starlightConstants::JPSI2S_mumu:
100 case starlightConstants::UPSILON:
101 case starlightConstants::UPSILON_ee:
102 case starlightConstants::UPSILON_mumu:
106 case starlightConstants::UPSILON2S:
107 case starlightConstants::UPSILON2S_ee:
108 case starlightConstants::UPSILON2S_mumu:
112 case starlightConstants::UPSILON3S:
113 case starlightConstants::UPSILON3S_ee:
114 case starlightConstants::UPSILON3S_mumu:
118 default: cout<<"No proper vector meson defined, gammaavectormeson::gammaavectormeson"<<endl;
123 //______________________________________________________________________________
124 Gammaavectormeson::~Gammaavectormeson()
127 delete _phaseSpaceGen;
131 //______________________________________________________________________________
132 void Gammaavectormeson::pickwy(double &W, double &Y)
134 double dW, dY, xw,xy,xtest;
137 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
138 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
142 xw = randyInstance.Rndom();// random()/(RAND_MAX+1.0);
143 W = _VMWmin + xw*(_VMWmax-_VMWmin);
145 if (W < 2 * starlightConstants::pionChargedMass)
148 IW = int((W-_VMWmin)/dW); //+ 1;
149 xy = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
150 Y = _VMYmin + xy*(_VMYmax-_VMYmin);
151 IY = int((Y-_VMYmin)/dY); //+ 1;
152 xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
154 if( xtest > _Farray[IW][IY] )
160 //______________________________________________________________________________
161 void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
162 double, // E (unused)
164 double px0, double py0, double pz0,
165 double& px1, double& py1, double& pz1,
166 double& px2, double& py2, double& pz2,
169 // This routine decays a particle into two particles of mass mdec,
170 // taking spin into account
173 // double anglelep[20001],xtest,ytest=0.,dndtheta;
174 double phi,theta,Ecm;
175 double betax,betay,betaz;
177 double E1=0.0,E2=0.0;
179 // set the mass of the daughter particles
180 mdec=getDaughterMass(ipid);
182 // calculate the magnitude of the momenta
184 cout<<" ERROR: W="<<W<<endl;
188 pmag = sqrt(W*W/4. - mdec*mdec);
190 // pick an orientation, based on the spin
191 // phi has a flat distribution in 2*pi
192 phi = randyInstance.Rndom()*2.*starlightConstants::pi;//(random()/(RAND_MAX+1.0))* 2.*pi;
194 // find theta, the angle between one of the outgoing particles and
195 // the beamline, in the frame of the two photons
197 theta=getTheta(ipid);
199 // compute unboosted momenta
200 px1 = sin(theta)*cos(phi)*pmag;
201 py1 = sin(theta)*sin(phi)*pmag;
202 pz1 = cos(theta)*pmag;
207 Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
208 E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
209 E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
215 transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
216 transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
224 //______________________________________________________________________________
225 // decays a particle into four particles with isotropic angular distribution
226 bool Gammaavectormeson::fourBodyDecay
227 (starlightConstants::particleTypeEnum& ipid,
228 const double , // E (unused)
229 const double W, // mass of produced particle
230 const double* p, // momentum of produced particle; expected to have size 3
231 lorentzVector* decayVecs, // array of Lorentz vectors of daughter particles; expected to have size 4
234 const double parentMass = W;
236 // set the mass of the daughter particles
237 const double daughterMass = getDaughterMass(ipid);
238 if (parentMass < 4 * daughterMass){
239 cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
244 // construct parent four-vector
245 const double parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
246 + parentMass * parentMass);
247 const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
249 // setup n-body phase-space generator
250 assert(_phaseSpaceGen);
251 static bool firstCall = true;
253 const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
254 _phaseSpaceGen->setDecay(4, m);
255 // estimate maximum phase-space weight
256 _phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
260 // generate phase-space event
261 if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
264 // set Lorentzvectors of decay daughters
265 for (unsigned int i = 0; i < 4; ++i)
266 decayVecs[i] = _phaseSpaceGen->daughter(i);
271 //______________________________________________________________________________
272 double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
274 //This will return the daughter particles mass, and the final particles outputed id...
275 double ytest=0.,mdec=0.;
278 case starlightConstants::RHO:
279 case starlightConstants::RHOZEUS:
280 case starlightConstants::FOURPRONG:
281 case starlightConstants::OMEGA:
282 mdec = starlightConstants::pionChargedMass;
283 ipid = starlightConstants::PION;
285 case starlightConstants::PHI:
286 mdec = starlightConstants::kaonChargedMass;
287 ipid = starlightConstants::KAONCHARGE;
289 case starlightConstants::JPSI:
290 mdec = starlightConstants::mel;
291 ipid = starlightConstants::ELECTRON;
293 case starlightConstants::JPSI_ee:
294 mdec = starlightConstants::mel;
295 ipid = starlightConstants::ELECTRON;
297 case starlightConstants::JPSI_mumu:
298 mdec = starlightConstants::muonMass;
299 ipid = starlightConstants::MUON;
301 case starlightConstants::JPSI2S_ee:
302 mdec = starlightConstants::mel;
303 ipid = starlightConstants::ELECTRON;
305 case starlightConstants::JPSI2S_mumu:
306 mdec = starlightConstants::muonMass;
307 ipid = starlightConstants::MUON;
310 case starlightConstants::JPSI2S:
311 case starlightConstants::UPSILON:
312 case starlightConstants::UPSILON2S:
313 case starlightConstants::UPSILON3S:
314 // decays 50% to e+/e-, 50% to mu+/mu-
315 ytest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
317 mdec = starlightConstants::muonMass;
318 ipid = starlightConstants::MUON;
320 case starlightConstants::UPSILON_ee:
321 case starlightConstants::UPSILON2S_ee:
322 case starlightConstants::UPSILON3S_ee:
323 mdec = starlightConstants::mel;
324 ipid = starlightConstants::ELECTRON;
326 case starlightConstants::UPSILON_mumu:
327 case starlightConstants::UPSILON2S_mumu:
328 case starlightConstants::UPSILON3S_mumu:
329 mdec = starlightConstants::muonMass;
330 ipid = starlightConstants::MUON;
332 default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
339 //______________________________________________________________________________
340 double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid)
342 //This depends on the decay angular distribution
343 //Valid for rho, phi, omega.
349 theta = starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
350 xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
351 // Follow distribution for helicity +/-1
352 // Eq. 19 of J. Breitweg et al., Eur. Phys. J. C2, 247 (1998)
357 case starlightConstants::MUON:
358 case starlightConstants::ELECTRON:
359 //primarily for upsilon/j/psi. VM->ee/mumu
360 dndtheta = sin(theta)*(1.+((cos(theta))*(cos(theta))));
363 case starlightConstants::PION:
364 case starlightConstants::KAONCHARGE:
366 dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
369 default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
380 //______________________________________________________________________________
381 double Gammaavectormeson::getWidth()
387 //______________________________________________________________________________
388 double Gammaavectormeson::getMass()
394 //______________________________________________________________________________
395 double Gammaavectormeson::getSpin()
397 return 1.0; //VM spins are the same
401 //______________________________________________________________________________
402 void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)
404 // This subroutine calculates momentum and energy of vector meson
405 // given W and Y, without interference. Subroutine vmpt.f handles
406 // production with interference
409 double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
410 double px1,py1,px2,py2;
411 double pt,xt,xtest,ytest;
412 // double photon_spectrum;
415 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
416 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
418 //Find Egam,Epom in CM frame
419 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
420 if( _ProductionMode == 2 ){
422 Epom = 0.5*W*exp(-Y);
424 Egam = 0.5*W*exp(-Y);
427 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
428 if( _ProductionMode == 2 ){
429 Egam = 0.5*W*exp(-Y);
433 Epom = 0.5*W*exp(-Y);
437 Epom = 0.5*W*exp(-Y);
440 // cout<<" Y: "<<Y<<" W: "<<W<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl;
442 phi1 = 2.*starlightConstants::pi*randyInstance.Rndom();
444 // cout<<" pt1: "<<pt1<<endl;
446 if( (_bbs.beam1().A()==1 && _bbs.beam2().A()==1) ||
447 (_ProductionMode == 4) ) {
448 if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
449 // Use dipole form factor for light VM
451 xtest = 2.0*randyInstance.Rndom();
452 double ttest = xtest*xtest;
453 ytest = randyInstance.Rndom();
455 double yprob = xtest*_bbs.beam1().dipoleFormFactor(ttest,t0)*_bbs.beam1().dipoleFormFactor(ttest,t0);
456 if( ytest > yprob ) goto L555vm;
460 //Use dsig/dt= exp(-_VMbslope*t) for heavy VM
461 xtest = randyInstance.Rndom();
462 t2 = (-1./_VMbslope)*log(xtest);
467 tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
470 cout<<" WARNING: tmin= "<<tmin<<endl;
471 cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl;
472 cout<<" Will pick a new W,Y "<<endl;
477 xt = randyInstance.Rndom();
478 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
479 if( _ProductionMode == 2 ){
480 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
482 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
484 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
485 if( _ProductionMode == 2 ){
486 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
488 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
491 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
494 xtest = randyInstance.Rndom();
497 if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
498 if(1.0 < _bbs.beam2().formFactor(t2)*pt2) cout <<"POMERON"<<endl;
499 if( xtest > _bbs.beam2().formFactor(t2)*pt2) goto L203vm;
504 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
505 if( _ProductionMode == 2 ){
506 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
508 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
510 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
511 if( _ProductionMode == 2 ){
512 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
514 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
517 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
520 if( xtest > comp ) goto L203vm;
524 if(_VMCoherence==0 && (!(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2))){
525 //dsig/dt= exp(-_VMbslope*t)
526 xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
527 t2 = (-1./_VMbslope)*log(xtest);
533 phi2 = 2.*starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
535 // cout<<" pt2: "<<pt2<<endl;
542 // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
545 pt = sqrt( px*px + py*py );
547 E = sqrt(W*W+pt*pt)*cosh(Y);
548 pz = sqrt(W*W+pt*pt)*sinh(Y);
550 // cout<< " Y = "<<Y<<" W = "<<W<<" Egam = "<<Egam<<" gamma = "<<_VMgamma_em<<endl;
552 // Randomly choose to make pz negative 50% of the time
553 if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
555 }else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
559 if (randyInstance.Rndom() >= 0.5) pz = -pz;
564 //______________________________________________________________________________
565 double Gammaavectormeson::pTgamma(double E)
567 // returns on random draw from pp(E) distribution
568 double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
569 double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
572 ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
573 //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
574 Cm = sqrt(3.)*E/_VMgamma_em;
575 // If E is very small, the drawing of a pT below is extremely slow.
576 // ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E.
577 // Should have no observable consequences (JN, SRK 11-Sep-2014)
578 if( E < 0.0005 )return Cm;
580 //the amplitude of the p_t spectrum at the maximum
582 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
583 if( _ProductionMode == 2 ){
584 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
586 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
588 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
589 if( _ProductionMode == 2 ){
590 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
592 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
595 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
598 Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
600 //pick a test value pp, and find the amplitude there
601 x = randyInstance.Rndom();
603 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
604 if( _ProductionMode == 2 ){
605 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
606 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
608 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
609 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
611 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
612 if( _ProductionMode == 2 ){
613 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
614 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
616 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
617 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
620 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
621 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
624 test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
627 u = randyInstance.Rndom();
633 x =randyInstance.Rndom();
634 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
635 if( _ProductionMode == 2 ){
636 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
637 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
639 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
640 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
642 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
643 if( _ProductionMode == 2 ){
644 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
645 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
647 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
648 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
651 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
652 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
654 test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
662 //______________________________________________________________________________
663 void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
664 int&) // tcheck (unused)
666 // This function calculates momentum and energy of vector meson
667 // given W and Y, including interference.
668 // It gets the pt distribution from a lookup table.
669 double dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
672 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
673 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
675 // Y is already fixed; choose a pt
676 // Follow the approavh in pickwy.f
677 // in _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
679 IY=int(fabs(Y)/dY);//+1;
680 if (IY > (_VMnumy/2)-1){
684 yleft=fabs(Y)-(IY)*dY;
687 xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
689 for(j=0;j<_VMNPT+1;j++){
690 if (xpt < _fptarray[IY][j]) goto L60;
694 // now do linear interpolation - start with extremes
697 pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
701 pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
705 // we're in the middle
707 ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
708 pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
710 // at an extreme in y?
711 if (IY == (_VMnumy/2)-1){
716 // interpolate in y repeat for next fractional y bin
718 for(j=0;j<_VMNPT+1;j++){
719 if (xpt < _fptarray[IY+1][j]) goto L90;
723 // now do linear interpolation - start with extremes
726 pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
730 pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
734 // we're in the middle
736 ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
737 pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
741 // now interpolate in y
743 pt=yfract*pt2+(1-yfract)*pt1;
749 theta=2.*starlightConstants::pi*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
753 // I guess W is the mass of the vector meson (not necessarily
754 // on-mass-shell), and E is the energy
756 E = sqrt(W*W+pt*pt)*cosh(Y);
757 pz = sqrt(W*W+pt*pt)*sinh(Y);
758 // randomly choose to make pz negative 50% of the time
759 if(randyInstance.Rndom()>=0.5) pz = -pz;
763 //______________________________________________________________________________
764 starlightConstants::event Gammaavectormeson::produceEvent(int&)
766 // Not used; return default event
767 return starlightConstants::event();
771 //______________________________________________________________________________
772 upcEvent Gammaavectormeson::produceEvent()
774 // The new event type
779 starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
780 starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
782 if (_VMpidtest == starlightConstants::FOURPRONG) {
783 double comenergy = 0;
784 double mom[3] = {0, 0, 0};
786 lorentzVector decayVecs[4];
789 pickwy(comenergy, rapidity);
790 if (_VMinterferencemode == 0)
791 momenta(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
792 else if (_VMinterferencemode==1)
793 vmpt(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
794 } while (!fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent));
795 if ((iFbadevent == 0) and (tcheck == 0))
796 for (unsigned int i = 0; i < 4; ++i) {
797 starlightParticle daughter(decayVecs[i].GetPx(),
798 decayVecs[i].GetPy(),
799 decayVecs[i].GetPz(),
800 starlightConstants::UNKNOWN, // energy
801 starlightConstants::UNKNOWN, // _mass
804 event.addParticle(daughter);
807 double comenergy = 0.;
808 double rapidity = 0.;
810 double momx=0.,momy=0.,momz=0.;
812 double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
813 bool accepted = false;
816 pickwy(comenergy,rapidity);
818 if (_VMinterferencemode==0){
819 momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
820 } else if (_VMinterferencemode==1){
821 vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck);
824 // cout << "_ptCutMin: " << _ptCutMin << " _ptCutMax: " << _ptCutMax << " _etaCutMin: " << _etaCutMin << " _etaCutMax: " << _etaCutMax << endl;
826 //cout << "n tries: " << _nmbAttempts<< endl;
828 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
829 double pt1chk = sqrt(px1*px1+py1*py1);
830 double pt2chk = sqrt(px2*px2+py2*py2);
832 //cout << "pt1: " << pt1chk << " pt2: " << pt2chk << endl;
833 double eta1 = pseudoRapidity(px1, py1, pz1);
834 double eta2 = pseudoRapidity(px2, py2, pz2);
835 //cout << "eta1: " << eta1 << " eta2: " << eta2 << endl;
836 if(_ptCutEnabled && !_etaCutEnabled){
837 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
842 else if(!_ptCutEnabled && _etaCutEnabled){
843 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
848 else if(_ptCutEnabled && _etaCutEnabled){
849 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
850 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
856 else if(!_ptCutEnabled && !_etaCutEnabled)
859 }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
861 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
863 if (iFbadevent==0&&tcheck==0) {
867 double xtest = randyInstance.Rndom();
878 if ( ipid == 11 || ipid == 13 ){
886 double md = getDaughterMass(vmpid);
887 double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1);
888 starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
889 event.addParticle(particle1);
891 double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2);
892 starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
893 event.addParticle(particle2);
894 // End of the new stuff
902 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
904 double pT = sqrt(px*px + py*py);
905 double p = sqrt(pz*pz + pT*pT);
906 double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
910 //______________________________________________________________________________
911 Gammaanarrowvm::Gammaanarrowvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
913 cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
915 cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
916 narrowResonanceCrossSection sigma(bbsystem);
917 sigma.crossSectionCalculation(_bwnormsave);
918 _VMbslope=sigma.slopeParameter();
922 //______________________________________________________________________________
923 Gammaanarrowvm::~Gammaanarrowvm()
927 //______________________________________________________________________________
928 Gammaaincoherentvm::Gammaaincoherentvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
930 cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
932 cout<<"Creating and calculating crosssection. Gammaainkoherentvm()"<<endl;
933 incoherentVMCrossSection sigma(bbsystem); sigma.crossSectionCalculation(_bwnormsave);
934 _VMbslope=sigma.slopeParameter();
938 //______________________________________________________________________________
939 Gammaaincoherentvm::~Gammaaincoherentvm()
943 //______________________________________________________________________________
944 Gammaawidevm::Gammaawidevm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
946 cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
948 cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
949 wideResonanceCrossSection sigma(bbsystem);
950 sigma.crossSectionCalculation(_bwnormsave);
951 _VMbslope=sigma.slopeParameter();
955 //______________________________________________________________________________
956 Gammaawidevm::~Gammaawidevm()