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:: 164 $: revision of last commit
24 // $Author:: odjuvsla $: author of last commit
25 // $Date:: 2013-10-06 16:18:08 +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;
576 //the amplitude of the p_t spectrum at the maximum
578 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
579 if( _ProductionMode == 2 ){
580 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
582 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
584 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
585 if( _ProductionMode == 2 ){
586 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
588 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
591 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
594 Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
596 //pick a test value pp, and find the amplitude there
597 x = randyInstance.Rndom();
599 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
600 if( _ProductionMode == 2 ){
601 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
602 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
604 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
605 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
607 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
608 if( _ProductionMode == 2 ){
609 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
610 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
612 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
613 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
616 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
617 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
620 test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
623 u = randyInstance.Rndom();
629 x =randyInstance.Rndom();
630 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
631 if( _ProductionMode == 2 ){
632 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
633 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
635 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
636 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
638 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
639 if( _ProductionMode == 2 ){
640 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
641 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
643 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
644 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
647 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
648 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
650 test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
658 //______________________________________________________________________________
659 void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
660 int&) // tcheck (unused)
662 // This function calculates momentum and energy of vector meson
663 // given W and Y, including interference.
664 // It gets the pt distribution from a lookup table.
665 double dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
668 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
669 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
671 // Y is already fixed; choose a pt
672 // Follow the approavh in pickwy.f
673 // in _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
675 IY=int(fabs(Y)/dY);//+1;
676 if (IY > (_VMnumy/2)-1){
680 yleft=fabs(Y)-(IY)*dY;
683 xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
685 for(j=0;j<_VMNPT+1;j++){
686 if (xpt < _fptarray[IY][j]) goto L60;
690 // now do linear interpolation - start with extremes
693 pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
697 pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
701 // we're in the middle
703 ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
704 pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
706 // at an extreme in y?
707 if (IY == (_VMnumy/2)-1){
712 // interpolate in y repeat for next fractional y bin
714 for(j=0;j<_VMNPT+1;j++){
715 if (xpt < _fptarray[IY+1][j]) goto L90;
719 // now do linear interpolation - start with extremes
722 pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
726 pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
730 // we're in the middle
732 ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
733 pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
737 // now interpolate in y
739 pt=yfract*pt2+(1-yfract)*pt1;
745 theta=2.*starlightConstants::pi*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
749 // I guess W is the mass of the vector meson (not necessarily
750 // on-mass-shell), and E is the energy
752 E = sqrt(W*W+pt*pt)*cosh(Y);
753 pz = sqrt(W*W+pt*pt)*sinh(Y);
754 // randomly choose to make pz negative 50% of the time
755 if(randyInstance.Rndom()>=0.5) pz = -pz;
759 //______________________________________________________________________________
760 starlightConstants::event Gammaavectormeson::produceEvent(int&)
762 // Not used; return default event
763 return starlightConstants::event();
767 //______________________________________________________________________________
768 upcEvent Gammaavectormeson::produceEvent()
770 // The new event type
775 starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
776 starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
778 if (_VMpidtest == starlightConstants::FOURPRONG) {
779 double comenergy = 0;
780 double mom[3] = {0, 0, 0};
782 lorentzVector decayVecs[4];
785 pickwy(comenergy, rapidity);
786 if (_VMinterferencemode == 0)
787 momenta(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
788 else if (_VMinterferencemode==1)
789 vmpt(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
790 } while (!fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent));
791 if ((iFbadevent == 0) and (tcheck == 0))
792 for (unsigned int i = 0; i < 4; ++i) {
793 starlightParticle daughter(decayVecs[i].GetPx(),
794 decayVecs[i].GetPy(),
795 decayVecs[i].GetPz(),
796 starlightConstants::UNKNOWN, // energy
797 starlightConstants::UNKNOWN, // _mass
800 event.addParticle(daughter);
803 double comenergy = 0.;
804 double rapidity = 0.;
806 double momx=0.,momy=0.,momz=0.;
808 double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
809 bool accepted = false;
812 pickwy(comenergy,rapidity);
814 if (_VMinterferencemode==0){
815 momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
816 } else if (_VMinterferencemode==1){
817 vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck);
820 // cout << "_ptCutMin: " << _ptCutMin << " _ptCutMax: " << _ptCutMax << " _etaCutMin: " << _etaCutMin << " _etaCutMax: " << _etaCutMax << endl;
822 //cout << "n tries: " << _nmbAttempts<< endl;
824 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
825 double pt1chk = sqrt(px1*px1+py1*py1);
826 double pt2chk = sqrt(px2*px2+py2*py2);
828 //cout << "pt1: " << pt1chk << " pt2: " << pt2chk << endl;
829 double eta1 = pseudoRapidity(px1, py1, pz1);
830 double eta2 = pseudoRapidity(px2, py2, pz2);
831 //cout << "eta1: " << eta1 << " eta2: " << eta2 << endl;
832 if(_ptCutEnabled && !_etaCutEnabled){
833 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
838 else if(!_ptCutEnabled && _etaCutEnabled){
839 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
844 else if(_ptCutEnabled && _etaCutEnabled){
845 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
846 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
852 else if(!_ptCutEnabled && !_etaCutEnabled)
855 }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
857 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
859 if (iFbadevent==0&&tcheck==0) {
863 double xtest = randyInstance.Rndom();
874 if ( ipid == 11 || ipid == 13 ){
882 double md = getDaughterMass(vmpid);
883 double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1);
884 starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
885 event.addParticle(particle1);
887 double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2);
888 starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
889 event.addParticle(particle2);
890 // End of the new stuff
898 double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
900 double pT = sqrt(px*px + py*py);
901 double p = sqrt(pz*pz + pT*pT);
902 double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
906 //______________________________________________________________________________
907 Gammaanarrowvm::Gammaanarrowvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
909 cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
911 cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
912 narrowResonanceCrossSection sigma(bbsystem);
913 sigma.crossSectionCalculation(_bwnormsave);
914 _VMbslope=sigma.slopeParameter();
918 //______________________________________________________________________________
919 Gammaanarrowvm::~Gammaanarrowvm()
923 //______________________________________________________________________________
924 Gammaaincoherentvm::Gammaaincoherentvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
926 cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
928 cout<<"Creating and calculating crosssection. Gammaainkoherentvm()"<<endl;
929 incoherentVMCrossSection sigma(bbsystem); sigma.crossSectionCalculation(_bwnormsave);
930 _VMbslope=sigma.slopeParameter();
934 //______________________________________________________________________________
935 Gammaaincoherentvm::~Gammaaincoherentvm()
939 //______________________________________________________________________________
940 Gammaawidevm::Gammaawidevm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
942 cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
944 cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
945 wideResonanceCrossSection sigma(bbsystem);
946 sigma.crossSectionCalculation(_bwnormsave);
947 _VMbslope=sigma.slopeParameter();
951 //______________________________________________________________________________
952 Gammaawidevm::~Gammaawidevm()