]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/.svn/text-base/gammaavm.cpp.svn-base
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
STARLIGHT code and interface
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / .svn / text-base / gammaavm.cpp.svn-base
CommitLineData
da32329d
AM		 1///////////////////////////////////////////////////////////////////////////
2//
3// Copyright 2010
4//
5// This file is part of starlight.
6//
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.
11//
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.
16//
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/>.
19//
20///////////////////////////////////////////////////////////////////////////
21//
22// File and Version Information:
23// $Rev:: $: revision of last commit
24// $Author:: $: author of last commit
25// $Date:: $: date of last commit
26//
27// Description:
28// Added incoherent t2-> pt2 selection. Following pp selection scheme
29//
30//
31///////////////////////////////////////////////////////////////////////////
32
33
34#include <iostream>
35#include <fstream>
36#include <cassert>
37#include <cmath>
38
39#include "gammaavm.h"
40#include "photonNucleusCrossSection.h"
41#include "wideResonanceCrossSection.h"
42#include "narrowResonanceCrossSection.h"
43#include "incoherentVMCrossSection.h"
44
45using namespace std;
46
47
48//______________________________________________________________________________
49Gammaavectormeson::Gammaavectormeson(beamBeamSystem& bbsystem):eventChannel(bbsystem), _phaseSpaceGen(0) //:readLuminosity(input),_bbs(bbsystem)
50{
51 _VMNPT=inputParametersInstance.nmbPtBinsInterference();
52 _VMWmax=inputParametersInstance.maxW();
53 _VMWmin=inputParametersInstance.minW();
54 _VMYmax=inputParametersInstance.maxRapidity();
55 _VMYmin=-1.*_VMYmax;
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();
67
68 switch(_VMpidtest){
69 case starlightConstants::RHO:
70 case starlightConstants::RHOZEUS:
71 _width=0.1507;
72 _mass=0.7685;
73 break;
74 case starlightConstants::FOURPRONG:
75 // create n-body phase-space generator instance
76 _phaseSpaceGen = new nBodyPhaseSpaceGen();
77 _width = 0.360;
78 _mass = 1.350;
79 break;
80 case starlightConstants::OMEGA:
81 _width=0.00843;
82 _mass=0.78194;
83 break;
84 case starlightConstants::PHI:
85 _width=0.00443;
86 _mass=1.019413;
87 break;
88 case starlightConstants::JPSI:
89 case starlightConstants::JPSI_ee:
90 case starlightConstants::JPSI_mumu:
91 _width=0.000091;
92 _mass=3.09692;
93 break;
94 case starlightConstants::JPSI2S:
95 case starlightConstants::JPSI2S_ee:
96 case starlightConstants::JPSI2S_mumu:
97 _width=0.000337;
98 _mass=3.686093;
99 break;
100 case starlightConstants::UPSILON:
101 case starlightConstants::UPSILON_ee:
102 case starlightConstants::UPSILON_mumu:
103 _width=0.00005402;
104 _mass=9.46030;
105 break;
106 case starlightConstants::UPSILON2S:
107 case starlightConstants::UPSILON2S_ee:
108 case starlightConstants::UPSILON2S_mumu:
109 _width=0.00003198;
110 _mass=10.02326;
111 break;
112 case starlightConstants::UPSILON3S:
113 case starlightConstants::UPSILON3S_ee:
114 case starlightConstants::UPSILON3S_mumu:
115 _width=0.00002032;
116 _mass=10.3552;
117 break;
118 default: cout<<"No proper vector meson defined, gammaavectormeson::gammaavectormeson"<<endl;
119 }
120}
121
122
123//______________________________________________________________________________
124Gammaavectormeson::~Gammaavectormeson()
125{
126 if (_phaseSpaceGen)
127 delete _phaseSpaceGen;
128}
129
130
131//______________________________________________________________________________
132void Gammaavectormeson::pickwy(double &W, double &Y)
133{
134 double dW, dY, xw,xy,xtest;
135 int IW,IY;
136
137 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
138 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
139
140 L201pwy:
141
142 xw = randyInstance.Rndom();// random()/(RAND_MAX+1.0);
143 W = _VMWmin + xw*(_VMWmax-_VMWmin);
144
145 if (W < 2 * starlightConstants::pionChargedMass)
146 goto L201pwy;
147
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);
153
154 if( xtest > _Farray[IW][IY] )
155 goto L201pwy;
156
157}
158
159
160//______________________________________________________________________________
161void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
162 double, // E (unused)
163 double W,
164 double px0, double py0, double pz0,
165 double& px1, double& py1, double& pz1,
166 double& px2, double& py2, double& pz2,
167 int& iFbadevent)
168{
169 // This routine decays a particle into two particles of mass mdec,
170 // taking spin into account
171
172 double pmag;
173 // double anglelep[20001],xtest,ytest=0.,dndtheta;
174 double phi,theta,Ecm;
175 double betax,betay,betaz;
176 double mdec=0.0;
177 double E1=0.0,E2=0.0;
178
179 // set the mass of the daughter particles
180 mdec=getDaughterMass(ipid);
181
182 // calculate the magnitude of the momenta
183 if(W < 2*mdec){
184 cout<<" ERROR: W="<<W<<endl;
185 iFbadevent = 1;
186 return;
187 }
188 pmag = sqrt(W*W/4. - mdec*mdec);
189
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;
193
194 // find theta, the angle between one of the outgoing particles and
195 // the beamline, in the frame of the two photons
196
197 theta=getTheta(ipid);
198
199 // compute unboosted momenta
200 px1 = sin(theta)*cos(phi)*pmag;
201 py1 = sin(theta)*sin(phi)*pmag;
202 pz1 = cos(theta)*pmag;
203 px2 = -px1;
204 py2 = -py1;
205 pz2 = -pz1;
206
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);
210
211 betax = -(px0/Ecm);
212 betay = -(py0/Ecm);
213 betaz = -(pz0/Ecm);
214
215 transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
216 transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
217
218 if(iFbadevent == 1)
219 return;
220
221}
222
223
224//______________________________________________________________________________
225// decays a particle into four particles with isotropic angular distribution
226bool 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
232 int& iFbadevent)
233{
234 const double parentMass = W;
235
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;
240 iFbadevent = 1;
241 return false;
242 }
243
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);
248
249 // setup n-body phase-space generator
250 assert(_phaseSpaceGen);
251 static bool firstCall = true;
252 if (firstCall) {
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));
257 firstCall = false;
258 }
259
260 // generate phase-space event
261 if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
262 return false;
263
264 // set Lorentzvectors of decay daughters
265 for (unsigned int i = 0; i < 4; ++i)
266 decayVecs[i] = _phaseSpaceGen->daughter(i);
267 return true;
268}
269
270
271//______________________________________________________________________________
272double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
273{
274 //This will return the daughter particles mass, and the final particles outputed id...
275 double ytest=0.,mdec=0.;
276
277 switch(_VMpidtest){
278 case starlightConstants::RHO:
279 case starlightConstants::RHOZEUS:
280 case starlightConstants::FOURPRONG:
281 case starlightConstants::OMEGA:
282 mdec = starlightConstants::pionChargedMass;
283 ipid = starlightConstants::PION;
284 break;
285 case starlightConstants::PHI:
286 mdec = starlightConstants::kaonChargedMass;
287 ipid = starlightConstants::KAONCHARGE;
288 break;
289 case starlightConstants::JPSI:
290 mdec = starlightConstants::mel;
291 ipid = starlightConstants::ELECTRON;
292 break;
293 case starlightConstants::JPSI_ee:
294 mdec = starlightConstants::mel;
295 ipid = starlightConstants::ELECTRON;
296 break;
297 case starlightConstants::JPSI_mumu:
298 mdec = starlightConstants::muonMass;
299 ipid = starlightConstants::MUON;
300 break;
301 case starlightConstants::JPSI2S_ee:
302 mdec = starlightConstants::mel;
303 ipid = starlightConstants::ELECTRON;
304 break;
305 case starlightConstants::JPSI2S_mumu:
306 mdec = starlightConstants::muonMass;
307 ipid = starlightConstants::MUON;
308 break;
309
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);
316
317 mdec = starlightConstants::muonMass;
318 ipid = starlightConstants::MUON;
319 break;
320 case starlightConstants::UPSILON_ee:
321 case starlightConstants::UPSILON2S_ee:
322 case starlightConstants::UPSILON3S_ee:
323 mdec = starlightConstants::mel;
324 ipid = starlightConstants::ELECTRON;
325 break;
326 case starlightConstants::UPSILON_mumu:
327 case starlightConstants::UPSILON2S_mumu:
328 case starlightConstants::UPSILON3S_mumu:
329 mdec = starlightConstants::muonMass;
330 ipid = starlightConstants::MUON;
331 break;
332 default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
333 }
334
335 return mdec;
336}
337
338
339//______________________________________________________________________________
340double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid)
341{
342 //This depends on the decay angular distribution
343 //Valid for rho, phi, omega.
344 double theta=0.;
345 double xtest=0.;
346 double dndtheta=0.;
347 L200td:
348
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)
353 // SRK 11/14/2000
354
355 switch(ipid){
356
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))));
361 break;
362
363 case starlightConstants::PION:
364 case starlightConstants::KAONCHARGE:
365 //rhos etc
366 dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
367 break;
368
369 default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
370 }//end of switch
371
372 if(xtest > dndtheta)
373 goto L200td;
374
375 return theta;
376
377}
378
379
380//______________________________________________________________________________
381double Gammaavectormeson::getWidth()
382{
383 return _width;
384}
385
386
387//______________________________________________________________________________
388double Gammaavectormeson::getMass()
389{
390 return _mass;
391}
392
393
394//______________________________________________________________________________
395double Gammaavectormeson::getSpin()
396{
397 return 1.0; //VM spins are the same
398}
399
400
401//______________________________________________________________________________
402void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)
403{
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
407
408 double dW,dY;
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;
413 double t2;
414
415 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
416 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
417
418 //Find Egam,Epom in CM frame
419 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
420 if( _ProductionMode == 2 ){
421 Egam = 0.5*W*exp(Y);
422 Epom = 0.5*W*exp(-Y);
423 }else{
424 Egam = 0.5*W*exp(-Y);
425 Epom = 0.5*W*exp(Y);
426 }
427 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
428 if( _ProductionMode == 2 ){
429 Egam = 0.5*W*exp(-Y);
430 Epom = 0.5*W*exp(Y);
431 }else{
432 Egam = 0.5*W*exp(Y);
433 Epom = 0.5*W*exp(-Y);
434 }
435 } else {
436 Egam = 0.5*W*exp(Y);
437 Epom = 0.5*W*exp(-Y);
438 }
439
440 // cout<<" Y: "<<Y<<" W: "<<W<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl;
441 pt1 = pTgamma(Egam);
442 phi1 = 2.*starlightConstants::pi*randyInstance.Rndom();
443
444 // cout<<" pt1: "<<pt1<<endl;
445
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
450 L555vm:
451 xtest = 2.0*randyInstance.Rndom();
452 double ttest = xtest*xtest;
453 ytest = randyInstance.Rndom();
454 double t0 = 1./2.23;
455 double yprob = xtest*_bbs.beam1().dipoleFormFactor(ttest,t0)*_bbs.beam1().dipoleFormFactor(ttest,t0);
456 if( ytest > yprob ) goto L555vm;
457 t2 = ttest;
458 pt2 = xtest;
459 }else{
460 //Use dsig/dt= exp(-_VMbslope*t) for heavy VM
461 xtest = randyInstance.Rndom();
462 t2 = (-1./_VMbslope)*log(xtest);
463 pt2 = sqrt(1.*t2);
464 }
465 } else {
466 // >> Check tmin
467 tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
468
469 if(tmin > 0.5){
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;
473 tcheck = 1;
474 return;
475 }
476 L203vm:
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();
481 }else{
482 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
483 }
484 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
485 if( _ProductionMode == 2 ){
486 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
487 }else{
488 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
489 }
490 } else {
491 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
492 }
493
494 xtest = randyInstance.Rndom();
495 t2 = tmin + pt2*pt2;
496
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;
500 }
501 else{
502
503 double comp=0.0;
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;
507 }else{
508 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
509 }
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;
513 }else{
514 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
515 }
516 } else {
517 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
518 }
519
520 if( xtest > comp ) goto L203vm;
521 }
522
523 /*
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);
528 pt2 = sqrt(1.*t2);
529 }
530 */
531
532 }//else end from pp
533 phi2 = 2.*starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
534
535 // cout<<" pt2: "<<pt2<<endl;
536
537 px1 = pt1*cos(phi1);
538 py1 = pt1*sin(phi1);
539 px2 = pt2*cos(phi2);
540 py2 = pt2*sin(phi2);
541
542 // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
543 px = px1 + px2;
544 py = py1 + py2;
545 pt = sqrt( px*px + py*py );
546
547 E = sqrt(W*W+pt*pt)*cosh(Y);
548 pz = sqrt(W*W+pt*pt)*sinh(Y);
549
550 // cout<< " Y = "<<Y<<" W = "<<W<<" Egam = "<<Egam<<" gamma = "<<_VMgamma_em<<endl;
551
552 // Randomly choose to make pz negative 50% of the time
553 if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
554 pz = -pz;
555 }else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
556 // Don't switch
557 }
558 else{
559 if (randyInstance.Rndom() >= 0.5) pz = -pz;
560 }
561
562}
563
564//______________________________________________________________________________
565double Gammaavectormeson::pTgamma(double E)
566{
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.;
570 int satisfy =0;
571
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
576 //the amplitude of the p_t spectrum at the maximum
577
578 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
579 if( _ProductionMode == 2 ){
580 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
581 }else{
582 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
583 }
584 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
585 if( _ProductionMode == 2 ){
586 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
587 }else{
588 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
589 }
590 } else {
591 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
592 }
593
594 Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
595
596 //pick a test value pp, and find the amplitude there
597 x = randyInstance.Rndom();
598
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);
603 }else{
604 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
605 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
606 }
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);
611 }else{
612 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
613 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
614 }
615 } else {
616 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
617 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
618 }
619
620 test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
621
622 while(satisfy==0){
623 u = randyInstance.Rndom();
624 if(u*Coef <= test)
625 {
626 satisfy =1;
627 }
628 else{
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);
634 }else{
635 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
636 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
637 }
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);
642 }else{
643 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
644 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
645 }
646 } else {
647 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
648 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
649 }
650 test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
651 }
652 }
653
654 return pp;
655}
656
657
658//______________________________________________________________________________
659void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
660 int&) // tcheck (unused)
661{
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.;
666 int IY=0,j=0;
667
668 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
669 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
670
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
674
675 IY=int(fabs(Y)/dY);//+1;
676 if (IY > (_VMnumy/2)-1){
677 IY=(_VMnumy/2)-1;
678 }
679
680 yleft=fabs(Y)-(IY)*dY;
681 yfract=yleft*dY;
682
683 xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
684
685 for(j=0;j<_VMNPT+1;j++){
686 if (xpt < _fptarray[IY][j]) goto L60;
687 }
688 L60:
689
690 // now do linear interpolation - start with extremes
691
692 if (j == 0){
693 pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
694 goto L80;
695 }
696 if (j == _VMNPT){
697 pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
698 goto L80;
699 }
700
701 // we're in the middle
702
703 ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
704 pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
705
706 // at an extreme in y?
707 if (IY == (_VMnumy/2)-1){
708 pt=pt1;
709 goto L120;
710 }
711 L80:
712 // interpolate in y repeat for next fractional y bin
713
714 for(j=0;j<_VMNPT+1;j++){
715 if (xpt < _fptarray[IY+1][j]) goto L90;
716 }
717 L90:
718
719 // now do linear interpolation - start with extremes
720
721 if (j == 0){
722 pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
723 goto L100;
724 }
725 if (j == _VMNPT){
726 pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
727 goto L100;
728 }
729
730 // we're in the middle
731
732 ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
733 pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
734
735 L100:
736
737 // now interpolate in y
738
739 pt=yfract*pt2+(1-yfract)*pt1;
740
741 L120:
742
743 // we have a pt
744
745 theta=2.*starlightConstants::pi*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
746 px=pt*cos(theta);
747 py=pt*sin(theta);
748
749 // I guess W is the mass of the vector meson (not necessarily
750 // on-mass-shell), and E is the energy
751
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;
756}
757
758
759//______________________________________________________________________________
760starlightConstants::event Gammaavectormeson::produceEvent(int&)
761{
762 // Not used; return default event
763 return starlightConstants::event();
764}
765
766
767//______________________________________________________________________________
768upcEvent Gammaavectormeson::produceEvent()
769{
770 // The new event type
771 upcEvent event;
772
773 int iFbadevent=0;
774 int tcheck=0;
775 starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
776 starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
777
778 if (_VMpidtest == starlightConstants::FOURPRONG) {
779 double comenergy = 0;
780 double mom[3] = {0, 0, 0};
781 double E = 0;
782 lorentzVector decayVecs[4];
783 do {
784 double rapidity = 0;
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
798 ipid,
799 (i < 2) ? -1 : +1);
800 event.addParticle(daughter);
801 }
802 } else {
803 double comenergy = 0.;
804 double rapidity = 0.;
805 double E = 0.;
806 double momx=0.,momy=0.,momz=0.;
807
808 double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
809 bool accepted = false;
810 // if(_accCut){
811 do{
812 pickwy(comenergy,rapidity);
813
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);
818 }
819
820 // cout << "_ptCutMin: " << _ptCutMin << " _ptCutMax: " << _ptCutMax << " _etaCutMin: " << _etaCutMin << " _etaCutMax: " << _etaCutMax << endl;
821 _nmbAttempts++;
822 //cout << "n tries: " << _nmbAttempts<< endl;
823 vmpid = ipid;
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);
827
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){
834 accepted = true;
835 _nmbAccepted++;
836 }
837 }
838 else if(!_ptCutEnabled && _etaCutEnabled){
839 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
840 accepted = true;
841 _nmbAccepted++;
842 }
843 }
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){
847 accepted = true;
848 _nmbAccepted++;
849 }
850 }
851 }
852 else if(!_ptCutEnabled && !_etaCutEnabled)
853 _nmbAccepted++;
854
855 }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
856 /* }else{
857 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
858 }*/
859 if (iFbadevent==0&&tcheck==0) {
860 int q1=0,q2=0;
861 int ipid1,ipid2=0;
862
863 double xtest = randyInstance.Rndom();
864 if (xtest<0.5)
865 {
866 q1=1;
867 q2=-1;
868 }
869 else {
870 q1=-1;
871 q2=1;
872 }
873
874 if ( ipid == 11 || ipid == 13 ){
875 ipid1 = -q1*ipid;
876 ipid2 = -q2*ipid;
877 } else {
878 ipid1 = q1*ipid;
879 ipid2 = q2*ipid;
880 }
881 // The new stuff
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);
886
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
891
892 }
893 }
894
895 return event;
896
897}
898double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
899{
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));}
903 return eta;
904}
905
906//______________________________________________________________________________
907Gammaanarrowvm::Gammaanarrowvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
908{
909 cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
910 read();
911 cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
912 narrowResonanceCrossSection sigma(bbsystem);
913 sigma.crossSectionCalculation(_bwnormsave);
914 _VMbslope=sigma.slopeParameter();
915}
916
917
918//______________________________________________________________________________
919Gammaanarrowvm::~Gammaanarrowvm()
920{ }
921
922
923//______________________________________________________________________________
924Gammaaincoherentvm::Gammaaincoherentvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
925{
926 cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
927 read();
928 cout<<"Creating and calculating crosssection. Gammaainkoherentvm()"<<endl;
929 incoherentVMCrossSection sigma(bbsystem); sigma.crossSectionCalculation(_bwnormsave);
930 _VMbslope=sigma.slopeParameter();
931}
932
933
934//______________________________________________________________________________
935Gammaaincoherentvm::~Gammaaincoherentvm()
936{ }
937
938
939//______________________________________________________________________________
940Gammaawidevm::Gammaawidevm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
941{
942 cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
943 read();
944 cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
945 wideResonanceCrossSection sigma(bbsystem);
946 sigma.crossSectionCalculation(_bwnormsave);
947 _VMbslope=sigma.slopeParameter();
948}
949
950
951//______________________________________________________________________________
952Gammaawidevm::~Gammaawidevm()
953{ }
954
955