]>
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 | // | |
29 | // | |
30 | // | |
31 | /////////////////////////////////////////////////////////////////////////// | |
32 | ||
33 | #include <iostream> | |
34 | #include <fstream> | |
35 | #include <cmath> | |
36 | #include <vector> | |
37 | ||
38 | ||
39 | #include "starlightconstants.h" | |
40 | #include "gammagammasingle.h" | |
41 | #include "starlightconfig.h" | |
42 | ||
43 | using namespace std; | |
44 | ||
45 | ||
46 | //______________________________________________________________________________ | |
47 | Gammagammasingle::Gammagammasingle(beamBeamSystem& bbsystem) | |
48 | : eventChannel(bbsystem) | |
49 | #ifdef ENABLE_PYTHIA | |
50 | ,_pyDecayer() | |
51 | #endif | |
52 | { | |
53 | ||
54 | #ifdef ENABLE_PYTHIA | |
55 | _pyDecayer.init(); | |
56 | #endif | |
57 | ||
58 | //Storing inputparameters into protected members for use | |
59 | _GGsingInputnumw=inputParametersInstance.nmbWBins(); | |
60 | _GGsingInputnumy=inputParametersInstance.nmbRapidityBins(); | |
61 | _GGsingInputpidtest=inputParametersInstance.prodParticleType(); | |
62 | _GGsingInputGamma_em=inputParametersInstance.beamLorentzGamma(); | |
63 | cout<<"SINGLE MESON pid test: "<<_GGsingInputpidtest<<endl; | |
64 | //reading in luminosity tables | |
65 | read(); | |
66 | //Now calculating crosssection | |
67 | singleCrossSection(); | |
68 | } | |
69 | ||
70 | ||
71 | //______________________________________________________________________________ | |
72 | Gammagammasingle::~Gammagammasingle() | |
73 | { } | |
74 | ||
75 | ||
76 | //______________________________________________________________________________ | |
77 | void Gammagammasingle::singleCrossSection() | |
78 | { | |
79 | //This function carries out a delta function cross-section calculation. For reference, see STAR Note 243, Eq. 8 | |
80 | //Multiply all _Farray[] by _f_max | |
81 | double _sigmaSum=0.,remainw=0.;//_remainwd=0.; | |
82 | int ivalw =0;//_ivalwd; | |
83 | //calculate the differential cross section and place in the sigma table | |
84 | _wdelta=getMass(); | |
85 | for(int i=0;i<_GGsingInputnumw;i++){ | |
86 | for(int j=0;j<_GGsingInputnumy;j++){ | |
87 | // Eq. 1 of starnote 347 | |
88 | _sigmax[i][j]=(getSpin()*2.+1.)*4*starlightConstants::pi*starlightConstants::pi*getWidth()/ | |
89 | (getMass()*getMass()*getMass())*_f_max*_Farray[i][j]*starlightConstants::hbarc*starlightConstants::hbarc/100.; | |
90 | } | |
91 | } | |
92 | //find the index, i,for the value of w just less than the mass because we want to use the value from the sigma table that has w=mass | |
93 | ||
94 | for(int i=0;i<_GGsingInputnumw;i++){ | |
95 | if(getMass()>_Warray[i]) ivalw=i; | |
96 | } | |
97 | ||
98 | remainw = (getMass()-_Warray[ivalw])/(_Warray[ivalw+1]-_Warray[ivalw+1]-_Warray[ivalw]); | |
99 | _ivalwd = ivalw; | |
100 | _remainwd = remainw; | |
101 | //if we are interested rho pairs at threshold, the just set sigma to 100nb | |
102 | switch(_GGsingInputpidtest){ | |
103 | case starlightConstants::ZOVERZ03: | |
104 | _sigmaSum =0.; | |
105 | for(int j=0;j<_GGsingInputnumy-1;j++){ | |
106 | _sigmaSum = _sigmaSum +(_Yarray[j+1]-_Yarray[j])* | |
107 | 100.0E-9*(.1/getMass())*((1.-remainw)*_f_max* | |
108 | (_Farray[ivalw][j]+_Farray[ivalw][j])/2.+remainw*_f_max* | |
109 | (_Farray[ivalw+1][j]+_Farray[ivalw+1][j+1])/2.); | |
110 | } | |
111 | break; | |
112 | default: | |
113 | //Sum to find the total cross-section | |
114 | _sigmaSum =0.; | |
115 | for(int j =0;j<_GGsingInputnumy-1;j++){ | |
116 | _sigmaSum = _sigmaSum+ | |
117 | (_Yarray[j+1]-_Yarray[j])*((1.-remainw)* | |
118 | (_sigmax[ivalw][j]+_sigmax[ivalw][j+1])/2.+remainw* | |
119 | (_sigmax[ivalw+1][j]+_sigmax[ivalw+1][j+1])/2.); | |
120 | } | |
121 | } | |
122 | if(_sigmaSum > 0.1) cout <<"The total cross-section is: "<<_sigmaSum<<" barns."<<endl; | |
123 | else if(_sigmaSum > 0.0001)cout <<"The total cross-section is: "<<_sigmaSum*1000<<" mb."<<endl; | |
124 | else cout <<"The total cross-section is: "<<_sigmaSum*1000000<<" ub."<<endl; | |
125 | ||
126 | ||
127 | return; | |
128 | } | |
129 | ||
130 | ||
131 | //______________________________________________________________________________ | |
132 | void Gammagammasingle::pickw(double &w) | |
133 | { | |
134 | //This function picks a w for the 2-photon calculation. | |
135 | double sgf=0.,signorm=0.,x=0.,remainarea=0.,remainw=0.,a=0.,b=0.,c=0.; | |
136 | int ivalw=0; | |
137 | ||
138 | double * _sigofw; | |
139 | double * sgfint; | |
140 | _sigofw = new double[starlightLimits::MAXWBINS]; | |
141 | sgfint = new double[starlightLimits::MAXYBINS]; | |
142 | ||
143 | if(_wdelta != 0){ | |
144 | w=_wdelta; | |
145 | ivalw=_ivalwd; | |
146 | remainw=_remainwd; | |
147 | } | |
148 | else{ | |
149 | //deal with the case where sigma is an array | |
150 | //_sigofw is simga integrated over y using a linear interpolation | |
151 | //sigint is the integral of sgfint, normalized | |
152 | ||
153 | //integrate sigma down to a function of just w | |
154 | for(int i=0;i<_GGsingInputnumw;i++){ | |
155 | _sigofw[i]=0.; | |
156 | for(int j=0;j<_GGsingInputnumy-1;j++){ | |
157 | _sigofw[i] = _sigofw[i]+(_Yarray[j+1]-_Yarray[j])*(_sigmax[i][j+1]+_sigmax[i][j])/2.; | |
158 | } | |
159 | } | |
160 | //calculate the unnormalized sgfint array | |
161 | sgfint[0]=0.; | |
162 | for(int i=0;i<_GGsingInputnumw-1;i++){ | |
163 | sgf=(_sigofw[i+1]+_sigofw[i])*(_Warray[i+1]-_Warray[i])/2.; | |
164 | sgfint[i+1]=sgfint[i]+sgf; | |
165 | } | |
166 | //normalize sgfint array | |
167 | signorm=sgfint[_GGsingInputnumw-1]; | |
168 | ||
169 | for(int i=0;i<_GGsingInputnumw;i++){ | |
170 | sgfint[i]=sgfint[i]/signorm; | |
171 | } | |
172 | //pick a random number | |
173 | x = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
174 | //compare x and sgfint to find the ivalue which is just less than the random number x | |
175 | for(int i=0;i<_GGsingInputnumw;i++){ | |
176 | if(x > sgfint[i]) ivalw=i; | |
177 | } | |
178 | //remainder above ivalw | |
179 | remainarea = x - sgfint[ivalw]; | |
180 | ||
181 | //figure out what point corresponds to the excess area in remainarea | |
182 | c = -remainarea*signorm/(_Warray[ivalw+1]-_Warray[ivalw]); | |
183 | b = _sigofw[ivalw]; | |
184 | a = (_sigofw[ivalw+1]-_sigofw[ivalw])/2.; | |
185 | if(a==0.){ | |
186 | remainw = -c/b; | |
187 | } | |
188 | else{ | |
189 | remainw = (-b+sqrt(b*b-4.*a*c))/(2.*a); | |
190 | } | |
191 | _ivalwd = ivalw; | |
192 | _remainwd = remainw; | |
193 | //calculate the w value | |
194 | w = _Warray[ivalw]+(_Warray[ivalw+1]-_Warray[ivalw])*remainw; | |
195 | } | |
196 | ||
197 | delete[] _sigofw; | |
198 | delete[] sgfint; | |
199 | } | |
200 | ||
201 | ||
202 | //______________________________________________________________________________ | |
203 | void Gammagammasingle::picky(double &y) | |
204 | { | |
205 | double * sigofy; | |
206 | double * sgfint; | |
207 | sigofy = new double[starlightLimits::MAXYBINS]; | |
208 | sgfint = new double[starlightLimits::MAXYBINS]; | |
209 | ||
210 | double remainw =0.,remainarea=0.,remainy=0.,a=0.,b=0.,c=0.,sgf=0.,signorm=0.,x=0.; | |
211 | int ivalw=0,ivaly=0; | |
212 | ||
213 | ivalw=_ivalwd; | |
214 | remainw=_remainwd; | |
215 | //average over two colums to get y array | |
216 | for(int j=0;j<_GGsingInputnumy;j++){ | |
217 | sigofy[j]=_sigmax[ivalw][j]+(_sigmax[ivalw+1][j]-_sigmax[ivalw][j])*remainw; | |
218 | } | |
219 | //calculate the unnormalized sgfint | |
220 | ||
221 | sgfint[0]=0.; | |
222 | for(int j=0;j<_GGsingInputnumy-1;j++){ | |
223 | sgf = (sigofy[j+1]+sigofy[j])/2.; | |
224 | sgfint[j+1]=sgfint[j]+sgf*(_Yarray[j+1]-_Yarray[j]); | |
225 | } | |
226 | ||
227 | //normalize the sgfint array | |
228 | signorm = sgfint[_GGsingInputnumy-1]; | |
229 | ||
230 | for(int j=0;j<_GGsingInputnumy;j++){ | |
231 | sgfint[j]=sgfint[j]/signorm; | |
232 | } | |
233 | //pick a random number | |
234 | x = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
235 | //compare x and sgfint to find the ivalue which is just less then the random number x | |
236 | for(int i=0;i<_GGsingInputnumy;i++){ | |
237 | if(x > sgfint[i]) | |
238 | ivaly = i; | |
239 | } | |
240 | //remainder above ivaly | |
241 | remainarea = x - sgfint[ivaly]; | |
242 | //figure what point corresponds to the leftover area in remainarea | |
243 | c = -remainarea*signorm/(_Yarray[ivaly+1]-_Yarray[ivaly]); | |
244 | b = sigofy[ivaly]; | |
245 | a = (sigofy[ivaly+1]-sigofy[ivaly])/2.; | |
246 | if(a==0.){ | |
247 | remainy = -c/b; | |
248 | } | |
249 | else{ | |
250 | remainy = (-b + sqrt(b*b-4.*a*c))/(2.*a); | |
251 | } | |
252 | //calculate the y value | |
253 | y = _Yarray[ivaly]+(_Yarray[ivaly+1]-_Yarray[ivaly])*remainy; | |
254 | delete[] sigofy; | |
255 | delete[] sgfint; | |
256 | } | |
257 | ||
258 | ||
259 | //______________________________________________________________________________ | |
260 | void Gammagammasingle::parentMomentum(double w,double y,double &E,double &px,double &py,double &pz) | |
261 | { | |
262 | //this function calculates px,py,pz,and E given w and y | |
263 | double anglepp1=0.,anglepp2=0.,pp1=0.,pp2=0.,E1=0.,E2=0.,signpx=0.,pt=0.; | |
264 | ||
265 | //E1 and E2 are for the 2 photons in the CM frame | |
266 | E1 = w*exp(y)/2.; | |
267 | E2 = w*exp(-y)/2.; | |
268 | //pz = E1-E2; | |
269 | //calculate px and py | |
270 | //to get x and y components-- phi is random between 0 and 2*pi | |
271 | anglepp1 = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
272 | anglepp2 = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
273 | ||
274 | pp1 = pp(E1); | |
275 | pp2 = pp(E2); | |
276 | px = pp1*cos(2.*starlightConstants::pi*anglepp1)+pp2*cos(2.*starlightConstants::pi*anglepp2); | |
277 | py = pp1*sin(2.*starlightConstants::pi*anglepp1)+pp2*sin(2.*starlightConstants::pi*anglepp2); | |
278 | //Compute vector sum Pt=Pt1+Pt2 to find pt for the produced particle | |
279 | pt = sqrt(px*px+py*py); | |
280 | //W is the mass of the produced particle (not necessarily on-mass-shell).Now compute its energy and pz | |
281 | E = sqrt(w*w+pt*pt)*cosh(y); | |
282 | pz= sqrt(w*w+pt*pt)*sinh(y); | |
283 | signpx = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
284 | //pick the z direction | |
285 | if(signpx > 0.5) | |
286 | pz = -pz; | |
287 | } | |
288 | ||
289 | ||
290 | //______________________________________________________________________________ | |
291 | double Gammagammasingle::pp(double E) | |
292 | { | |
293 | // will probably have to pass in beambeamsys? that way we can do beam1.formFactor(t) or beam2..., careful with the way sergey did it for asymmetry | |
294 | // returns on random draw from pp(E) distribution | |
295 | ||
296 | double ereds =0.,Cm=0.,Coef=0.,x=0.,pp=0.,test=0.,u=0.; | |
297 | double singleformfactorCm=0.,singleformfactorpp1=0.,singleformfactorpp2=0.; | |
298 | int satisfy =0; | |
299 | ||
300 | ereds = (E/_GGsingInputGamma_em)*(E/_GGsingInputGamma_em); | |
301 | Cm = sqrt(3.)*E/_GGsingInputGamma_em; | |
302 | //the amplitude of the p_t spectrum at the maximum | |
303 | singleformfactorCm=_bbs.beam1().formFactor(Cm*Cm+ereds); | |
304 | //Doing this once and then storing it as a double, which we square later...SYMMETRY?using beam1 for now. | |
305 | Coef = 3.0*(singleformfactorCm*singleformfactorCm*Cm*Cm*Cm)/((2.*(starlightConstants::pi)*(ereds+Cm*Cm))*(2.*(starlightConstants::pi)*(ereds+Cm*Cm))); | |
306 | ||
307 | //pick a test value pp, and find the amplitude there | |
308 | x = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
309 | pp = x*5.*starlightConstants::hbarc/_bbs.beam1().nuclearRadius(); //Will use nucleus #1, there should be two for symmetry//nextline | |
310 | singleformfactorpp1=_bbs.beam1().formFactor(pp*pp+ereds); | |
311 | test = (singleformfactorpp1*singleformfactorpp1)*pp*pp*pp/((2.*starlightConstants::pi*(ereds+pp*pp))*(2.*starlightConstants::pi*(ereds+pp*pp))); | |
312 | ||
313 | while(satisfy==0){ | |
314 | u = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
315 | if(u*Coef <= test){ | |
316 | satisfy =1; | |
317 | } | |
318 | else{ | |
319 | x =randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
320 | pp = 5*starlightConstants::hbarc/_bbs.beam1().nuclearRadius()*x; | |
321 | singleformfactorpp2=_bbs.beam1().formFactor(pp*pp+ereds);//Symmetry | |
322 | test = (singleformfactorpp2*singleformfactorpp2)*pp*pp*pp/(2.*starlightConstants::pi*(ereds+pp*pp)*2.*starlightConstants::pi*(ereds+pp*pp)); | |
323 | } | |
324 | } | |
325 | return pp; | |
326 | } | |
327 | ||
328 | ||
329 | //______________________________________________________________________________ | |
330 | void Gammagammasingle::twoBodyDecay(starlightConstants::particleTypeEnum &ipid,double /*E*/,double W,double px0,double py0,double pz0,double &px1,double &py1,double &pz1,double &px2,double &py2,double &pz2,int &iFbadevent) | |
331 | { | |
332 | // This routine decays a particle into two particles of mass mdec, | |
333 | // taking spin into account | |
334 | ||
335 | double mdec=0.,E1=0.,E2=0.; | |
336 | double pmag,ytest=0.; | |
337 | double phi,theta,xtest,dndtheta,Ecm; | |
338 | double betax,betay,betaz; | |
339 | ||
340 | // set the mass of the daughter particles | |
341 | switch(_GGsingInputpidtest){ | |
342 | case starlightConstants::ZOVERZ03: | |
343 | case starlightConstants::F2: | |
344 | mdec = starlightConstants::pionChargedMass; | |
345 | break; | |
346 | case starlightConstants::F2PRIME: | |
347 | // decays 50% to K+/K-, 50% to K_0's | |
348 | ytest = randyInstance.Rndom(); | |
349 | if(ytest >= 0.5){ | |
350 | mdec = starlightConstants::kaonChargedMass; | |
351 | } | |
352 | else{ | |
353 | mdec = 0.493677; | |
354 | } | |
355 | break; | |
356 | default : | |
357 | cout<<"No default mass selected for single photon-photon particle, expect errant results"<<endl; | |
358 | } | |
359 | ||
360 | //Calculating the momentum's magnitude | |
361 | //add switch for rho pairs at threshold and everything else. | |
362 | switch(_GGsingInputpidtest){ | |
363 | case starlightConstants::ZOVERZ03: //the rho pairs produced at threshold | |
364 | pmag = sqrt(getMass()*getMass()/4. - mdec*mdec); | |
365 | break; | |
366 | default : | |
367 | if(W < 2*mdec){ | |
368 | cout<<" ERROR: W="<<W<<endl; | |
369 | iFbadevent = 1; | |
370 | return; | |
371 | } | |
372 | pmag = sqrt(W*W/4. - mdec*mdec); | |
373 | } | |
374 | // pick an orientation, based on the spin | |
375 | // phi has a flat distribution in 2*pi | |
376 | phi = randyInstance.Rndom()*2.*starlightConstants::pi; //(random()/(RAND_MAX+1.0))* 2.*starlightConstants::pi; | |
377 | ||
378 | // find theta, the angle between one of the outgoing particles and | |
379 | // the beamline, in the frame of the two photons | |
380 | //this will depend on spin, F2,F2' and z/z03 all have spin 2, all other photonphoton-single mesons are handled by jetset | |
381 | //Applies to spin2 mesons. | |
382 | L300td: | |
383 | theta = starlightConstants::pi*randyInstance.Rndom(); | |
384 | xtest = randyInstance.Rndom(); | |
385 | dndtheta = sin(theta)*sin(theta)*sin(theta)*sin(theta)*sin(theta); | |
386 | if(xtest > dndtheta) | |
387 | goto L300td; | |
388 | ||
389 | // compute unboosted momenta | |
390 | px1 = sin(theta)*cos(phi)*pmag; | |
391 | py1 = sin(theta)*sin(phi)*pmag; | |
392 | pz1 = cos(theta)*pmag; | |
393 | px2 = -px1; | |
394 | py2 = -py1; | |
395 | pz2 = -pz1; | |
396 | // compute energies | |
397 | //Changed mass to W 11/9/2000 SRK | |
398 | Ecm = sqrt(W*W+px0*px0+py0*py0+pz0*pz0); | |
399 | E1 = sqrt(mdec*mdec+px1*px1+py1*py1+pz1*pz1); | |
400 | E2 = sqrt(mdec*mdec+px2*px2+py2*py2+pz2*pz2); | |
401 | ||
402 | // Lorentz transform into the lab frame | |
403 | // betax,betay,betaz are the boost of the complete system | |
404 | betax = -(px0/Ecm); | |
405 | betay = -(py0/Ecm); | |
406 | betaz = -(pz0/Ecm); | |
407 | ||
408 | transform (betax,betay,betaz,E1,px1,py1,pz1,iFbadevent); | |
409 | transform (betax,betay,betaz,E2,px2,py2,pz2,iFbadevent); | |
410 | ||
411 | ||
412 | if(iFbadevent == 1) | |
413 | return; | |
414 | ||
415 | // change particle id from that of parent to that of daughters | |
416 | ||
417 | switch(_GGsingInputpidtest){ | |
418 | //These decay into a pi+ pi- pair | |
419 | case starlightConstants::ZOVERZ03: | |
420 | case starlightConstants::F2: | |
421 | ipid=starlightConstants::PION; | |
422 | break; | |
423 | case starlightConstants::F2PRIME: | |
424 | if( ytest >= 0.5 ) | |
425 | { | |
426 | //Decays 50/50 into k+ k- or k_s k_l | |
427 | ipid=starlightConstants::KAONCHARGE; | |
428 | } | |
429 | else | |
430 | { | |
431 | ipid=starlightConstants::KAONNEUTRAL; | |
432 | } | |
433 | break; | |
434 | default: | |
435 | cout<<"Rethink the daughter particles"<<endl; | |
436 | } | |
437 | } | |
438 | ||
439 | ||
440 | //______________________________________________________________________________ | |
441 | starlightConstants::event Gammagammasingle::produceEvent(int &/*ievent*/) | |
442 | { | |
443 | // Not in use anymore, default event struct returned | |
444 | return starlightConstants::event(); | |
445 | } | |
446 | ||
447 | ||
448 | //______________________________________________________________________________ | |
449 | // fix it ... lost functionality | |
450 | //starlightConstants::event Gammagammasingle::produceEvent(int &ievent) | |
451 | upcEvent Gammagammasingle::produceEvent() | |
452 | { | |
453 | // cout << "NOT IMPLEMENTED!" << endl; | |
454 | ||
455 | // return upcEvent(); | |
456 | ||
457 | // returns the vector with the decay particles inside. | |
458 | // onedecayparticle single; | |
459 | starlightConstants::event single; | |
460 | double comenergy = 0.; | |
461 | double rapidity = 0.; | |
462 | double parentE = 0.; | |
463 | double parentmomx=0.,parentmomy=0.,parentmomz=0.; | |
464 | ||
465 | //this function decays particles and writes events to a file | |
466 | //zeroing out the event structure | |
467 | single._numberOfTracks=0; | |
468 | for(int i=0;i<4;i++){ | |
469 | single.px[i]=0.; | |
470 | single.py[i]=0.; | |
471 | single.pz[i]=0.; | |
472 | single._fsParticle[i]=starlightConstants::UNKNOWN; | |
473 | single._charge[i]=0; | |
474 | } | |
475 | ||
476 | pickw(comenergy); | |
477 | picky(rapidity); | |
478 | parentMomentum(comenergy,rapidity,parentE,parentmomx,parentmomy,parentmomz); | |
479 | ||
480 | ||
481 | if(_GGsingInputpidtest != starlightConstants::F2 && _GGsingInputpidtest != starlightConstants::F2PRIME) | |
482 | { | |
483 | #ifdef ENABLE_PYTHIA | |
484 | starlightParticle particle(parentmomx,parentmomy,parentmomz, parentE, getMass(),_GGsingInputpidtest , 0); | |
485 | ||
486 | _pyDecayer.addParticle(particle); | |
487 | ||
488 | return _pyDecayer.execute(); | |
489 | #endif | |
490 | } | |
491 | ||
492 | ||
493 | int ievent = 0; | |
494 | int iFbadevent=0; | |
495 | starlightConstants::particleTypeEnum ipid = starlightConstants::UNKNOWN; | |
496 | double px2=0.,px1=0.,py2=0.,py1=0.,pz2=0.,pz1=0.; | |
497 | double px3=0.,px4=0.,py3=0.,py4=0.,pz3=0.,pz4=0.; | |
498 | // double theta=0.,phi=0.;//angles from jetset | |
499 | double xtest=0.,ztest=0.; | |
500 | switch(_GGsingInputpidtest){ | |
501 | case starlightConstants::ZOVERZ03: | |
502 | //Decays into two pairs. | |
503 | parentmomx=parentmomx/2.; | |
504 | parentmomy=parentmomy/2.; | |
505 | parentmomz=parentmomz/2.; | |
506 | //Pair #1 | |
507 | twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent); | |
508 | //Pair #2 | |
509 | twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px3,py3,pz3,px4,py4,pz4,iFbadevent); | |
510 | //Now add them to vectors to be written out later. | |
511 | ||
512 | single._numberOfTracks=4;//number of tracks per event | |
513 | if (iFbadevent==0){ | |
514 | xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
515 | ztest = randyInstance.Rndom(); | |
516 | //Assigning charges randomly. | |
517 | if (xtest<0.5){ | |
518 | single._charge[0]=1;//q1=1; | |
519 | single._charge[1]=-1;//q2=-1; | |
520 | } | |
521 | else{ | |
522 | single._charge[0]=1;//q1=-1; | |
523 | single._charge[1]=-1;//q2=1; | |
524 | } | |
525 | if (ztest<0.5){ | |
526 | single._charge[2]=1;//q3=1; | |
527 | single._charge[3]=-1;//q4=-1; | |
528 | } | |
529 | else{ | |
530 | single._charge[2]=-1;//q3=-1; | |
531 | single._charge[3]=1;//q4=1; | |
532 | } | |
533 | //Track #1 | |
534 | single.px[0]=px1; | |
535 | single.py[0]=py1; | |
536 | single.pz[0]=pz1; | |
537 | single._fsParticle[0]=ipid; | |
538 | //Track #2 | |
539 | single.px[1]=px2; | |
540 | single.py[1]=py2; | |
541 | single.pz[1]=pz2; | |
542 | single._fsParticle[1]=ipid; | |
543 | //Track #3 | |
544 | single.px[2]=px3; | |
545 | single.py[2]=py3; | |
546 | single.pz[2]=pz3; | |
547 | single._fsParticle[2]=ipid; | |
548 | //Track #4 | |
549 | single.px[3]=px4; | |
550 | single.py[3]=py4; | |
551 | single.pz[3]=pz4; | |
552 | single._fsParticle[3]=ipid; | |
553 | ||
554 | ievent=ievent+1; | |
555 | } | |
556 | ||
557 | break; | |
558 | case starlightConstants::F2: | |
559 | case starlightConstants::F2PRIME: | |
560 | twoBodyDecay(ipid,parentE,comenergy,parentmomx,parentmomy,parentmomz,px1,py1,pz1,px2,py2,pz2,iFbadevent); | |
561 | ||
562 | single._numberOfTracks=2; | |
563 | if (iFbadevent==0){ | |
564 | xtest = randyInstance.Rndom();//random()/(RAND_MAX+1.0); | |
565 | if (xtest<0.5){ | |
566 | single._charge[0]=1;//q1=1; | |
567 | single._charge[1]=-1;//q2=-1; | |
568 | } | |
569 | else{ | |
570 | single._charge[0]=-1;//q1=-1; | |
571 | single._charge[1]=1;//q2=1; | |
572 | } | |
573 | //Track #1 | |
574 | single.px[0]=px1; | |
575 | single.py[0]=py1; | |
576 | single.pz[0]=pz1; | |
577 | single._fsParticle[0]=ipid*single._charge[0]; | |
578 | //Track #2 | |
579 | single.px[1]=px2; | |
580 | single.py[1]=py2; | |
581 | single.pz[1]=pz2; | |
582 | single._fsParticle[1]=ipid*single._charge[1]; | |
583 | ievent=ievent+1; | |
584 | } | |
585 | break; | |
586 | default: | |
587 | break; | |
588 | } | |
589 | ||
590 | return upcEvent(single); | |
591 | } | |
592 | ||
593 | ||
594 | //______________________________________________________________________________ | |
595 | void Gammagammasingle::thephi(double W,double px,double py,double pz,double E,double &theta,double &phi) | |
596 | { | |
597 | // This subroutine calculates angles for channels decayed by jetset. | |
598 | // subroutine thephi(W,px,py,pz,E,theta,phi) | |
599 | E = sqrt (W*W+px*px+py*py+pz*pz); | |
600 | ||
601 | theta = acos(pz/sqrt(px*px+py*py+pz*pz)); | |
602 | phi = acos(px/sqrt(px*px+py*py)); | |
603 | ||
604 | if ((px == 0) && (py == 0)) | |
605 | phi = 0.; | |
606 | if (py < 0) | |
607 | phi = 2*starlightConstants::pi - phi; | |
608 | } | |
609 | ||
610 | ||
611 | //______________________________________________________________________________ | |
612 | double Gammagammasingle::getMass() | |
613 | { | |
614 | using namespace starlightConstants; | |
615 | double singlemass=0.; | |
616 | switch(_GGsingInputpidtest){ | |
617 | case starlightConstants::ETA: | |
618 | singlemass= etaMass; | |
619 | break; | |
620 | case starlightConstants::ETAPRIME: | |
621 | singlemass=etaPrimeMass; | |
622 | break; | |
623 | case starlightConstants::ETAC: | |
624 | singlemass=etaCMass; | |
625 | break; | |
626 | case starlightConstants::F0: | |
627 | singlemass=f0Mass; | |
628 | break; | |
629 | case starlightConstants::F2: | |
630 | singlemass=f2Mass; | |
631 | break; | |
632 | case starlightConstants::A2: | |
633 | singlemass=a2Mass; | |
634 | break; | |
635 | case starlightConstants::F2PRIME: | |
636 | singlemass=f2PrimeMass; | |
637 | break; | |
638 | case starlightConstants::ZOVERZ03: | |
639 | singlemass=1.540; | |
640 | break; | |
641 | default: | |
642 | cout<<"Not a recognized single particle, Gammagammasingle::getmass(), mass = 0."<<endl; | |
643 | } | |
644 | return singlemass; | |
645 | } | |
646 | ||
647 | ||
648 | //______________________________________________________________________________ | |
649 | double Gammagammasingle::getWidth() | |
650 | { | |
651 | double singlewidth=0.; | |
652 | switch(_GGsingInputpidtest){ | |
653 | case starlightConstants::ETA: | |
654 | singlewidth=1.E-6; | |
655 | break; | |
656 | case starlightConstants::ETAPRIME: | |
657 | singlewidth=5.E-6; | |
658 | break; | |
659 | case starlightConstants::ETAC: | |
660 | singlewidth=6.4E-6; | |
661 | break; | |
662 | case starlightConstants::F0: | |
663 | singlewidth=0.56E-6; | |
664 | break; | |
665 | case starlightConstants::F2: | |
666 | singlewidth=2.6E-6; | |
667 | break; | |
668 | case starlightConstants::A2: | |
669 | singlewidth=1.04E-6; | |
670 | break; | |
671 | case starlightConstants::F2PRIME: | |
672 | singlewidth=0.1E-6; | |
673 | break; | |
674 | case starlightConstants::ZOVERZ03: | |
675 | singlewidth=0.1E-6; | |
676 | break; | |
677 | default: | |
678 | cout<<"Not a recognized single particle, Gammagammasingle::getwidth(), width = 0."<<endl; | |
679 | } | |
680 | return singlewidth; | |
681 | } | |
682 | ||
683 | ||
684 | //______________________________________________________________________________ | |
685 | double Gammagammasingle::getSpin() | |
686 | { | |
687 | double singlespin=0.5; | |
688 | switch(_GGsingInputpidtest){ | |
689 | case starlightConstants::ETA: | |
690 | singlespin=0.0; | |
691 | break; | |
692 | case starlightConstants::ETAPRIME: | |
693 | singlespin=0.0; | |
694 | break; | |
695 | case starlightConstants::ETAC: | |
696 | singlespin=0.0; | |
697 | break; | |
698 | case starlightConstants::F0: | |
699 | singlespin=0.0; | |
700 | break; | |
701 | case starlightConstants::F2: | |
702 | singlespin=2.0; | |
703 | break; | |
704 | case starlightConstants::A2: | |
705 | singlespin=2.0; | |
706 | break; | |
707 | case starlightConstants::F2PRIME: | |
708 | singlespin=2.0; | |
709 | break; | |
710 | case starlightConstants::ZOVERZ03: | |
711 | singlespin=2.0; | |
712 | break; | |
713 | default: | |
714 | cout<<"Not a recognized single particle, Gammagammasingle::getspin(), spin = 0."<<endl; | |
715 | } | |
716 | return singlespin; | |
717 | } | |
718 | ||
719 | ||
720 |