]>
Commit | Line | Data |
---|---|---|
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 | // Nystrand 220710 | |
29 | // Fixed bug which gave incorrect minv distribution in gammagammaleptonpair. | |
30 | // Moved loop over W and Y from pickw to twoLeptonCrossSection in | |
31 | // gammagammaleptonpair to speed up event generation. | |
32 | // Changed to Woods-Saxon radius in twophotonluminosity to be consistent | |
33 | // with old starligt. | |
34 | // | |
35 | // | |
36 | /////////////////////////////////////////////////////////////////////////// | |
37 | ||
38 | ||
39 | #include <iostream> | |
40 | #include <fstream> | |
41 | #include <cmath> | |
42 | #include <vector> | |
43 | ||
44 | #include "starlightconstants.h" | |
45 | #include "gammagammaleptonpair.h" | |
46 | ||
47 | ||
48 | using namespace std; | |
49 | ||
50 | ||
51 | //_____________________________________________________________________________ | |
52 | Gammagammaleptonpair::Gammagammaleptonpair(beamBeamSystem& bbsystem) | |
53 | : eventChannel(bbsystem) | |
54 | { | |
55 | //Initialize randomgenerator with our seed. | |
56 | cout<<"Randy in leptonpair construction: "<<randyInstance.Rndom()<<endl; | |
57 | //Storing inputparameters into protected members for use | |
58 | _GGlepInputnumw=inputParametersInstance.nmbWBins(); | |
59 | _GGlepInputnumy=inputParametersInstance.nmbRapidityBins(); | |
60 | _GGlepInputpidtest=inputParametersInstance.prodParticleType(); | |
61 | _GGlepInputGamma_em=inputParametersInstance.beamLorentzGamma(); | |
62 | //Let us read in the luminosity tables | |
63 | read(); | |
64 | //Now we will calculate the crosssection | |
65 | twoLeptonCrossSection(); | |
66 | //If it is a tauon, calculate its tables | |
67 | if(inputParametersInstance.prodParticleId()==starlightConstants::TAUON) calculateTable(); | |
68 | } | |
69 | ||
70 | ||
71 | //______________________________________________________________________________ | |
72 | Gammagammaleptonpair::~Gammagammaleptonpair() | |
73 | { } | |
74 | ||
75 | ||
76 | //______________________________________________________________________________ | |
77 | void Gammagammaleptonpair::twoLeptonCrossSection() | |
78 | { | |
79 | //This function calculates section for 2-particle decay. For reference, see STAR Note 243, Eq. 9. | |
80 | //calculate the 2-lepton differential cross section | |
81 | //the 100 is to convert to barns | |
82 | //the 2 is to account for the fact that we integrate only over one half of the rapidity range | |
83 | //Multiply all _Farray[] by _f_max | |
84 | ||
85 | for(int i=0;i<_GGlepInputnumw;i++) | |
86 | { | |
87 | for(int j=0;j<_GGlepInputnumy;j++) | |
88 | { | |
89 | // _sigmax[i][j]=2.*Gammagammaleptonpair::twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/100.; | |
90 | _sigmax[i][j]=twoMuonCrossSection(_Warray[i])*_f_max*_Farray[i][j]/(100.*_Warray[i]); | |
91 | } | |
92 | } | |
93 | //calculate the total two-lepton cross section | |
94 | double sigmasum =0.; | |
95 | for(int i=0;i<_GGlepInputnumw-1;i++) | |
96 | { | |
97 | for(int j=0;j<_GGlepInputnumy-1;j++) | |
98 | { | |
99 | // _sigmaSum = _sigmaSum +2.*((_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i])/((_Warray[i+1]+_Warray[i])/2.)); | |
100 | // _sigmaSum = _sigmaSum +((_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i])/((_Warray[i+1]+_Warray[i])/2.)); | |
101 | sigmasum = sigmasum +(_sigmax[i][j]+_sigmax[i+1][j]+_sigmax[i][j+1]+_sigmax[i+1][j+1])/4.*(_Yarray[j+1]-_Yarray[j])*(_Warray[i+1]-_Warray[i]); | |
102 | } | |
103 | } | |
104 | cout << "The total "<<_GGlepInputpidtest<<" cross-section is: "<<sigmasum<<" barns."<<endl; | |
105 | ||
106 | // Do this integration here, once per run rather than once per event (JN 220710) | |
107 | //integrate sigma down to a function of just w | |
108 | ||
109 | double sgf=0.; | |
110 | ||
111 | for(int i=0;i<_ReadInputnumw;i++) | |
112 | { | |
113 | _sigofw[i]=0.; | |
114 | for(int j=0;j<_ReadInputnumy-1;j++) | |
115 | { | |
116 | _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.; | |
117 | } | |
118 | } | |
119 | ||
120 | //calculate the unnormalized sgfint array | |
121 | _sigfint[0]=0.; | |
122 | for(int i=0;i<_ReadInputnumw-1;i++) | |
123 | { | |
124 | sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.; | |
125 | _sigfint[i+1]=_sigfint[i]+sgf; | |
126 | } | |
127 | ||
128 | //normalize sgfint array | |
129 | _signormw=_sigfint[_ReadInputnumw-1]; | |
130 | for(int i=0;i<_ReadInputnumw;i++) | |
131 | { | |
132 | _sigfint[i]=_sigfint[i]/_signormw; | |
133 | } | |
134 | return; | |
135 | } | |
136 | ||
137 | ||
138 | //______________________________________________________________________________ | |
139 | double Gammagammaleptonpair::twoMuonCrossSection(double w) | |
140 | { | |
141 | //This function gives the two muon cross section as a function of Y&W. | |
142 | //Using the formula given in G.Soff et. al Nuclear Equation of State, part B, 579 | |
143 | double s=0.,Etest=0.,deltat=0.,xL=0.,sigmuix=0.,alphasquared=0.,hbarcsquared=0.; | |
144 | s = w*w; | |
145 | Etest = 4.*getMass()*getMass()/s; | |
146 | deltat = s * sqrt(1.-Etest); | |
147 | xL = 2.*log(sqrt(s)/(2.*getMass())+sqrt(1./Etest-1.)); | |
148 | alphasquared = starlightConstants::alpha*starlightConstants::alpha; | |
149 | hbarcsquared = starlightConstants::hbarc*starlightConstants::hbarc; | |
150 | sigmuix = 4.*starlightConstants::pi*alphasquared/s*hbarcsquared*((1+Etest-0.5*Etest*Etest)*xL-(1./s+Etest/s)*deltat); | |
151 | if(Etest > 1.) | |
152 | sigmuix = 0.; | |
153 | return sigmuix; | |
154 | } | |
155 | ||
156 | ||
157 | //______________________________________________________________________________ | |
158 | void Gammagammaleptonpair::pickw(double &w) | |
159 | { | |
160 | // This function picks a w for the 2-photon calculation. | |
161 | ||
162 | double x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.; | |
163 | int ivalw=0; | |
164 | ||
165 | if(_wdelta != 0) | |
166 | { | |
167 | w=_wdelta; | |
168 | ivalw=_ivalwd; | |
169 | remainw=_remainwd; | |
170 | } | |
171 | else{ | |
172 | //deal with the case where sigma is an array | |
173 | //_sigofw is simga integrated over y using a linear interpolation | |
174 | //sigint is the integral of sgfint, normalized | |
175 | ||
176 | //pick a random number | |
177 | x = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
178 | //compare x and sgfint to find the ivalue which is just less than the random number x | |
179 | for(int i=0;i<_GGlepInputnumw;i++) | |
180 | { | |
181 | if(x > _sigfint[i]) ivalw=i; | |
182 | } | |
183 | //remainder above ivalw | |
184 | remainarea = x - _sigfint[ivalw]; | |
185 | ||
186 | //figure out what point corresponds to the excess area in remainarea | |
187 | c = -remainarea*_signormw/(_Warray[ivalw+1]-_Warray[ivalw]); | |
188 | b = _sigofw[ivalw]; | |
189 | a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.; | |
190 | if(a==0.) | |
191 | { | |
192 | remainw = -c/b; | |
193 | } | |
194 | else{ | |
195 | remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a); | |
196 | } | |
197 | _ivalwd = ivalw; | |
198 | _remainwd = remainw; | |
199 | //calculate the w value | |
200 | w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw; | |
201 | ||
202 | } | |
203 | } | |
204 | ||
205 | ||
206 | //______________________________________________________________________________ | |
207 | void Gammagammaleptonpair::picky(double &y) | |
208 | { | |
209 | // This function picks a y given a W | |
210 | ||
211 | double * sigofy; | |
212 | double * sgfint; | |
213 | sigofy = new double[starlightLimits::MAXYBINS]; | |
214 | sgfint = new double[starlightLimits::MAXWBINS]; | |
215 | ||
216 | double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.; | |
217 | int ivalw=0,ivaly=0; | |
218 | ||
219 | ivalw=_ivalwd; | |
220 | remainw=_remainwd; | |
221 | //average over two colums to get y array | |
222 | for(int j=0;j<_GGlepInputnumy;j++) | |
223 | { | |
224 | sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw; | |
225 | } | |
226 | ||
227 | //calculate the unnormalized sgfint | |
228 | sgfint[0]=0.; | |
229 | for(int j=0;j<_GGlepInputnumy-1;j++) | |
230 | { | |
231 | sgf = (sigofy[j+1]+sigofy[j])/2.; | |
232 | sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]); | |
233 | } | |
234 | ||
235 | //normalize the sgfint array | |
236 | signorm = sgfint[_GGlepInputnumy-1]; | |
237 | ||
238 | for(int j=0;j<_GGlepInputnumy;j++) | |
239 | { | |
240 | sgfint[j]=sgfint[j]/signorm; | |
241 | } | |
242 | ||
243 | //pick a random number | |
244 | x = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
245 | //compare x and sgfint to find the ivalue which is just less then the random number x | |
246 | for(int i=0;i<_GGlepInputnumy;i++) | |
247 | { | |
248 | if(x > sgfint[i]) ivaly = i; | |
249 | } | |
250 | //remainder above ivaly | |
251 | remainarea = x - sgfint[ivaly]; | |
252 | ||
253 | //figure what point corresponds to the leftover area in remainarea | |
254 | c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]); | |
255 | b = sigofy[ivaly]; | |
256 | a = (sigofy[ivaly+1]-sigofy[ivaly])/2.; | |
257 | if(a==0.) | |
258 | { | |
259 | remainy = -c/b; | |
260 | } | |
261 | else{ | |
262 | remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a); | |
263 | } | |
264 | //calculate the y value | |
265 | y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy; | |
266 | delete[] sigofy; | |
267 | delete[] sgfint; | |
268 | } | |
269 | ||
270 | ||
271 | //______________________________________________________________________________ | |
272 | void Gammagammaleptonpair::pairMomentum(double w,double y,double &E,double &px,double &py,double &pz) | |
273 | { | |
274 | //this function calculates px,py,pz,and E given w and y | |
275 | ||
276 | double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.; | |
277 | ||
278 | //E1 and E2 are for the 2 photons in the CM frame | |
279 | E1 = w*exp(y)/2.; | |
280 | E2 = w*exp(-y)/2.; | |
281 | ||
282 | //calculate px and py | |
283 | //to get x and y components-- phi is random between 0 and 2*pi | |
284 | anglepp1 = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
285 | anglepp2 = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
286 | ||
287 | pp1 = pp_1(E1); | |
288 | pp2 = pp_2(E2); | |
289 | px = pp1*cos(2.*starlightConstants::pi*anglepp1)+pp2*cos(2.*starlightConstants::pi*anglepp2); | |
290 | py = pp1*sin(2.*starlightConstants::pi*anglepp1)+pp2*sin(2.*starlightConstants::pi*anglepp2); | |
291 | ||
292 | //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle | |
293 | pt = sqrt(px*px+py*py); | |
294 | ||
295 | //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz | |
296 | E = sqrt(w*w+pt*pt)*cosh(y); | |
297 | pz= sqrt(w*w+pt*pt)*sinh(y); | |
298 | signpx = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
299 | ||
300 | //pick the z direction | |
301 | //Don't do this anymore since y goes from -ymax to +ymax (JN 15-02-2013) | |
302 | //if(signpx > 0.5) pz = -pz; | |
303 | } | |
304 | ||
305 | ||
306 | //______________________________________________________________________________ | |
307 | double Gammagammaleptonpair::pp_1(double E) | |
308 | { | |
309 | // This is for beam 1 | |
310 | // returns on random draw from pp(E) distribution | |
311 | double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.; | |
312 | double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.; | |
313 | int satisfy =0; | |
314 | ||
315 | ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em); | |
316 | //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum | |
317 | Cm = sqrt(3.)*E/_GGlepInputGamma_em; | |
318 | //the amplitude of the p_t spectrum at the maximum | |
319 | singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds); | |
320 | Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm))); | |
321 | ||
322 | //pick a test value pp, and find the amplitude there | |
323 | x = randyInstance.Rndom(); | |
324 | pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); | |
325 | singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds); | |
326 | test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp))); | |
327 | ||
328 | while(satisfy==0){ | |
329 | u = randyInstance.Rndom(); | |
330 | if(u*Coef <= test) | |
331 | { | |
332 | satisfy =1; | |
333 | } | |
334 | else{ | |
335 | x =randyInstance.Rndom(); | |
336 | pp = 5*starlightConstants::hbarc/_bbs.beam1().nuclearRadius()*x; | |
337 | singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds); | |
338 | test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp)); | |
339 | } | |
340 | } | |
341 | ||
342 | return pp; | |
343 | } | |
344 | ||
345 | //______________________________________________________________________________ | |
346 | double Gammagammaleptonpair::pp_2(double E) | |
347 | { | |
348 | ||
349 | // This is for beam 2 | |
350 | //returns on random draw from pp(E) distribution | |
351 | double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.; | |
352 | double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.; | |
353 | int satisfy =0; | |
354 | ||
355 | ereds = (E/_GGlepInputGamma_em)*(E/_GGlepInputGamma_em); | |
356 | //sqrt(3)*E/gamma_em is p_t where the distribution is a maximum | |
357 | Cm = sqrt(3.)*E/_GGlepInputGamma_em; | |
358 | //the amplitude of the p_t spectrum at the maximum | |
359 | singleformfactorCm=_bbs.beam2().formFactor(Cm*Cm+ereds); | |
360 | Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm))); | |
361 | ||
362 | //pick a test value pp, and find the amplitude there | |
363 | x = randyInstance.Rndom(); | |
364 | pp = x*5.*starlightConstants::hbarc/_bbs.beam2().nuclearRadius(); //Will use nucleus #1 | |
365 | singleformfactorpp1=_bbs.beam2().formFactor(pp*pp+ereds); | |
366 | test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp))); | |
367 | ||
368 | while(satisfy==0){ | |
369 | u = randyInstance.Rndom(); | |
370 | if(u*Coef <= test) | |
371 | { | |
372 | satisfy =1; | |
373 | } | |
374 | else{ | |
375 | x =randyInstance.Rndom(); | |
376 | pp = 5*starlightConstants::hbarc/_bbs.beam2().nuclearRadius()*x; | |
377 | singleformfactorpp2=_bbs.beam2().formFactor(pp*pp+ereds); | |
378 | test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp)); | |
379 | } | |
380 | } | |
381 | ||
382 | return pp; | |
383 | } | |
384 | ||
385 | ||
386 | ||
387 | //______________________________________________________________________________ | |
388 | void Gammagammaleptonpair::twoBodyDecay(starlightConstants::particleTypeEnum &ipid, | |
389 | double , // E (unused) | |
390 | double W, | |
391 | double px0, double py0, double pz0, | |
392 | double& px1, double& py1, double& pz1, | |
393 | double& px2, double& py2, double& pz2, | |
394 | int& iFbadevent) | |
395 | { | |
396 | // This routine decays a particle into two particles of mass mdec, | |
397 | // taking spin into account | |
398 | ||
399 | double mdec=0.,E1=0.,E2=0.; | |
400 | double pmag, anglelep[20001]; | |
401 | // double ytest=0.,dndtheta; | |
402 | double phi,theta,xtest,Ecm; | |
403 | double betax,betay,betaz; | |
404 | double hirestheta,hirestest,hiresw; //added from JN->needed precision | |
405 | ||
406 | // set the mass of the daughter particles | |
407 | ||
408 | mdec = getMass(); | |
409 | if(W < 2*mdec) | |
410 | { | |
411 | cout<<" ERROR: W="<<W<<endl; | |
412 | iFbadevent = 1; | |
413 | return; | |
414 | } | |
415 | pmag = sqrt(W*W/4. - mdec*mdec); | |
416 | ||
417 | // pick an orientation, based on the spin | |
418 | // phi has a flat distribution in 2*pi | |
419 | phi = randyInstance.Rndom()*2.*starlightConstants::pi; //(random()/(RAND_MAX+1.0))* 2.*starlightConstants::pi; | |
420 | ||
421 | // find theta, the angle between one of the outgoing particles and | |
422 | // the beamline, in the frame of the two photons | |
423 | ||
424 | if(getSpin()==0.5){ | |
425 | // calculate a table of integrated angle values for leptons | |
426 | // JN05: Go from 1000->20000bins, use double precision for anglelep and thetalep. needed when W >6Gev. | |
427 | hiresw = W; | |
428 | ||
429 | anglelep[0] = 0.; | |
430 | ||
431 | for(int i =1;i<=20000;i++) | |
432 | { | |
433 | hirestheta = starlightConstants::pi * double(i) /20000.; | |
434 | ||
435 | // Added sin(theta) phase space factor (not in thetalep) and changed E to W in thetalep call | |
436 | // 11/9/2000 SRK | |
437 | // Note that thetalep is form invariant, so it could be called for E, theta_lab just | |
438 | // as well as W,theta_cm. Since there is a boost from cm to lab below, the former is fine. | |
439 | ||
440 | anglelep[i] = anglelep[i-1] + thetalep(hiresw,hirestheta)*sin(hirestheta); | |
441 | } | |
442 | ||
443 | hirestheta = 0.; | |
444 | xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
445 | hirestest = xtest; | |
446 | for(int i =1;i<=20000;i++) | |
447 | { | |
448 | if(xtest > (anglelep[i]/anglelep[20000])) | |
449 | hirestheta = starlightConstants::pi * double(i) / 20000.; | |
450 | } | |
451 | theta=hirestheta; | |
452 | ||
453 | } | |
454 | ||
455 | if(getSpin() != 0.5) | |
456 | cout<<" This model cannot yet handle this spin value for lepton pairs: "<<getSpin()<<endl; | |
457 | ||
458 | ||
459 | // compute unboosted momenta | |
460 | px1 = sin(theta)*cos(phi)*pmag; | |
461 | py1 = sin(theta)*sin(phi)*pmag; | |
462 | pz1 = cos(theta)*pmag; | |
463 | px2 = -px1; | |
464 | py2 = -py1; | |
465 | pz2 = -pz1; | |
466 | ||
467 | // compute energies | |
468 | //Changed mass to W 11/9/2000 SRK | |
469 | Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0); | |
470 | ||
471 | E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1); | |
472 | E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2); | |
473 | // decay tau to electrons | |
474 | // note that after this routine px1, etc., refer to the electrons | |
475 | if(_GGlepInputpidtest == starlightConstants::TAUON) | |
476 | tauDecay(px1,py1,pz1,E1,px2,py2,pz2,E2); | |
477 | ||
478 | // Lorentz transform into the lab frame | |
479 | // betax,betay,betaz are the boost of the complete system | |
480 | betax = -(px0/Ecm); | |
481 | betay = -(py0/Ecm); | |
482 | betaz = -(pz0/Ecm); | |
483 | ||
484 | transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent); | |
485 | transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent); | |
486 | ||
487 | ||
488 | if(iFbadevent == 1) | |
489 | return; | |
490 | ||
491 | // change particle id from that of parent to that of daughters | |
492 | // change taoun id into electron id, already switched particles in tau decay | |
493 | if(_GGlepInputpidtest == starlightConstants::TAUON) | |
494 | ipid = starlightConstants::ELECTRON; | |
495 | // electrons remain electrons; muons remain muons | |
496 | if ((_GGlepInputpidtest == starlightConstants::ELECTRON) || (_GGlepInputpidtest == starlightConstants::MUON)) | |
497 | ipid = _GGlepInputpidtest; | |
498 | } | |
499 | ||
500 | ||
501 | //______________________________________________________________________________ | |
502 | double Gammagammaleptonpair::thetalep(double W,double theta) | |
503 | { | |
504 | // This function calculates the cross-section as a function of | |
505 | // angle for a given W and Y, for the production of two muons. | |
506 | // (or tauons) | |
507 | // expression is taken from Brodsky et al. PRD 1971, 1532 | |
508 | // equation 5.7 | |
509 | // factor that are not dependant on theta are scrapped, so the | |
510 | // the absolute crosssections given by this function are inaccurate | |
511 | // here we are working in the CM frame of the photons and the last | |
512 | // term is 0 | |
513 | ||
514 | // real function thetalep (W,theta) | |
515 | ||
516 | double moverw=0., W1sq=0.; | |
517 | double thetalep_r=0.,denominator=0.; | |
518 | ||
519 | W1sq = (W / 2.)*(W/2.); | |
520 | moverw = getMass()*getMass() / W1sq; | |
521 | denominator = (1. - (1. - moverw)*(cos(theta)*cos(theta))); | |
522 | ||
523 | thetalep_r = 2. + 4.*(1.-moverw)*((1.-moverw)*sin(theta)*sin(theta)*cos(theta)*cos(theta) + moverw) / (denominator*denominator); | |
524 | ||
525 | return thetalep_r; | |
526 | ||
527 | } | |
528 | ||
529 | ||
530 | //______________________________________________________________________________ | |
531 | starlightConstants::event Gammagammaleptonpair::produceEvent(int &ievent) | |
532 | {//returns the vector with the decay particles inside. | |
533 | starlightConstants::event leptonpair; //This object will store all the tracks for a single event | |
534 | double comenergy = 0.; | |
535 | double rapidity = 0.; | |
536 | double pairE = 0.; | |
537 | double pairmomx=0.,pairmomy=0.,pairmomz=0.; | |
538 | int iFbadevent=0; | |
539 | starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN; | |
540 | ||
541 | double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.; | |
542 | //this function decays particles and writes events to a file | |
543 | //zero out the event structure | |
544 | leptonpair._numberOfTracks=0; | |
545 | for(int i=0;i<4;i++) | |
546 | { | |
547 | leptonpair.px[i]=0.; | |
548 | leptonpair.py[i]=0.; | |
549 | leptonpair.pz[i]=0.; | |
550 | leptonpair._fsParticle[i]=starlightConstants::UNKNOWN; | |
551 | leptonpair._charge[i]=0; | |
552 | } | |
553 | ||
554 | pickw(comenergy); | |
555 | ||
556 | picky(rapidity); | |
557 | ||
558 | pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz); | |
559 | twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent); | |
560 | if (iFbadevent==0){ | |
561 | int q1=0,q2=0; | |
562 | ||
563 | double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
564 | if (xtest<0.5) | |
565 | { | |
566 | q1=1; | |
567 | q2=-1; | |
568 | } | |
569 | else{ | |
570 | q1=-1; | |
571 | q2=1; | |
572 | } | |
573 | leptonpair._numberOfTracks=2;//leptonpairs are two tracks... | |
574 | leptonpair.px[0]=px1; | |
575 | leptonpair.py[0]=py1; | |
576 | leptonpair.pz[0]=pz1; | |
577 | leptonpair._fsParticle[0]=ipid; | |
578 | leptonpair._charge[0]=q1; | |
579 | ||
580 | leptonpair.px[1]=px2; | |
581 | leptonpair.py[1]=py2; | |
582 | leptonpair.pz[1]=pz2; | |
583 | leptonpair._fsParticle[1]=ipid; | |
584 | leptonpair._charge[1]=q2; | |
585 | ||
586 | ievent=ievent+1; | |
587 | } | |
588 | ||
589 | return leptonpair; | |
590 | } | |
591 | ||
592 | ||
593 | //______________________________________________________________________________ | |
594 | upcEvent Gammagammaleptonpair::produceEvent() | |
595 | { | |
596 | //returns the vector with the decay particles inside. | |
597 | ||
598 | upcEvent event; | |
599 | ||
600 | double comenergy = 0.; | |
601 | double rapidity = 0.; | |
602 | double pairE = 0.; | |
603 | double pairmomx=0.,pairmomy=0.,pairmomz=0.; | |
604 | int iFbadevent=0; | |
605 | starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN; | |
606 | ||
607 | double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.,E2=0.,E1=0.; | |
608 | bool accepted = false; | |
609 | do{ | |
610 | //this function decays particles and writes events to a file | |
611 | //zero out the event structure | |
612 | pickw(comenergy); | |
613 | ||
614 | picky(rapidity); | |
615 | ||
616 | pairMomentum(comenergy,rapidity,pairE,pairmomx,pairmomy,pairmomz); | |
617 | ||
618 | ||
619 | _nmbAttempts++; | |
620 | twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent); | |
621 | double pt1chk = sqrt(px1*px1+py1*py1); | |
622 | double pt2chk = sqrt(px2*px2+py2*py2); | |
623 | ||
624 | double eta1 = pseudoRapidity(px1, py1, pz1); | |
625 | double eta2 = pseudoRapidity(px2, py2, pz2); | |
626 | ||
627 | if(_ptCutEnabled && !_etaCutEnabled){ | |
628 | if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){ | |
629 | accepted = true; | |
630 | _nmbAccepted++; | |
631 | } | |
632 | } | |
633 | else if(!_ptCutEnabled && _etaCutEnabled){ | |
634 | if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){ | |
635 | accepted = true; | |
636 | _nmbAccepted++; | |
637 | } | |
638 | } | |
639 | else if(_ptCutEnabled && _etaCutEnabled){ | |
640 | if(pt1chk > _ptCutMin && pt1chk < _ptCutMax && pt2chk > _ptCutMin && pt2chk < _ptCutMax){ | |
641 | if(eta1 > _etaCutMin && eta1 < _etaCutMax && eta2 > _etaCutMin && eta2 < _etaCutMax){ | |
642 | accepted = true; | |
643 | _nmbAccepted++; | |
644 | } | |
645 | } | |
646 | } | |
647 | ||
648 | }while((_ptCutEnabled || _etaCutEnabled) && !accepted); | |
649 | //twoBodyDecay(ipid,pairE,comenergy,pairmomx,pairmomy,pairmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent); | |
650 | ||
651 | if (iFbadevent==0){ | |
652 | int q1=0,q2=0; | |
653 | ||
654 | double xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
655 | if (xtest<0.5) | |
656 | { | |
657 | q1=1; | |
658 | q2=-1; | |
659 | } | |
660 | else{ | |
661 | q1=-1; | |
662 | q2=1; | |
663 | } | |
664 | ||
665 | // The new stuff | |
666 | double mlepton = getMass(); | |
667 | E1 = sqrt( mlepton*mlepton + px1*px1 + py1*py1 + pz1*pz1 ); | |
668 | E2 = sqrt( mlepton*mlepton + px2*px2 + py2*py2 + pz2*pz2 ); | |
669 | ||
670 | starlightParticle particle1(px1, py1, pz1, E1, starlightConstants::UNKNOWN, -q1*ipid, q1); | |
671 | event.addParticle(particle1); | |
672 | ||
673 | starlightParticle particle2(px2, py2, pz2, E2, starlightConstants::UNKNOWN, -q2*ipid, q2); | |
674 | event.addParticle(particle2); | |
675 | ||
676 | } | |
677 | return event; | |
678 | } | |
679 | ||
680 | ||
681 | //______________________________________________________________________________ | |
682 | void Gammagammaleptonpair::calculateTable() | |
683 | { | |
684 | // this subroutine calculates the tables that are used | |
685 | // elsewhere in the montecarlo | |
686 | // the tauon decay is taken from V-A theory, 1 - 1/3 cos(theta) | |
687 | // the energy of each of the two leptons in tau decay | |
688 | // is calculated using formula 10.35 given | |
689 | // in introduction to elementary particles by D. griffiths | |
690 | // which assmes that the mass of the electron is 0. | |
691 | // the highest energy attainable by the electron in such a system is | |
692 | // .5 * mass of the tau | |
693 | ||
694 | // subroutine calculateTable | |
695 | ||
696 | ||
697 | double E,theta; | |
698 | ||
699 | _tautolangle[0] = 0.; | |
700 | _dgammade[0] = 0.; | |
701 | ||
702 | ||
703 | for(int i =1;i<=100;i++) | |
704 | { | |
705 | // calculate energy of tau decay | |
706 | E = double(i)/100. * .5 * starlightConstants::tauMass; | |
707 | _dgammade[i] = _dgammade[i-1] + E*E * (1. - 4.*E/(3.*starlightConstants::tauMass)); | |
708 | ||
709 | // calculate angles for tau | |
710 | theta = starlightConstants::pi * double(i) / 100.; | |
711 | _tautolangle[i] = _tautolangle[i-1] + (1 + 0.3333333 * cos(theta)); | |
712 | } | |
713 | ||
714 | ||
715 | } | |
716 | ||
717 | ||
718 | //______________________________________________________________________________ | |
719 | void Gammagammaleptonpair::tauDecay(double &px1,double &py1,double &pz1,double &E1,double &px2,double &py2,double &pz2,double &E2) | |
720 | { | |
721 | // this routine assumes that the tauons decay to electrons and | |
722 | // calculates the directons of the decays | |
723 | ||
724 | double Ee1,Ee2,theta1,theta2,phi1,phi2, ran1, ran2 ; | |
725 | double pmag1,pex1,pey1,pez1,pex2,pey2,pez2,pmag2; | |
726 | double betax,betay,betaz,dir; | |
727 | ||
728 | int Tmp_Par=0; // temp variable for the transform function .. kind of starnge - being called with 7 parameter instead of 8 | |
729 | ||
730 | // the highest energy attainable by the electron in this system is | |
731 | // .5 * mass of the tau | |
732 | ||
733 | // get two random numbers to compare with | |
734 | ||
735 | ||
736 | ran1 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100]; | |
737 | ran2 = randyInstance.Rndom()*_dgammade[100];//(random()/(RAND_MAX+1.0)) * _dgammade[100]; | |
738 | ||
739 | // compute the energies that correspond to those numbers | |
740 | Ee1 = 0.; | |
741 | Ee2 = 0.; | |
742 | ||
743 | for( int i =1;i<=100;i++) | |
744 | { | |
745 | if (ran1 > _dgammade[i]) | |
746 | Ee1 = double(i) /100. * .5 * getMass(); | |
747 | if (ran2 > _dgammade[i]) | |
748 | Ee2 = double(i) /100. * .5 * getMass(); | |
749 | } | |
750 | ||
751 | // to find the direction of emmission, first | |
752 | // we determine if the tauons have spin of +1 or -1 along the | |
753 | // direction of the beam line | |
754 | dir = 1.; | |
755 | if ( randyInstance.Rndom() < 0.5 )//(random()/(RAND_MAX+1.0)) < 0.5) | |
756 | dir = -1.; | |
757 | ||
758 | // get two random numbers to compare with | |
759 | ran1 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0)) * _tautolangle[100]; | |
760 | ran2 = randyInstance.Rndom()*_tautolangle[99];//(random()/(RAND_MAX+1.0)) * _tautolangle[100]; | |
761 | ||
762 | // find the angles corrsponding to those numbers | |
763 | theta1 = 0.; | |
764 | theta2 = 0.; | |
765 | for( int i =1;i<=100;i++) | |
766 | { | |
767 | if (ran1 > _tautolangle[i]) theta1 = starlightConstants::pi * double(i) /100.; | |
768 | if (ran2 > _tautolangle[i]) theta2 = starlightConstants::pi * double(i) /100.; | |
769 | } | |
770 | ||
771 | // grab another two random numbers to determine phi's | |
772 | phi1 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi; | |
773 | phi2 = randyInstance.Rndom()*2.*starlightConstants::pi;// (random()/(RAND_MAX+1.0))* 2. * starlightConstants::pi; | |
774 | // figure out the momenta of the electron in the frames of the | |
775 | // tauons from which they decayed, that is electron1 is in the | |
776 | // rest frame of tauon1 and e2 is in the rest fram of tau2 | |
777 | // again the electrons are assumed to be massless | |
778 | pmag1 = Ee1; | |
779 | pex1 = cos(phi1)*sin(theta1)*pmag1; | |
780 | pey1 = sin(phi1)*sin(theta1)*pmag1; | |
781 | pez1 = cos(theta1)*pmag1*dir; | |
782 | pmag2 = Ee2; | |
783 | pex2 = cos(phi2)*sin(theta2)*pmag2; | |
784 | pey2 = sin(phi2)*sin(theta2)*pmag2; | |
785 | pez2 = cos(theta2)*pmag2*(-1.*dir); | |
786 | // now Lorentz transform into the frame of each of the particles | |
787 | // do particle one first | |
788 | betax = -(px1/E1); | |
789 | betay = -(py1/E1); | |
790 | betaz = -(pz1/E1); | |
791 | //cout<<"2decay betax,pex1= "<<betax<<" "<<pex1<<endl; | |
792 | transform (betax,betay,betaz,Ee1,pex1,pey1,pez1,Tmp_Par); | |
793 | // then do particle two | |
794 | betax = -(px1/E2); | |
795 | betay = -(py1/E2); | |
796 | betaz = -(pz1/E2); | |
797 | ||
798 | transform (betax,betay,betaz,Ee2,pex2,pey2,pez2, Tmp_Par); | |
799 | // finally dump the electron values into the approriate | |
800 | // variables | |
801 | E1 = Ee1; | |
802 | E2 = Ee2; | |
803 | px1 = pex1; | |
804 | px2 = pex2; | |
805 | py1 = pey1; | |
806 | py2 = pey2; | |
807 | pz1 = pez1; | |
808 | pz2 = pez2; | |
809 | } | |
810 | ||
811 | ||
812 | //______________________________________________________________________________ | |
813 | double Gammagammaleptonpair::getMass() | |
814 | { | |
815 | double leptonmass=0.; | |
816 | switch(_GGlepInputpidtest){ | |
817 | case starlightConstants::ELECTRON: | |
818 | leptonmass=starlightConstants::mel; | |
819 | break; | |
820 | case starlightConstants::MUON: | |
821 | leptonmass=starlightConstants::muonMass; | |
822 | break; | |
823 | case starlightConstants::TAUON: | |
824 | leptonmass=starlightConstants::tauMass; | |
825 | break; | |
826 | default: | |
827 | cout<<"Not a recognized lepton, Gammagammaleptonpair::getmass(), mass = 0."<<endl; | |
828 | } | |
829 | ||
830 | return leptonmass; | |
831 | } | |
832 | ||
833 | ||
834 | //______________________________________________________________________________ | |
835 | double Gammagammaleptonpair::getWidth() | |
836 | { | |
837 | double leptonwidth=0.; | |
838 | return leptonwidth; | |
839 | ||
840 | } | |
841 | ||
842 | ||
843 | //______________________________________________________________________________ | |
844 | double Gammagammaleptonpair::getSpin() | |
845 | { | |
846 | double leptonspin=0.5; | |
847 | return leptonspin; | |
848 | } |