]> git.uio.no Git - u/mrichter/AliRoot.git/blame - STARLIGHT/starlight/src/gammaavm.cpp
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
STARLIGHT update:
[u/mrichter/AliRoot.git] / STARLIGHT / starlight / src / gammaavm.cpp
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:
75ce6a3a 23// $Rev:: 193 $: revision of last commit
0e1bd874 24// $Author:: jnystrand $: author of last commit
75ce6a3a 25// $Date:: 2014-12-01 20:39:46 +0100 #$: date of last commit
da32329d
AM		 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
75ce6a3a 68 N1 = 0; N2 = 0;
69
da32329d
AM		 70 switch(_VMpidtest){
71 case starlightConstants::RHO:
72 case starlightConstants::RHOZEUS:
73 _width=0.1507;
74 _mass=0.7685;
75 break;
76 case starlightConstants::FOURPRONG:
77 // create n-body phase-space generator instance
78 _phaseSpaceGen = new nBodyPhaseSpaceGen();
79 _width = 0.360;
80 _mass = 1.350;
81 break;
82 case starlightConstants::OMEGA:
83 _width=0.00843;
84 _mass=0.78194;
85 break;
86 case starlightConstants::PHI:
87 _width=0.00443;
88 _mass=1.019413;
89 break;
90 case starlightConstants::JPSI:
91 case starlightConstants::JPSI_ee:
92 case starlightConstants::JPSI_mumu:
93 _width=0.000091;
94 _mass=3.09692;
95 break;
96 case starlightConstants::JPSI2S:
97 case starlightConstants::JPSI2S_ee:
98 case starlightConstants::JPSI2S_mumu:
99 _width=0.000337;
100 _mass=3.686093;
101 break;
102 case starlightConstants::UPSILON:
103 case starlightConstants::UPSILON_ee:
104 case starlightConstants::UPSILON_mumu:
105 _width=0.00005402;
106 _mass=9.46030;
107 break;
108 case starlightConstants::UPSILON2S:
109 case starlightConstants::UPSILON2S_ee:
110 case starlightConstants::UPSILON2S_mumu:
111 _width=0.00003198;
112 _mass=10.02326;
113 break;
114 case starlightConstants::UPSILON3S:
115 case starlightConstants::UPSILON3S_ee:
116 case starlightConstants::UPSILON3S_mumu:
117 _width=0.00002032;
118 _mass=10.3552;
119 break;
120 default: cout<<"No proper vector meson defined, gammaavectormeson::gammaavectormeson"<<endl;
121 }
122}
123
124
125//______________________________________________________________________________
126Gammaavectormeson::~Gammaavectormeson()
127{
128 if (_phaseSpaceGen)
129 delete _phaseSpaceGen;
130}
131
132
133//______________________________________________________________________________
134void Gammaavectormeson::pickwy(double &W, double &Y)
135{
75ce6a3a 136 double dW, dY, xw,xy,xtest,btest;
da32329d
AM		 137 int IW,IY;
138
139 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
140 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
141
142 L201pwy:
143
144 xw = randyInstance.Rndom();// random()/(RAND_MAX+1.0);
145 W = _VMWmin + xw*(_VMWmax-_VMWmin);
146
147 if (W < 2 * starlightConstants::pionChargedMass)
148 goto L201pwy;
149
150 IW = int((W-_VMWmin)/dW); //+ 1;
151 xy = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
152 Y = _VMYmin + xy*(_VMYmax-_VMYmin);
153 IY = int((Y-_VMYmin)/dY); //+ 1;
154 xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
155
156 if( xtest > _Farray[IW][IY] )
157 goto L201pwy;
75ce6a3a 158
159 // Determine the target nucleus
160 // For pA this is given, for all other cases use the relative probabilities in _Farray1 and _Farray2
161 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
162 if( _ProductionMode == 2 ){
163 _TargetBeam = 2;
164 } else {
165 _TargetBeam = 1;
166 }
167 } else if( _bbs.beam1().A() != 1 && _bbs.beam2().A()==1 ){
168 if( _ProductionMode == 2 ){
169 _TargetBeam = 1;
170 } else {
171 _TargetBeam = 2;
172 }
173 } else {
174 btest = randyInstance.Rndom();
175 if ( btest < _Farray1[IW][IY]/_Farray[IW][IY] ){
176 _TargetBeam = 2;
177 N2++;
178 } else {
179 _TargetBeam = 1;
180 N1++;
181 }
182 }
183 // cout<<"N1: "<<N1<<" N2; "<<N2<<endl;
da32329d
AM		 184}
185
186
187//______________________________________________________________________________
188void Gammaavectormeson::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,
189 double, // E (unused)
190 double W,
191 double px0, double py0, double pz0,
192 double& px1, double& py1, double& pz1,
193 double& px2, double& py2, double& pz2,
194 int& iFbadevent)
195{
196 // This routine decays a particle into two particles of mass mdec,
197 // taking spin into account
198
199 double pmag;
200 // double anglelep[20001],xtest,ytest=0.,dndtheta;
201 double phi,theta,Ecm;
202 double betax,betay,betaz;
203 double mdec=0.0;
204 double E1=0.0,E2=0.0;
205
206 // set the mass of the daughter particles
207 mdec=getDaughterMass(ipid);
208
209 // calculate the magnitude of the momenta
210 if(W < 2*mdec){
211 cout<<" ERROR: W="<<W<<endl;
212 iFbadevent = 1;
213 return;
214 }
215 pmag = sqrt(W*W/4. - mdec*mdec);
216
217 // pick an orientation, based on the spin
218 // phi has a flat distribution in 2*pi
219 phi = randyInstance.Rndom()*2.*starlightConstants::pi;//(random()/(RAND_MAX+1.0))* 2.*pi;
220
221 // find theta, the angle between one of the outgoing particles and
222 // the beamline, in the frame of the two photons
223
224 theta=getTheta(ipid);
225
226 // compute unboosted momenta
227 px1 = sin(theta)*cos(phi)*pmag;
228 py1 = sin(theta)*sin(phi)*pmag;
229 pz1 = cos(theta)*pmag;
230 px2 = -px1;
231 py2 = -py1;
232 pz2 = -pz1;
233
234 Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0);
235 E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1);
236 E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2);
237
238 betax = -(px0/Ecm);
239 betay = -(py0/Ecm);
240 betaz = -(pz0/Ecm);
241
242 transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent);
243 transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent);
244
245 if(iFbadevent == 1)
246 return;
247
248}
249
250
251//______________________________________________________________________________
252// decays a particle into four particles with isotropic angular distribution
253bool Gammaavectormeson::fourBodyDecay
254(starlightConstants::particleTypeEnum& ipid,
255 const double , // E (unused)
256 const double W, // mass of produced particle
257 const double* p, // momentum of produced particle; expected to have size 3
258 lorentzVector* decayVecs, // array of Lorentz vectors of daughter particles; expected to have size 4
259 int& iFbadevent)
260{
261 const double parentMass = W;
262
263 // set the mass of the daughter particles
264 const double daughterMass = getDaughterMass(ipid);
265 if (parentMass < 4 * daughterMass){
266 cout << " ERROR: W=" << parentMass << " GeV too small" << endl;
267 iFbadevent = 1;
268 return false;
269 }
270
271 // construct parent four-vector
272 const double parentEnergy = sqrt(p[0] * p[0] + p[1] * p[1] + p[2] * p[2]
273 + parentMass * parentMass);
274 const lorentzVector parentVec(p[0], p[1], p[2], parentEnergy);
275
276 // setup n-body phase-space generator
277 assert(_phaseSpaceGen);
278 static bool firstCall = true;
279 if (firstCall) {
280 const double m[4] = {daughterMass, daughterMass, daughterMass, daughterMass};
281 _phaseSpaceGen->setDecay(4, m);
282 // estimate maximum phase-space weight
283 _phaseSpaceGen->setMaxWeight(1.01 * _phaseSpaceGen->estimateMaxWeight(_VMWmax));
284 firstCall = false;
285 }
286
287 // generate phase-space event
288 if (!_phaseSpaceGen->generateDecayAccepted(parentVec))
289 return false;
290
291 // set Lorentzvectors of decay daughters
292 for (unsigned int i = 0; i < 4; ++i)
293 decayVecs[i] = _phaseSpaceGen->daughter(i);
294 return true;
295}
296
297
298//______________________________________________________________________________
299double Gammaavectormeson::getDaughterMass(starlightConstants::particleTypeEnum &ipid)
300{
301 //This will return the daughter particles mass, and the final particles outputed id...
302 double ytest=0.,mdec=0.;
303
304 switch(_VMpidtest){
305 case starlightConstants::RHO:
306 case starlightConstants::RHOZEUS:
307 case starlightConstants::FOURPRONG:
308 case starlightConstants::OMEGA:
309 mdec = starlightConstants::pionChargedMass;
310 ipid = starlightConstants::PION;
311 break;
312 case starlightConstants::PHI:
313 mdec = starlightConstants::kaonChargedMass;
314 ipid = starlightConstants::KAONCHARGE;
315 break;
316 case starlightConstants::JPSI:
317 mdec = starlightConstants::mel;
318 ipid = starlightConstants::ELECTRON;
319 break;
320 case starlightConstants::JPSI_ee:
321 mdec = starlightConstants::mel;
322 ipid = starlightConstants::ELECTRON;
323 break;
324 case starlightConstants::JPSI_mumu:
325 mdec = starlightConstants::muonMass;
326 ipid = starlightConstants::MUON;
327 break;
328 case starlightConstants::JPSI2S_ee:
329 mdec = starlightConstants::mel;
330 ipid = starlightConstants::ELECTRON;
331 break;
332 case starlightConstants::JPSI2S_mumu:
333 mdec = starlightConstants::muonMass;
334 ipid = starlightConstants::MUON;
335 break;
336
337 case starlightConstants::JPSI2S:
338 case starlightConstants::UPSILON:
339 case starlightConstants::UPSILON2S:
340 case starlightConstants::UPSILON3S:
341 // decays 50% to e+/e-, 50% to mu+/mu-
342 ytest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
343
344 mdec = starlightConstants::muonMass;
345 ipid = starlightConstants::MUON;
346 break;
347 case starlightConstants::UPSILON_ee:
348 case starlightConstants::UPSILON2S_ee:
349 case starlightConstants::UPSILON3S_ee:
350 mdec = starlightConstants::mel;
351 ipid = starlightConstants::ELECTRON;
352 break;
353 case starlightConstants::UPSILON_mumu:
354 case starlightConstants::UPSILON2S_mumu:
355 case starlightConstants::UPSILON3S_mumu:
356 mdec = starlightConstants::muonMass;
357 ipid = starlightConstants::MUON;
358 break;
359 default: cout<<"No daughtermass defined, gammaavectormeson::getdaughtermass"<<endl;
360 }
361
362 return mdec;
363}
364
365
366//______________________________________________________________________________
367double Gammaavectormeson::getTheta(starlightConstants::particleTypeEnum ipid)
368{
369 //This depends on the decay angular distribution
370 //Valid for rho, phi, omega.
371 double theta=0.;
372 double xtest=0.;
373 double dndtheta=0.;
374 L200td:
375
376 theta = starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
377 xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0);
378 // Follow distribution for helicity +/-1
379 // Eq. 19 of J. Breitweg et al., Eur. Phys. J. C2, 247 (1998)
380 // SRK 11/14/2000
381
382 switch(ipid){
383
384 case starlightConstants::MUON:
385 case starlightConstants::ELECTRON:
386 //primarily for upsilon/j/psi. VM->ee/mumu
387 dndtheta = sin(theta)*(1.+((cos(theta))*(cos(theta))));
388 break;
389
390 case starlightConstants::PION:
391 case starlightConstants::KAONCHARGE:
392 //rhos etc
393 dndtheta= sin(theta)*(1.-((cos(theta))*(cos(theta))));
394 break;
395
396 default: cout<<"No proper theta dependence defined, check gammaavectormeson::gettheta"<<endl;
397 }//end of switch
398
399 if(xtest > dndtheta)
400 goto L200td;
401
402 return theta;
403
404}
405
406
407//______________________________________________________________________________
408double Gammaavectormeson::getWidth()
409{
410 return _width;
411}
412
413
414//______________________________________________________________________________
415double Gammaavectormeson::getMass()
416{
417 return _mass;
418}
419
420
421//______________________________________________________________________________
422double Gammaavectormeson::getSpin()
423{
424 return 1.0; //VM spins are the same
425}
426
427
428//______________________________________________________________________________
429void Gammaavectormeson::momenta(double W,double Y,double &E,double &px,double &py,double &pz,int &tcheck)
430{
431 // This subroutine calculates momentum and energy of vector meson
432 // given W and Y, without interference. Subroutine vmpt.f handles
433 // production with interference
434
435 double dW,dY;
436 double Egam,Epom,tmin,pt1,pt2,phi1,phi2;
437 double px1,py1,px2,py2;
438 double pt,xt,xtest,ytest;
439 // double photon_spectrum;
440 double t2;
441
442 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
443 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
444
445 //Find Egam,Epom in CM frame
446 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
447 if( _ProductionMode == 2 ){
448 Egam = 0.5*W*exp(Y);
449 Epom = 0.5*W*exp(-Y);
450 }else{
451 Egam = 0.5*W*exp(-Y);
452 Epom = 0.5*W*exp(Y);
453 }
454 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
455 if( _ProductionMode == 2 ){
456 Egam = 0.5*W*exp(-Y);
457 Epom = 0.5*W*exp(Y);
458 }else{
459 Egam = 0.5*W*exp(Y);
460 Epom = 0.5*W*exp(-Y);
461 }
75ce6a3a 462 } else if ( _TargetBeam == 1 ) {
463 Egam = 0.5*W*exp(-Y);
464 Epom = 0.5*W*exp(Y);
465 } else {
466 Egam = 0.5*W*exp(Y);
da32329d 467 Epom = 0.5*W*exp(-Y);
75ce6a3a 468 }
da32329d 469
75ce6a3a 470 // cout<<" Y: "<<Y<<" W: "<<W<<" TargetBeam; "<<_TargetBeam<<" Egam: "<<Egam<<" Epom: "<<Epom<<endl;
da32329d
AM		 471 pt1 = pTgamma(Egam);
472 phi1 = 2.*starlightConstants::pi*randyInstance.Rndom();
473
da32329d
AM		 474 if( (_bbs.beam1().A()==1 && _bbs.beam2().A()==1) ||
475 (_ProductionMode == 4) ) {
476 if( (_VMpidtest == starlightConstants::RHO) || (_VMpidtest == starlightConstants::RHOZEUS) || (_VMpidtest == starlightConstants::OMEGA)){
477 // Use dipole form factor for light VM
478 L555vm:
479 xtest = 2.0*randyInstance.Rndom();
480 double ttest = xtest*xtest;
481 ytest = randyInstance.Rndom();
482 double t0 = 1./2.23;
483 double yprob = xtest*_bbs.beam1().dipoleFormFactor(ttest,t0)*_bbs.beam1().dipoleFormFactor(ttest,t0);
484 if( ytest > yprob ) goto L555vm;
485 t2 = ttest;
486 pt2 = xtest;
487 }else{
488 //Use dsig/dt= exp(-_VMbslope*t) for heavy VM
489 xtest = randyInstance.Rndom();
490 t2 = (-1./_VMbslope)*log(xtest);
491 pt2 = sqrt(1.*t2);
492 }
493 } else {
494 // >> Check tmin
495 tmin = ((Epom/_VMgamma_em)*(Epom/_VMgamma_em));
496
497 if(tmin > 0.5){
498 cout<<" WARNING: tmin= "<<tmin<<endl;
499 cout<< " Y = "<<Y<<" W = "<<W<<" Epom = "<<Epom<<" gamma = "<<_VMgamma_em<<endl;
500 cout<<" Will pick a new W,Y "<<endl;
501 tcheck = 1;
502 return;
503 }
504 L203vm:
505 xt = randyInstance.Rndom();
506 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
507 if( _ProductionMode == 2 ){
508 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
509 }else{
510 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
511 }
512 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
513 if( _ProductionMode == 2 ){
514 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
515 }else{
516 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
517 }
75ce6a3a 518 } else if (_TargetBeam==1) {
519 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
da32329d
AM		 520 } else {
521 pt2 = 8.*xt*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
522 }
523
524 xtest = randyInstance.Rndom();
525 t2 = tmin + pt2*pt2;
526
527 if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
528 if(1.0 < _bbs.beam2().formFactor(t2)*pt2) cout <<"POMERON"<<endl;
529 if( xtest > _bbs.beam2().formFactor(t2)*pt2) goto L203vm;
530 }
531 else{
532
533 double comp=0.0;
534 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
535 if( _ProductionMode == 2 ){
536 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
537 }else{
538 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
539 }
540 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
541 if( _ProductionMode == 2 ){
542 comp = _bbs.beam1().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
543 }else{
544 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
545 }
75ce6a3a 546 } else if (_TargetBeam==1) {
547 comp = _bbs.beam2().formFactor(t2)*_bbs.beam1().formFactor(t2)*pt2;
da32329d 548 } else {
75ce6a3a 549 comp = _bbs.beam2().formFactor(t2)*_bbs.beam2().formFactor(t2)*pt2;
da32329d
AM		 550 }
551
552 if( xtest > comp ) goto L203vm;
553 }
554
da32329d
AM		 555 }//else end from pp
556 phi2 = 2.*starlightConstants::pi*randyInstance.Rndom();//random()/(RAND_MAX+1.0);
557
da32329d
AM		 558 px1 = pt1*cos(phi1);
559 py1 = pt1*sin(phi1);
560 px2 = pt2*cos(phi2);
561 py2 = pt2*sin(phi2);
562
563 // Compute vector sum Pt = Pt1 + Pt2 to find pt for the vector meson
564 px = px1 + px2;
565 py = py1 + py2;
566 pt = sqrt( px*px + py*py );
567
568 E = sqrt(W*W+pt*pt)*cosh(Y);
569 pz = sqrt(W*W+pt*pt)*sinh(Y);
570
571 // cout<< " Y = "<<Y<<" W = "<<W<<" Egam = "<<Egam<<" gamma = "<<_VMgamma_em<<endl;
572
75ce6a3a 573 // Keep this special case for d+A
da32329d
AM		 574 if(_bbs.beam2().Z()==1&&_bbs.beam2().A()==2){
575 pz = -pz;
75ce6a3a 576 }
577
578 /*
579 else if( (_bbs.beam1().A()==1 && _bbs.beam2().A() != 1) || (_bbs.beam2().A()==1 && _bbs.beam1().A() != 1) ){
da32329d
AM		 580 // Don't switch
581 }
582 else{
583 if (randyInstance.Rndom() >= 0.5) pz = -pz;
584 }
75ce6a3a 585 */
da32329d
AM		 586}
587
588//______________________________________________________________________________
589double Gammaavectormeson::pTgamma(double E)
590{
591 // returns on random draw from pp(E) distribution
592 double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.;
593 double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.;
594 int satisfy =0;
595
596 ereds = (E/_VMgamma_em)*(E/_VMgamma_em);
597 //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum
598 Cm = sqrt(3.)*E/_VMgamma_em;
0e1bd874 599 // If E is very small, the drawing of a pT below is extremely slow.
600 // ==> Set pT = sqrt(3.)*E/_VMgamma_em for very small E.
601 // Should have no observable consequences (JN, SRK 11-Sep-2014)
602 if( E < 0.0005 )return Cm;
603
da32329d
AM		 604 //the amplitude of the p_t spectrum at the maximum
605
606 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
607 if( _ProductionMode == 2 ){
608 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
609 }else{
610 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
611 }
612 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
613 if( _ProductionMode == 2 ){
614 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
615 }else{
616 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
617 }
75ce6a3a 618 } else if (_TargetBeam == 1) {
619 singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds);
da32329d
AM		 620 } else {
621 singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds);
622 }
623
624 Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm)));
625
626 //pick a test value pp, and find the amplitude there
627 x = randyInstance.Rndom();
628
629 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
630 if( _ProductionMode == 2 ){
631 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
632 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
633 }else{
634 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
635 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
636 }
637 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
638 if( _ProductionMode == 2 ){
639 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
640 singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds);
641 }else{
642 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
643 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
644 }
75ce6a3a 645 } else if (_TargetBeam == 1) {
646 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
647 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
da32329d
AM		 648 } else {
649 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
650 singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds);
651 }
652
653 test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp)));
654
655 while(satisfy==0){
656 u = randyInstance.Rndom();
657 if(u*Coef <= test)
658 {
659 satisfy =1;
660 }
661 else{
662 x =randyInstance.Rndom();
663 if( _bbs.beam1().A()==1 && _bbs.beam2().A() != 1){
664 if( _ProductionMode == 2 ){
665 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
666 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
667 }else{
668 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
669 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
670 }
671 } else if( _bbs.beam2().A()==1 && _bbs.beam1().A() != 1 ){
672 if( _ProductionMode == 2 ){
673 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
674 singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds);
675 }else{
676 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
677 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
678 }
75ce6a3a 679 } else if (_TargetBeam == 1) {
680 pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius();
681 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
da32329d
AM		 682 } else {
683 pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius();
684 singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);
685 }
686 test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp));
687 }
688 }
689
690 return pp;
691}
692
693
694//______________________________________________________________________________
695void Gammaavectormeson::vmpt(double W,double Y,double &E,double &px,double &py, double &pz,
696 int&) // tcheck (unused)
697{
698 // This function calculates momentum and energy of vector meson
75ce6a3a 699 // given W and Y, including interference.
700 // It gets the pt distribution from a lookup table.
da32329d
AM		 701 double dW=0.,dY=0.,yleft=0.,yfract=0.,xpt=0.,pt1=0.,ptfract=0.,pt=0.,pt2=0.,theta=0.;
702 int IY=0,j=0;
703
704 dW = (_VMWmax-_VMWmin)/double(_VMnumw);
705 dY = (_VMYmax-_VMYmin)/double(_VMnumy);
706
707 // Y is already fixed; choose a pt
708 // Follow the approavh in pickwy.f
75ce6a3a 709 // in _fptarray(IY,pt) IY=1 corresponds to Y=0, IY=numy/2 corresponds to +y
da32329d
AM		 710
711 IY=int(fabs(Y)/dY);//+1;
712 if (IY > (_VMnumy/2)-1){
713 IY=(_VMnumy/2)-1;
714 }
715
716 yleft=fabs(Y)-(IY)*dY;
717 yfract=yleft*dY;
75ce6a3a 718
da32329d 719 xpt=randyInstance.Rndom(); //random()/(RAND_MAX+1.0);
75ce6a3a 720
da32329d
AM		 721 for(j=0;j<_VMNPT+1;j++){
722 if (xpt < _fptarray[IY][j]) goto L60;
723 }
724 L60:
725
726 // now do linear interpolation - start with extremes
75ce6a3a 727 if (j == 0){
da32329d
AM		 728 pt1=xpt/_fptarray[IY][j]*_VMdpt/2.;
729 goto L80;
730 }
731 if (j == _VMNPT){
732 pt1=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY][j])/(1.-_fptarray[IY][j]);
733 goto L80;
734 }
735
736 // we're in the middle
75ce6a3a 737 ptfract=(xpt-_fptarray[IY][j])/(_fptarray[IY][j+1]-_fptarray[IY][j]);
da32329d
AM		 738 pt1=(j+1)*_VMdpt+ptfract*_VMdpt;
739
740 // at an extreme in y?
741 if (IY == (_VMnumy/2)-1){
742 pt=pt1;
743 goto L120;
744 }
745 L80:
75ce6a3a 746
747 // interpolate in y repeat for next fractional y bin
da32329d
AM		 748 for(j=0;j<_VMNPT+1;j++){
749 if (xpt < _fptarray[IY+1][j]) goto L90;
750 }
751 L90:
752
753 // now do linear interpolation - start with extremes
da32329d
AM		 754 if (j == 0){
755 pt2=xpt/_fptarray[IY+1][j]*_VMdpt/2.;
756 goto L100;
757 }
758 if (j == _VMNPT){
759 pt2=(_VMptmax-_VMdpt/2.) + _VMdpt/2.*(xpt-_fptarray[IY+1][j])/(1.-_fptarray[IY+1][j]);
760 goto L100;
761 }
762
763 // we're in the middle
da32329d
AM		 764 ptfract=(xpt-_fptarray[IY+1][j])/(_fptarray[IY+1][j+1]-_fptarray[IY+1][j]);
765 pt2=(j+1)*_VMdpt+ptfract*_VMdpt;
da32329d 766 L100:
75ce6a3a 767
768 // now interpolate in y
da32329d 769 pt=yfract*pt2+(1-yfract)*pt1;
75ce6a3a 770
da32329d 771 L120:
75ce6a3a 772
773 // we have a pt
da32329d
AM		 774 theta=2.*starlightConstants::pi*randyInstance.Rndom();//(random()/(RAND_MAX+1.0))*2.*pi;
775 px=pt*cos(theta);
776 py=pt*sin(theta);
75ce6a3a 777
da32329d
AM		 778 E = sqrt(W*W+pt*pt)*cosh(Y);
779 pz = sqrt(W*W+pt*pt)*sinh(Y);
780 // randomly choose to make pz negative 50% of the time
781 if(randyInstance.Rndom()>=0.5) pz = -pz;
782}
783
784
785//______________________________________________________________________________
786starlightConstants::event Gammaavectormeson::produceEvent(int&)
787{
788 // Not used; return default event
789 return starlightConstants::event();
790}
791
792
793//______________________________________________________________________________
794upcEvent Gammaavectormeson::produceEvent()
795{
796 // The new event type
797 upcEvent event;
798
799 int iFbadevent=0;
800 int tcheck=0;
801 starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN;
802 starlightConstants::particleTypeEnum vmpid = starlightConstants::UNKNOWN;
803
804 if (_VMpidtest == starlightConstants::FOURPRONG) {
805 double comenergy = 0;
806 double mom[3] = {0, 0, 0};
807 double E = 0;
808 lorentzVector decayVecs[4];
809 do {
810 double rapidity = 0;
811 pickwy(comenergy, rapidity);
812 if (_VMinterferencemode == 0)
813 momenta(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
814 else if (_VMinterferencemode==1)
815 vmpt(comenergy, rapidity, E, mom[0], mom[1], mom[2], tcheck);
816 } while (!fourBodyDecay(ipid, E, comenergy, mom, decayVecs, iFbadevent));
817 if ((iFbadevent == 0) and (tcheck == 0))
818 for (unsigned int i = 0; i < 4; ++i) {
819 starlightParticle daughter(decayVecs[i].GetPx(),
820 decayVecs[i].GetPy(),
821 decayVecs[i].GetPz(),
822 starlightConstants::UNKNOWN, // energy
823 starlightConstants::UNKNOWN, // _mass
824 ipid,
825 (i < 2) ? -1 : +1);
826 event.addParticle(daughter);
827 }
828 } else {
829 double comenergy = 0.;
830 double rapidity = 0.;
831 double E = 0.;
832 double momx=0.,momy=0.,momz=0.;
833
834 double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.;
835 bool accepted = false;
836 // if(_accCut){
837 do{
838 pickwy(comenergy,rapidity);
839
840 if (_VMinterferencemode==0){
841 momenta(comenergy,rapidity,E,momx,momy,momz,tcheck);
842 } else if (_VMinterferencemode==1){
843 vmpt(comenergy,rapidity,E,momx,momy,momz,tcheck);
844 }
845
846 // cout << "_ptCutMin: " << _ptCutMin << " _ptCutMax: " << _ptCutMax << " _etaCutMin: " << _etaCutMin << " _etaCutMax: " << _etaCutMax << endl;
847 _nmbAttempts++;
848 //cout << "n tries: " << _nmbAttempts<< endl;
849 vmpid = ipid;
850 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
851 double pt1chk = sqrt(px1*px1+py1*py1);
852 double pt2chk = sqrt(px2*px2+py2*py2);
853
854 //cout << "pt1: " << pt1chk << " pt2: " << pt2chk << endl;
855 double eta1 = pseudoRapidity(px1, py1, pz1);
856 double eta2 = pseudoRapidity(px2, py2, pz2);
857 //cout << "eta1: " << eta1 << " eta2: " << eta2 << endl;
858 if(_ptCutEnabled && !_etaCutEnabled){
859 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
860 accepted = true;
861 _nmbAccepted++;
862 }
863 }
864 else if(!_ptCutEnabled && _etaCutEnabled){
865 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
866 accepted = true;
867 _nmbAccepted++;
868 }
869 }
870 else if(_ptCutEnabled && _etaCutEnabled){
871 if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){
872 if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){
873 accepted = true;
874 _nmbAccepted++;
875 }
876 }
877 }
878 else if(!_ptCutEnabled && !_etaCutEnabled)
879 _nmbAccepted++;
880
881 }while((_ptCutEnabled || _etaCutEnabled) && !accepted);
882 /* }else{
883 twoBodyDecay(ipid,E,comenergy,momx,momy,momz,px1,py1,pz1,px2,py2,pz2,iFbadevent);
884 }*/
885 if (iFbadevent==0&&tcheck==0) {
886 int q1=0,q2=0;
887 int ipid1,ipid2=0;
888
889 double xtest = randyInstance.Rndom();
890 if (xtest<0.5)
891 {
892 q1=1;
893 q2=-1;
894 }
895 else {
896 q1=-1;
897 q2=1;
898 }
899
900 if ( ipid == 11 || ipid == 13 ){
901 ipid1 = -q1*ipid;
902 ipid2 = -q2*ipid;
903 } else {
904 ipid1 = q1*ipid;
905 ipid2 = q2*ipid;
906 }
907 // The new stuff
908 double md = getDaughterMass(vmpid);
909 double Ed1 = sqrt(md*md+px1*px1+py1*py1+pz1*pz1);
910 starlightParticle particle1(px1, py1, pz1, Ed1, starlightConstants::UNKNOWN, ipid1, q1);
911 event.addParticle(particle1);
912
913 double Ed2 = sqrt(md*md+px2*px2+py2*py2+pz2*pz2);
914 starlightParticle particle2(px2, py2, pz2, Ed2, starlightConstants::UNKNOWN, ipid2, q2);
915 event.addParticle(particle2);
916 // End of the new stuff
917
918 }
919 }
920
921 return event;
922
923}
924double Gammaavectormeson::pseudoRapidity(double px, double py, double pz)
925{
926 double pT = sqrt(px*px + py*py);
927 double p = sqrt(pz*pz + pT*pT);
928 double eta = -99.9; if((p-pz) != 0){eta = 0.5*log((p+pz)/(p-pz));}
929 return eta;
930}
931
932//______________________________________________________________________________
933Gammaanarrowvm::Gammaanarrowvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
934{
935 cout<<"Reading in luminosity tables. Gammaanarrowvm()"<<endl;
936 read();
937 cout<<"Creating and calculating crosssection. Gammaanarrowvm()"<<endl;
938 narrowResonanceCrossSection sigma(bbsystem);
939 sigma.crossSectionCalculation(_bwnormsave);
940 _VMbslope=sigma.slopeParameter();
941}
942
943
944//______________________________________________________________________________
945Gammaanarrowvm::~Gammaanarrowvm()
946{ }
947
948
949//______________________________________________________________________________
950Gammaaincoherentvm::Gammaaincoherentvm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
951{
952 cout<<"Reading in luminosity tables. Gammaainkoherentvm()"<<endl;
953 read();
954 cout<<"Creating and calculating crosssection. Gammaainkoherentvm()"<<endl;
955 incoherentVMCrossSection sigma(bbsystem); sigma.crossSectionCalculation(_bwnormsave);
956 _VMbslope=sigma.slopeParameter();
957}
958
959
960//______________________________________________________________________________
961Gammaaincoherentvm::~Gammaaincoherentvm()
962{ }
963
964
965//______________________________________________________________________________
966Gammaawidevm::Gammaawidevm(beamBeamSystem& bbsystem):Gammaavectormeson(bbsystem)
967{
968 cout<<"Reading in luminosity tables. Gammaawidevm()"<<endl;
969 read();
970 cout<<"Creating and calculating crosssection. Gammaawidevm()"<<endl;
971 wideResonanceCrossSection sigma(bbsystem);
972 sigma.crossSectionCalculation(_bwnormsave);
973 _VMbslope=sigma.slopeParameter();
974}
975
976
977//______________________________________________________________________________
978Gammaawidevm::~Gammaawidevm()
979{ }
980
981