Pluto v6.3.0
Loading...
Searching...
No Matches
PHUrDilep.cc
Go to the documentation of this file.
1
2// dileptons (D0,D+,w,rho0,phi,eta,eta' and pi0)
3// can be decayed via shining (translated for original UrQMD
4// fortran routines, authors: Katharina Schmidt, Sascha Vogel,
5// Christian Sturm 2007-2009) into e+/e-.
6//
7// USAGE: This classs is called via PHUrReader, see urqmd_f15_input.C
8//
10
11#include "PHUrDilep.h"
12#include "TMath.h"
13#include "PParticle.h"
14#include "PHUrReader.h"
15
16ClassImp(PHUrDilep)
17
18
19PHUrDilep::PHUrDilep() {
20 omega = 103;
21 rho = 104;
22 phi = 109;
23 pion = 101;
24 etaprime = 107;
25 eta = 102;
26 delta = 17;
27 omegadir = 1031;
28 omegadal = 1032;
29
30 alpha_em = (1./137.036);
31 vacmass_pi0 = 0.1349764; // pi0
32 mass_electron = 0.00051099906; // e-
33 mass_nucleon = 0.93827231; // proton
34 vacmass_eta = 0.54745; // eta
35 vacmass_etaprime = 0.9577; // eta'
36 vacmass_rho = 0.7699; // rho0
37 vacmass_omega = 0.78194; // omega
38 vacmass_phi = 1.019413; // phi
39 vacmass_delta = 1.2320; // D+
40
41 gev = 0.2;
42 g = 5.44; // delta
43
44 lambda_omega = 0.65;
45 gamma_omega = 0.04;
46 gamma_photon = 0.00075561;
47
48 br_pi0 = 0.988;
49 b_pi0 = 5.5;
50
51 lambda_eta = 0.72;
52 br_eta = 0.3933;
53
54 lambda_etaprime = 0.76;
55 gamma_etaprime = 0.1;
56 br_etaprime = 0.0212;
57
58 out = NULL;
59 outputLeptons = kTRUE;
60 reader = NULL;
61}
62
63PHUrDilep::~PHUrDilep() {
64 if(out) fclose(out);
65 out = 0;
66}
67
68void PHUrDilep::Output(TString infile,TString outfile) {
69 if(outfile.CompareTo("") != 0) {
70 out = fopen(outfile.Data(), "w");
71 fprintf(out, "# InputFile : %s\n", infile.Data());
72 fprintf(out, "# OutputFile : %s\n", outfile.Data());
73 fprintf(out, "# ityp weight mres p0_el_lab px_el_lab py_el_lab pz_el_lab p0_po_lab px_po_lab py_po_lab pz_po_lab tau time_start time_end flag dens_cre dens_abs noe bev acce accp\n");
74 }
75}
76
77void PHUrDilep::Dilep() {
78
79 Double_t tau = (particleIn->t - particleOut->t)/particleIn->Gamma();
80
81 // multi: multiplier for better statistics for rare resonances
82 Int_t multi = 1;
83
84 // data input
85
86 // tbeam_input is a char ... in GeV ... quick and dirty work around !!! Take care !!!
87 // calculate beta from input tbeam
88
89 Double_t ebeam = (evtheader->E_beam+mass_nucleon);
90 Double_t pbeam = sqrt(pow(ebeam,2) - pow(mass_nucleon,2));
91
92 Double_t beta = pbeam/(ebeam+mass_nucleon);
93 Double_t gamma = 1./sqrt(1.-beta*beta);
94
95 Double_t weight = 0; // changed in dgamma_sum()
96 Double_t smin, smax, tmax; // changed in dgamma_sum()
97 smin = smax = tmax = 0;
98
99 Double_t s, t; // changed in gamma_star()
100 Double_t tgamma; // changed in gamma_star()
101 Double_t xgstar, ygstar, zgstar, tgstar; // changed in gamma_star()
102 Double_t xgamma, ygamma, zgamma; // changed in gamma_star()
103
104 Double_t dgamma; // changed in dalpi(),daleta(),daletaprime(),dalomega()
105 Double_t mgstar; // mgstar=s s from gamma_star()
106 Double_t p0_gstar, p0_nucleon, p0_pion, p0_gamma; // local variables
107
108 // ccccc calculation of decay widths ccccc
109 dgamma_sum(particleIn->id, tau, multi,
110 weight, smin, smax, tmax); // output
111
112 if(!outputLeptons) {
113 if( particleIn->id == omega ||
114 (particleIn->id == rho && particleIn->chrg == 0) ||
115 particleIn->id == phi ||
116 (particleIn->id == pion && particleIn->chrg == 0) ||
117 particleIn->id == eta ||
118 particleIn->id == etaprime ||
119 (particleIn->id == delta && (particleIn->chrg == 0 || particleIn->chrg == 1) )
120 ) OutputROOT(particleIn->id,weight);
121 return;
122 }
123
124
125 if (particleIn->id == omega) { // decay of omega meson (direct and Dalitz)
126 Double_t weight_omegadir = 0;
127 for(Int_t n = 1; n <= multi; n++) {
128 diromega(tau, particleOut->mass, multi,
129 weight_omegadir); // output
130
131 lobo_dir(beta, gamma, particleOut->mass); // em,ep output
132
133 Output(omegadir, weight_omegadir, particleOut->mass, multi,tau);
134
135 } // for multi omega
136 Int_t validm = 0;
137
138 while(validm < multi){
139 t_omega(tau, multi,
140 smin, smax, tmax); // output
141
142 gamma_star(smin, smax, tmax,
143 s, t, // output
144 xgstar, ygstar, zgstar, tgstar, // output
145 xgamma, ygamma, zgamma, tgamma); // output
146
147 dalomega(tau, s,
148 dgamma); // output
149
150 if (t <= dgamma) {
151 mgstar = s;
152
153 // energy of the gamma* and the pion due to energy and
154 // momentum conservation
155
156 p0_gstar = particleOut->mass/2. -(pow(vacmass_pi0,2) - pow(mgstar,2))/(2.*particleOut->mass);
157 p0_pion = particleOut->mass/2. +(pow(vacmass_pi0,2) - pow(mgstar,2))/(2.*particleOut->mass);
158
159 // check if energy is larger than mass (only for omega, for
160 // low mass omegas it might happen that p0_gstar < m_gstar
161
162 if(p0_gstar > mgstar) {
163 validm++;
164
165 lobo_dal(p0_gstar, p0_pion,
166 mgstar, beta, gamma,
167 xgstar, ygstar, zgstar, tgstar, xgamma, ygamma, zgamma, particleOut->mass);// em,ep output
168
169 Output(omegadal, weight, mgstar, multi, tau);
170
171 } // if p0
172 }
173 } // while vaidm < multi
174
175 //----------------------------------------------------------------------
176 } else if (particleIn->id == rho && particleIn->chrg == 0) { // direct decay of the rho meson
177
178 for(Int_t n = 1; n <= multi; n++){
179 dirrho(tau, particleOut->mass, multi,
180 weight); // output
181
182 lobo_dir(beta, gamma, particleOut->mass);// em,ep output
183
184 Output(rho, weight, particleOut->mass, multi, tau);
185
186 } // for multi rho
187
188 //----------------------------------------------------------------------
189 } else if (particleIn->id == phi) { // direct decay of the phi meson
190
191 for(Int_t n = 1; n <= multi; n++){
192 dirphi(tau, particleOut->mass,multi,
193 weight); // output
194
195 lobo_dir(beta, gamma, particleOut->mass);// em,ep output
196
197 Output(phi, weight, particleOut->mass, multi, tau);
198
199 } // for multi phi
200
201 //----------------------------------------------------------------------
202 } else if (particleIn->id == pion && particleIn->chrg == 0) { // Dalitz decay of the pi0 meson
203
204 Int_t validm = 0;
205
206 while(validm < multi) {
207 gamma_star(smin, smax, tmax,
208 s, t, // output
209 xgstar, ygstar, zgstar, tgstar, // output
210 xgamma, ygamma, zgamma, tgamma);// output
211
212 dalpi(s,
213 dgamma); // output
214
215 if (t <= dgamma) {
216 Double_t mass_gamma = 0;
217 mgstar = s;
218
219 p0_gstar = particleOut->mass/2 - (pow(mass_gamma,2) - pow(mgstar,2))/(2*particleOut->mass);
220 p0_gamma = particleOut->mass/2 + (pow(mass_gamma,2) - pow(mgstar,2))/(2*particleOut->mass);
221
222 validm = validm+1;
223
224 lobo_dal(p0_gstar,p0_gamma,
225 mgstar,beta,gamma,
226 xgstar,ygstar,zgstar,tgstar,xgamma,ygamma,zgamma,particleOut->mass);// em,ep output
227
228 Output (pion,weight,mgstar,multi,tau);
229 } // if t <= dgamma
230 } // end while validm < multi
231
232 //----------------------------------------------------------------------
233 } else if (particleIn->id == etaprime) { // Dalitz decay of the eta prime meson
234
235 Int_t validm = 0;
236
237 while (validm < multi) {
238 gamma_star(smin, smax, tmax,
239 s, t, // output
240 xgstar, ygstar, zgstar, tgstar, // output
241 xgamma, ygamma, zgamma, tgamma); // output
242
243 daletaprime(s,
244 dgamma); // output
245
246 if (t <= dgamma) {
247 Double_t mass_gamma = 0;
248 mgstar = s;
249
250 p0_gstar = particleOut->mass/2. - (pow(mass_gamma,2) - pow(mgstar,2))/(2.*particleOut->mass);
251 p0_gamma = particleOut->mass/2. + (pow(mass_gamma,2) - pow(mgstar,2))/(2.*particleOut->mass);
252
253 validm++;
254
255 lobo_dal(p0_gstar, p0_gamma,
256 mgstar, beta, gamma,
257 xgstar, ygstar, zgstar, tgstar, xgamma, ygamma, zgamma, particleOut->mass);// em,ep output
258
259 Output(etaprime, weight, mgstar, multi, tau);
260
261 } // if t <= dgamma
262 } // end while validm < multi
263
264
265 //----------------------------------------------------------------------
266 } else if (particleIn->id == eta) { // Dalitz decay of the eta meson
267 Int_t validm = 0;
268
269 while(validm < multi) {
270 gamma_star(smin, smax, tmax,
271 s, t, // output
272 xgstar, ygstar, zgstar, tgstar, // output
273 xgamma, ygamma, zgamma, tgamma); // output
274
275 daleta(s,
276 dgamma); // output
277
278 if (t <= dgamma){
279 Double_t mass_gamma = 0;
280 mgstar = s;
281
282 p0_gstar = particleOut->mass/2. - (pow(mass_gamma,2) - pow(mgstar,2))/(2.*particleOut->mass);
283 p0_gamma = particleOut->mass/2. + (pow(mass_gamma,2) - pow(mgstar,2))/(2.*particleOut->mass);
284
285 validm++;
286
287 lobo_dal(p0_gstar, p0_gamma, mgstar, beta, gamma,
288 xgstar, ygstar, zgstar, tgstar, xgamma, ygamma, zgamma, particleOut->mass);// em,ep output
289
290 Output (eta, weight, mgstar, multi, tau);
291
292 } // if t <= dgamma
293 } // end validm < multi
294
295 //----------------------------------------------------------------------
296 } else if((particleIn->id == delta && particleIn->chrg == 0) ||
297 (particleIn->id == delta && particleIn->chrg == 1)) { // Dalitz decay of the Delta(1232) baryon
298
299 Int_t validm = 0;
300
301 dgamma_sum_delta(tau, particleOut->mass, multi,
302 weight, // output
303 smin, smax, tmax); // output
304
305 while(validm < multi) {
306
307 t_delta(tau, particleOut->mass, multi,
308 tmax); // output
309
310 gamma_star(smin, smax, tmax,
311 s, t, // output
312 xgstar, ygstar, zgstar, tgstar, // output
313 xgamma, ygamma, zgamma, tgamma); // output
314
315 daldelta(tau, s, particleOut->mass,
316 dgamma); // output
317
318
319 if (t <= dgamma) {
320
321 mgstar = s ;
322
323 p0_gstar = particleOut->mass/2. - (pow(mass_nucleon,2) - pow(mgstar,2))/(2.*particleOut->mass);
324 p0_nucleon = particleOut->mass/2. + (pow(mass_nucleon,2) - pow(mgstar,2))/(2.*particleOut->mass);
325
326 validm++;
327
328 lobo_dal(p0_gstar, p0_nucleon,
329 mgstar, beta, gamma,
330 xgstar, ygstar, zgstar, tgstar, xgamma, ygamma, zgamma, particleOut->mass);// em,ep output
331
332 Output(delta, weight, mgstar, multi, tau);
333
334 } // if t <= dgamma
335 } // end validm < multi
336 } // end if case
337}
338
339void PHUrDilep::dalpi(Double_t mx, Double_t &dgamma) {
340 // form factor from Landsberg 1985 (Phys.Rept.128,301-376)
341
342 Double_t t1 = ((4.*alpha_em)/(3.*TMath::Pi()*mx))*
343 sqrt(1.-(4.*pow(mass_electron,2))/(pow(mx,2)))*
344 (1.+(2.*pow(mass_electron,2))/(pow(mx,2)));
345 Double_t t2 = pow((1.-(pow(mx,2))/(pow(vacmass_pi0,2))),3);
346 Double_t f1 = pow((1.+b_pi0*pow(mx,2)),2);
347 dgamma = t1*t2*f1*br_pi0;
348}
349
350void PHUrDilep::daleta(Double_t mx, Double_t &dgamma) {
351
352 //### integration dgamma, for more details please refer to (REFERENCE, MANUAL) ####
353
354 Double_t t1 =((4.*alpha_em)/(3.*TMath::Pi()*mx))*
355 sqrt(1.-(4.*pow(mass_electron,2))/(pow(mx,2)))*(1.+
356 (2.*pow(mass_electron,2))/(pow(mx,2)));
357 Double_t t2 = pow((1.-(pow(mx,2))/(pow(vacmass_eta,2))),3);
358
359 //### form factor from Phys.Rept.308:65-233,1999 with a modification
360 // from the diploma thesis of C.Ernst to ensure that F(0)=1 ###
361 Double_t f1 = pow((1./(1.-(pow(mx,2))/(pow(lambda_eta,2)))),2);
362
363 dgamma = t1*t2*f1*br_eta;
364}
365
366void PHUrDilep::daletaprime(Double_t mx, Double_t &dgamma) {
367
368 //#### integration dgamma, for more details please refer to (REFERENCE, MANUAL) ####
369
370 // ######### Landsberg 85
371 Double_t t1 = ((4.*alpha_em)/(3.*TMath::Pi()*mx))
372 *sqrt(1.-(4.*pow(mass_electron,2))/(pow(mx,2)))*(1.+
373 (2.*pow(mass_electron,2))/(pow(mx,2)));
374 Double_t t2 = pow((1.-(pow(mx,2))/(pow(vacmass_etaprime,2))),3);
375
376 //### form factor from Phys.Rept.308:65-233,1999 with a modification from the diploma thesis of C.Ernst to ensure that F(0)=1 ###
377
378 Double_t f1 = (pow(lambda_etaprime,2)*
379 (pow(lambda_etaprime,2)+pow(gamma_etaprime,2)))/
380 (pow((pow(lambda_etaprime,2)-pow(mx,2)),2)+
381 (pow(lambda_etaprime,2)*pow(gamma_etaprime,2)));
382
383 dgamma = t1*t2*f1*br_etaprime;
384}
385
386void PHUrDilep::daldelta(Double_t tau, Double_t mx, Double_t mres, Double_t &dgamma) {
387 //
388 // masses are in Gev/c2
389
390 // ######### Landsberg 85
391
392 Double_t pf = sqrt((pow(mres,2)-pow((mass_nucleon+mx),2))*
393 (pow(mres,2)-pow((mass_nucleon-mx),2)))/2./mres;
394
395 Double_t q0 = sqrt(pow(mx,2)+pow(pf,2));
396
397
398 Double_t f = (-3./2.)*(mres+mass_nucleon)/(mass_nucleon*
399 (pow((mres+mass_nucleon),2)-pow(mx,2)));
400
401 Double_t e = sqrt(4.*TMath::Pi()*alpha_em);
402
403 Double_t mt = pow(e*f*g,2)*pow(mres,2)/9./mass_nucleon *(q0*q0*(5.*mres-3.*
404 (q0+mass_nucleon))-pow(mx,2)*(mres+mass_nucleon+q0));
405
406 Double_t ml = pow(e*f*g,2)*pow(mres,2)/9./mass_nucleon*pow(mx,2)*4.*
407 (mres-mass_nucleon-q0);
408
409 Double_t lambda = pow(mx,4)+pow(mass_nucleon,4)+pow(mres,4)-2.*
410 (pow(mx,2)*pow(mass_nucleon,2)+pow(mx,2)*pow(mres,2)+pow(mass_nucleon,2)*pow(mres,2));
411
412 Double_t gamma0 = sqrt(lambda)/(16.*TMath::Pi()*pow(mres,2))*mass_nucleon*(2.*mt+ml);
413 dgamma = (2.*alpha_em)/(3.*TMath::Pi()*mx)*gamma0*(tau/gev);
414
415}
416
417void PHUrDilep::dalomega(Double_t tau, Double_t mx, Double_t &dgamma) {
418
419 //######### Landsberg 85
420
421 Double_t t1 = ((2.*alpha_em)/(3.*TMath::Pi()*mx))*
422 sqrt(1.-(4.*pow(mass_electron,2))/(pow(mx,2)))*(1.+
423 (2.*pow(mass_electron,2))/(pow(mx,2)));
424
425 Double_t t2 = sqrt(pow((1.+(pow(mx,2))/(pow(vacmass_omega,2)-pow(vacmass_pi0,2))),2)-
426 pow(((2.*vacmass_omega*mx)/(pow(vacmass_omega,2)-pow(vacmass_pi0,2))),2));
427
428 Double_t f1 =(pow(lambda_omega,2)*(pow(lambda_omega,2)+pow(gamma_omega,2)))/
429 (pow((pow(lambda_omega,2)-pow(mx,2)),2)+(pow(lambda_omega,2)*pow(gamma_omega,2)));
430
431 dgamma = t1*pow(t2,3)*f1*gamma_photon*(tau/gev);
432}
433
434void PHUrDilep::diromega(Double_t tau, Double_t mres, Int_t multi, Double_t &weight) {
435
436 const Double_t cv = 0.000000778869;
437
438 Double_t dwidth = cv/pow(mres,3)*pow(vacmass_omega,4);
439 Double_t br = dwidth*(tau/gev);
440
441 weight = br/multi;
442}
443
444void PHUrDilep::dirphi(Double_t, Double_t mres, Int_t multi, Double_t &weight) {
445
446 const Double_t gamma_tot = 0.00426;
447 const Double_t cv = 0.0000012411;
448
449 Double_t dwidth = cv/pow(mres,3)*pow(vacmass_phi,4);
450 Double_t br = dwidth/gamma_tot;
451
452 weight = br/multi;
453}
454
455void PHUrDilep::dirrho(Double_t tau, Double_t mres, Int_t multi, Double_t &weight) {
456
457 const Double_t cv = 0.0000090545;
458
459 Double_t dwidth = (cv/ pow(mres,3))*pow(vacmass_rho,4);
460 Double_t br = dwidth*(tau/gev);
461
462 weight = br/multi;
463}
464
465void PHUrDilep::lobo_dal(Double_t p0_gstar, Double_t,
466 Double_t m_gstar, Double_t beta, Double_t gamma,
467 Double_t x_gstar, Double_t y_gstar, Double_t z_gstar, Double_t,
468 Double_t, Double_t, Double_t,
469 Double_t) {
470 Double_t p0_res = particleOut->E();
471 Double_t px_res = particleOut->Px();
472 Double_t py_res = particleOut->Py();
473 Double_t pz_res = particleOut->Pz();
474
475// #### energies of the gamma* and the particle follow from energy and momentan conservation in the decay process ####
476
477 Double_t p_gstar = sqrt(p0_gstar*p0_gstar - m_gstar*m_gstar);
478 Double_t px_gstar = x_gstar*p_gstar;
479 Double_t py_gstar = y_gstar*p_gstar;
480 Double_t pz_gstar = z_gstar*p_gstar;
481
482 Double_t phi_el = 2.*TMath::Pi()*rndfunc();
483 Double_t x = rndfunc();
484 Double_t sin_theta_el = (x-0.5)*2.;
485 Double_t theta_el = asin(sin_theta_el)+TMath::PiOver2();
486
487 Double_t phi_po = phi_el + TMath::Pi() ;
488
489 if(phi_po > 2*TMath::Pi()){
490 phi_po = phi_po-2*TMath::Pi();
491 }
492 Double_t theta_po = 0;
493 if(theta_el <= TMath::Pi()){
494 theta_po = TMath::Pi()-theta_el;
495 } else {
496 theta_po = 0;
497 }
498
499 Double_t x_el = sin(theta_el)*cos(phi_el);
500 Double_t y_el = sin(theta_el)*sin(phi_el);
501 Double_t z_el = cos(theta_el);
502
503 Double_t x_po = sin(theta_po)*cos(phi_po);
504 Double_t y_po = sin(theta_po)*sin(phi_po);
505 Double_t z_po = cos(theta_po);
506
507
508 Double_t p0_el = m_gstar/2.;
509 Double_t p0_po = m_gstar/2.;
510
511 Double_t p_el = sqrt(p0_el*p0_el-mass_electron*mass_electron);
512 Double_t p_po = sqrt(p0_po*p0_po-mass_electron*mass_electron);
513
514
515 Double_t px_el = x_el*p_el;
516 Double_t py_el = y_el*p_el;
517 Double_t pz_el = z_el*p_el;
518
519 Double_t px_po = x_po*p_po;
520 Double_t py_po = y_po*p_po;
521 Double_t pz_po = z_po*p_po;
522
523
524 Double_t beta_x_gstar = -px_gstar/p0_gstar;
525 Double_t beta_y_gstar = -py_gstar/p0_gstar;
526 Double_t beta_z_gstar = -pz_gstar/p0_gstar;
527
528 Double_t beta_gstar = sqrt(beta_x_gstar*beta_x_gstar+beta_y_gstar*beta_y_gstar+beta_z_gstar*beta_z_gstar);
529 Double_t gamma_gstar = 1./sqrt(1-beta_gstar*beta_gstar);
530
531 Double_t p0_el_res = gamma_gstar*p0_el - beta_x_gstar*gamma_gstar*px_el -
532 beta_y_gstar*gamma_gstar*py_el -beta_z_gstar*gamma_gstar*pz_el;
533
534 Double_t px_el_res = (-beta_x_gstar*gamma_gstar)*p0_el +
535 (1+(gamma_gstar-1)*((beta_x_gstar*beta_x_gstar)/
536 (beta_gstar*beta_gstar)))*px_el +
537 (gamma_gstar-1)*((beta_x_gstar*beta_y_gstar)/
538 (beta_gstar*beta_gstar))*py_el +
539 (gamma_gstar-1)*((beta_x_gstar*beta_z_gstar)/
540 (beta_gstar*beta_gstar))*pz_el;
541
542 Double_t py_el_res = (-beta_y_gstar*gamma_gstar)*p0_el +
543 (gamma_gstar-1)*((beta_x_gstar*beta_y_gstar)/
544 (beta_gstar*beta_gstar))*px_el +
545 (1+(gamma_gstar-1)*((beta_y_gstar*beta_y_gstar)/
546 (beta_gstar*beta_gstar)))*py_el +
547 (gamma_gstar-1)*((beta_y_gstar*beta_z_gstar)/
548 (beta_gstar*beta_gstar))*pz_el;
549
550 Double_t pz_el_res = (-beta_z_gstar*gamma_gstar)*p0_el +
551 (gamma_gstar-1)*((beta_x_gstar*beta_z_gstar)/
552 (beta_gstar*beta_gstar))*px_el +
553 (gamma_gstar-1)*((beta_y_gstar*beta_z_gstar)/
554 (beta_gstar*beta_gstar))*py_el +
555 (1+(gamma_gstar-1)*((beta_z_gstar*beta_z_gstar)/
556 (beta_gstar*beta_gstar)))*pz_el;
557
558 Double_t p0_po_res = gamma_gstar*p0_po - beta_x_gstar*gamma_gstar*px_po -
559 beta_y_gstar*gamma_gstar*py_po -beta_z_gstar*gamma_gstar*pz_po;
560
561 Double_t px_po_res = (-beta_x_gstar*gamma_gstar)*p0_po +
562 (1+(gamma_gstar-1)*((beta_x_gstar*beta_x_gstar)/
563 (beta_gstar*beta_gstar)))*px_po +
564 (gamma_gstar-1)*((beta_x_gstar*beta_y_gstar)/
565 (beta_gstar*beta_gstar))*py_po +
566 (gamma_gstar-1)*((beta_x_gstar*beta_z_gstar)/
567 (beta_gstar*beta_gstar))*pz_po;
568
569 Double_t py_po_res = (-beta_y_gstar*gamma_gstar)*p0_po +
570 (gamma_gstar-1)*((beta_x_gstar*beta_y_gstar)/
571 (beta_gstar*beta_gstar))*px_po +
572 (1+(gamma_gstar-1)*((beta_y_gstar*beta_y_gstar)/
573 (beta_gstar*beta_gstar)))*py_po +
574 (gamma_gstar-1)*((beta_y_gstar*beta_z_gstar)/
575 (beta_gstar*beta_gstar))*pz_po;
576
577 Double_t pz_po_res = (-beta_z_gstar*gamma_gstar)*p0_po +
578 (gamma_gstar-1)*((beta_x_gstar*beta_z_gstar)/
579 (beta_gstar*beta_gstar))*px_po +
580 (gamma_gstar-1)*((beta_y_gstar*beta_z_gstar)/
581 (beta_gstar*beta_gstar))*py_po +
582 (1+(gamma_gstar-1)*((beta_z_gstar*beta_z_gstar)/
583 (beta_gstar*beta_gstar)))*pz_po;
584
585
586
587 Double_t beta_x_res = -px_res/p0_res;
588 Double_t beta_y_res = -py_res/p0_res;
589 Double_t beta_z_res = -pz_res/p0_res;
590
591 Double_t beta_res = sqrt(beta_x_res*beta_x_res+beta_y_res*beta_y_res+beta_z_res*beta_z_res);
592
593 Double_t gamma_res = 1./sqrt(1-beta_res*beta_res);
594
595 Double_t p0_el_cm = gamma_res*p0_el_res-beta_x_res*gamma_res*px_el_res-
596 beta_y_res*gamma_res*py_el_res-beta_z_res*gamma_res*pz_el_res;
597
598 Double_t px_el_cm = (-beta_x_res*gamma_res)*p0_el_res+
599 (1+(gamma_res-1)*((beta_x_res*beta_x_res)/
600 (beta_res*beta_res)))*px_el_res +
601 (gamma_res-1)*((beta_x_res*beta_y_res)/
602 (beta_res*beta_res))*py_el_res +
603 (gamma_res-1)*((beta_x_res*beta_z_res)/
604 (beta_res*beta_res))*pz_el_res;
605
606 Double_t py_el_cm = (-beta_y_res*gamma_res)*p0_el_res+
607 (gamma_res-1)*((beta_x_res*beta_y_res)/
608 (beta_res*beta_res))*px_el_res+
609 (1+(gamma_res-1)*((beta_y_res*beta_y_res)/
610 (beta_res*beta_res)))*py_el_res+
611 (gamma_res-1)*((beta_y_res*beta_z_res)/
612 (beta_res*beta_res))*pz_el_res;
613
614 Double_t pz_el_cm = (-beta_z_res*gamma_res)*p0_el_res+
615 (gamma_res-1)*((beta_x_res*beta_z_res)/
616 (beta_res*beta_res))*px_el_res+
617 (gamma_res-1)*((beta_y_res*beta_z_res)/
618 (beta_res*beta_res))*py_el_res+
619 (1+(gamma_res-1)*((beta_z_res*beta_z_res)/
620 (beta_res*beta_res)))*pz_el_res;
621
622
623 Double_t p0_po_cm = gamma_res*p0_po_res-beta_x_res*gamma_res*px_po_res-
624 beta_y_res*gamma_res*py_po_res -beta_z_res*gamma_res*pz_po_res;
625
626 Double_t px_po_cm = (-beta_x_res*gamma_res)*p0_po_res+
627 (1+(gamma_res-1)*((beta_x_res*beta_x_res)/
628 (beta_res*beta_res)))*px_po_res+
629 (gamma_res-1)*((beta_x_res*beta_y_res)/
630 (beta_res*beta_res))*py_po_res+
631 (gamma_res-1)*((beta_x_res*beta_z_res)/
632 (beta_res*beta_res))*pz_po_res;
633
634 Double_t py_po_cm = (-beta_y_res*gamma_res)*p0_po_res+
635 (gamma_res-1)*((beta_x_res*beta_y_res)/
636 (beta_res*beta_res))*px_po_res+
637 (1+(gamma_res-1)*((beta_y_res*beta_y_res)/
638 (beta_res*beta_res)))*py_po_res+
639 (gamma_res-1)*((beta_y_res*beta_z_res)/
640 (beta_res*beta_res))*pz_po_res;
641
642 Double_t pz_po_cm = (-beta_z_res*gamma_res)*p0_po_res+
643 (gamma_res-1)*((beta_x_res*beta_z_res)/
644 (beta_res*beta_res))*px_po_res+
645 (gamma_res-1)*((beta_y_res*beta_z_res)/
646 (beta_res*beta_res))*py_po_res+
647 (1+(gamma_res-1)*((beta_z_res*beta_z_res)/
648 (beta_res*beta_res)))*pz_po_res;
649
650
651 Double_t beta_x_urqmd = 0;
652 Double_t beta_y_urqmd = 0;
653 Double_t beta_z_urqmd = -beta;
654
655
656 Double_t beta_urqmd = beta;
657 Double_t gamma_urqmd = gamma;
658
659 // output
660
661 Double_t p0_el_lab = gamma_urqmd*p0_el_cm-beta_x_urqmd*gamma_urqmd*px_el_cm-
662 beta_y_urqmd*gamma_urqmd*py_el_cm-beta_z_urqmd*gamma_urqmd*
663 pz_el_cm;
664
665 Double_t px_el_lab = (-beta_x_urqmd*gamma)*p0_el_cm +
666 (1+(gamma_urqmd-1)*((beta_x_urqmd*beta_x_urqmd)/
667 (beta_urqmd*beta_urqmd)))*px_el_cm +
668 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
669 (beta_urqmd*beta_urqmd))*py_el_cm +
670 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
671 (beta_urqmd*beta_urqmd))*pz_el_cm;
672
673 Double_t py_el_lab = (-beta_y_urqmd*gamma)*p0_el_cm+
674 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
675 (beta_urqmd*beta_urqmd))*px_el_cm+
676 (1+(gamma_urqmd-1)*((beta_y_urqmd*beta_y_urqmd)/
677 (beta_urqmd*beta_urqmd)))*py_el_cm+
678 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
679 (beta_urqmd*beta_urqmd))*pz_el_cm;
680
681 Double_t pz_el_lab = (-beta_z_urqmd*gamma)*p0_el_cm+
682 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
683 (beta_urqmd*beta_urqmd))*px_el_cm+
684 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
685 (beta_urqmd*beta_urqmd))*py_el_cm+
686 (1+(gamma_urqmd-1)*((beta_z_urqmd*beta_z_urqmd)/
687 (beta_urqmd*beta_urqmd)))*pz_el_cm;
688
689 Double_t p0_po_lab = gamma_urqmd*p0_po_cm-beta_x_urqmd*gamma_urqmd*px_po_cm-
690 beta_y_urqmd*gamma_urqmd*py_po_cm -beta_z_urqmd*gamma_urqmd*
691 pz_po_cm;
692
693 Double_t px_po_lab = (-beta_x_urqmd*gamma)*p0_po_cm+
694 (1+(gamma_urqmd-1)*((beta_x_urqmd*beta_x_urqmd)/
695 (beta_urqmd*beta_urqmd)))*px_po_cm+
696 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
697 (beta_urqmd*beta_urqmd))*py_po_cm+
698 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
699 (beta_urqmd*beta_urqmd))*pz_po_cm;
700
701 Double_t py_po_lab = (-beta_y_urqmd*gamma)*p0_po_cm+
702 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
703 (beta_urqmd*beta_urqmd))*px_po_cm+
704 (1+(gamma_urqmd-1)*((beta_y_urqmd*beta_y_urqmd)/
705 (beta_urqmd*beta_urqmd)))*py_po_cm+
706 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
707 (beta_urqmd*beta_urqmd))*pz_po_cm;
708
709 Double_t pz_po_lab = (-beta_z_urqmd*gamma)*p0_po_cm+
710 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
711 (beta_urqmd*beta_urqmd))*px_po_cm+
712 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
713 (beta_urqmd*beta_urqmd))*py_po_cm+
714 (1+(gamma_urqmd-1)*((beta_z_urqmd*beta_z_urqmd)/
715 (beta_urqmd*beta_urqmd)))*pz_po_cm;
716
717 em.SetPxPyPzE(px_el_lab,py_el_lab,pz_el_lab,p0_el_lab);
718 ep.SetPxPyPzE(px_po_lab,py_po_lab,pz_po_lab,p0_po_lab);
719}
720
721void PHUrDilep::lobo_dir(Double_t beta, Double_t gamma, Double_t m_res) {
722
723 Double_t p0_res = particleOut->E();
724 Double_t px_res = particleOut->Px();
725 Double_t py_res = particleOut->Py();
726 Double_t pz_res = particleOut->Pz();
727
728 Double_t phi_el = 2.*TMath::Pi()*rndfunc();
729 Double_t x = rndfunc();
730 Double_t sin_theta_el = (x-0.5)*2.;
731 Double_t theta_el = asin(sin_theta_el)+TMath::PiOver2();
732
733 Double_t phi_po = phi_el + TMath::Pi();
734 if(phi_po > 2*TMath::Pi()) {
735 phi_po = phi_po-2*TMath::Pi();
736 }
737 Double_t theta_po = 0;
738 if(theta_el <= TMath::Pi()) {
739 theta_po = TMath::Pi()-theta_el;
740 } else {
741 theta_po = 0;
742 }
743
744 Double_t x_el = sin(theta_el)*cos(phi_el);
745 Double_t y_el = sin(theta_el)*sin(phi_el);
746 Double_t z_el = cos(theta_el);
747
748 Double_t x_po = sin(theta_po)*cos(phi_po);
749 Double_t y_po = sin(theta_po)*sin(phi_po);
750 Double_t z_po = cos(theta_po);
751
752
753 Double_t p0_el = m_res/2.;
754 Double_t p0_po = m_res/2.;
755
756 Double_t p_el = sqrt(pow(p0_el,2)-pow(mass_electron,2));
757 Double_t p_po = sqrt(pow(p0_po,2)-pow(mass_electron,2));
758
759 Double_t px_el = x_el*p_el;
760 Double_t py_el = y_el*p_el;
761 Double_t pz_el = z_el*p_el;
762
763 Double_t px_po = x_po*p_po;
764 Double_t py_po = y_po*p_po;
765 Double_t pz_po = z_po*p_po;
766
767
768 Double_t beta_x_res = -px_res/p0_res;
769 Double_t beta_y_res = -py_res/p0_res;
770 Double_t beta_z_res = -pz_res/p0_res;
771
772 Double_t beta_res = sqrt(beta_x_res*beta_x_res+beta_y_res*beta_y_res+beta_z_res*beta_z_res);
773
774 Double_t gamma_res = 1./sqrt(1-beta_res*beta_res);
775
776 Double_t p0_el_cm = gamma_res*p0_el-beta_x_res*gamma_res*px_el-
777 beta_y_res*gamma_res*py_el-beta_z_res*gamma_res*pz_el;
778
779 Double_t px_el_cm = (-beta_x_res*gamma_res)*p0_el+
780 (1+(gamma_res-1)*((beta_x_res*beta_x_res)/
781 (beta_res*beta_res)))*px_el+
782 (gamma_res-1)*((beta_x_res*beta_y_res)/
783 (beta_res*beta_res))*py_el+
784 (gamma_res-1)*((beta_x_res*beta_z_res)/
785 (beta_res*beta_res))*pz_el;
786
787 Double_t py_el_cm = (-beta_y_res*gamma_res)*p0_el+
788 (gamma_res-1)*((beta_x_res*beta_y_res)/
789 (beta_res*beta_res))*px_el+
790 (1+(gamma_res-1)*((beta_y_res*beta_y_res)/
791 (beta_res*beta_res)))*py_el+
792 (gamma_res-1)*((beta_y_res*beta_z_res)/
793 (beta_res*beta_res))*pz_el;
794
795 Double_t pz_el_cm = (-beta_z_res*gamma_res)*p0_el+
796 (gamma_res-1)*((beta_x_res*beta_z_res)/
797 (beta_res*beta_res))*px_el+
798 (gamma_res-1)*((beta_y_res*beta_z_res)/
799 (beta_res*beta_res))*py_el+
800 (1+(gamma_res-1)*((beta_z_res*beta_z_res)/
801 (beta_res*beta_res)))*pz_el;
802
803 Double_t p0_po_cm = gamma_res*p0_po-beta_x_res*gamma_res*px_po-
804 beta_y_res*gamma_res*py_po -beta_z_res*gamma_res*pz_po;
805
806 Double_t px_po_cm = (-beta_x_res*gamma_res)*p0_po+
807 (1+(gamma_res-1)*((beta_x_res*beta_x_res)/
808 (beta_res*beta_res)))*px_po+
809 (gamma_res-1)*((beta_x_res*beta_y_res)/
810 (beta_res*beta_res))*py_po+
811 (gamma_res-1)*((beta_x_res*beta_z_res)/
812 (beta_res*beta_res))*pz_po;
813
814 Double_t py_po_cm = (-beta_y_res*gamma_res)*p0_po+
815 (gamma_res-1)*((beta_x_res*beta_y_res)/
816 (beta_res*beta_res))*px_po+
817 (1+(gamma_res-1)*((beta_y_res*beta_y_res)/
818 (beta_res*beta_res)))*py_po+
819 (gamma_res-1)*((beta_y_res*beta_z_res)/
820 (beta_res*beta_res))*pz_po;
821
822 Double_t pz_po_cm = (-beta_z_res*gamma_res)*p0_po+
823 (gamma_res-1)*((beta_x_res*beta_z_res)/
824 (beta_res*beta_res))*px_po+
825 (gamma_res-1)*((beta_y_res*beta_z_res)/
826 (beta_res*beta_res))*py_po+
827 (1+(gamma_res-1)*((beta_z_res*beta_z_res)/
828 (beta_res*beta_res)))*pz_po;
829
830
831 Double_t beta_x_urqmd = 0;
832 Double_t beta_y_urqmd = 0;
833 Double_t beta_z_urqmd = -beta;
834
835 Double_t beta_urqmd = beta;
836 Double_t gamma_urqmd = gamma;
837
838 // output
839
840 Double_t p0_el_lab = gamma_urqmd*p0_el_cm-beta_x_urqmd*gamma_urqmd*px_el_cm-
841 beta_y_urqmd*gamma_urqmd*py_el_cm-beta_z_urqmd*gamma_urqmd*
842 pz_el_cm;
843
844 Double_t px_el_lab = (-beta_x_urqmd*gamma)*p0_el_cm +
845 (1+(gamma_urqmd-1)*((beta_x_urqmd*beta_x_urqmd)/
846 (beta_urqmd*beta_urqmd)))*px_el_cm +
847 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
848 (beta_urqmd*beta_urqmd))*py_el_cm +
849 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
850 (beta_urqmd*beta_urqmd))*pz_el_cm;
851
852 Double_t py_el_lab = (-beta_y_urqmd*gamma)*p0_el_cm+
853 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
854 (beta_urqmd*beta_urqmd))*px_el_cm+
855 (1+(gamma_urqmd-1)*((beta_y_urqmd*beta_y_urqmd)/
856 (beta_urqmd*beta_urqmd)))*py_el_cm+
857 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
858 (beta_urqmd*beta_urqmd))*pz_el_cm;
859
860 Double_t pz_el_lab = (-beta_z_urqmd*gamma)*p0_el_cm+
861 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
862 (beta_urqmd*beta_urqmd))*px_el_cm+
863 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
864 (beta_urqmd*beta_urqmd))*py_el_cm+
865 (1+(gamma_urqmd-1)*((beta_z_urqmd*beta_z_urqmd)/
866 (beta_urqmd*beta_urqmd)))*pz_el_cm;
867
868 Double_t p0_po_lab = gamma_urqmd*p0_po_cm-beta_x_urqmd*gamma_urqmd*px_po_cm-
869 beta_y_urqmd*gamma_urqmd*py_po_cm -beta_z_urqmd*gamma_urqmd*
870 pz_po_cm;
871
872 Double_t px_po_lab = (-beta_x_urqmd*gamma)*p0_po_cm+
873 (1+(gamma_urqmd-1)*((beta_x_urqmd*beta_x_urqmd)/
874 (beta_urqmd*beta_urqmd)))*px_po_cm+
875 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
876 (beta_urqmd*beta_urqmd))*py_po_cm+
877 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
878 (beta_urqmd*beta_urqmd))*pz_po_cm;
879
880 Double_t py_po_lab = (-beta_y_urqmd*gamma)*p0_po_cm+
881 (gamma_urqmd-1)*((beta_x_urqmd*beta_y_urqmd)/
882 (beta_urqmd*beta_urqmd))*px_po_cm+
883 (1+(gamma_urqmd-1)*((beta_y_urqmd*beta_y_urqmd)/
884 (beta_urqmd*beta_urqmd)))*py_po_cm+
885 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
886 (beta_urqmd*beta_urqmd))*pz_po_cm;
887
888 Double_t pz_po_lab = (-beta_z_urqmd*gamma)*p0_po_cm+
889 (gamma_urqmd-1)*((beta_x_urqmd*beta_z_urqmd)/
890 (beta_urqmd*beta_urqmd))*px_po_cm+
891 (gamma_urqmd-1)*((beta_y_urqmd*beta_z_urqmd)/
892 (beta_urqmd*beta_urqmd))*py_po_cm+
893 (1+(gamma_urqmd-1)*((beta_z_urqmd*beta_z_urqmd)/
894 (beta_urqmd*beta_urqmd)))*pz_po_cm;
895
896 em.SetPxPyPzE(px_el_lab,py_el_lab,pz_el_lab,p0_el_lab);
897 ep.SetPxPyPzE(px_po_lab,py_po_lab,pz_po_lab,p0_po_lab);
898
899}
900
901void PHUrDilep::dgamma_sum(Int_t ityp, Double_t tau, Int_t multi,
902 Double_t &weight,
903 Double_t &smin, Double_t &smax, Double_t &tmax) {
904 const Double_t dx = 0.001;
905
906 Double_t dgamma = 0, dgamma_sum = 0;
907 Double_t mx = 0, t1 = 0, t2 = 0, f1 = 0;
908 Int_t z = 0;
909 if(ityp == pion) { // pi0
910
911 smin = 0.0013;
912 smax = 0.1348;
913 tmax = 1.9065;
914
915 z = (Int_t)((smax-smin)/dx);
916
917 for(Int_t i = 1; i <= z; i++) {
918 mx = smin+dx*i;
919
920 t1 = ((4*alpha_em)/(3.*TMath::Pi()*mx)) *
921 sqrt(1-(4*pow(mass_electron,2))/(pow(mx,2)))*(1+(2*pow(mass_electron,2))/(pow(mx,2)));
922
923 t2 = pow((1-(pow(mx,2))/(pow(vacmass_pi0,2))),3);
924 f1 = pow((1+b_pi0*pow(mx,2)),2);
925
926 dgamma = t1*t2*f1*br_pi0;
927 dgamma_sum = dgamma_sum + dgamma;
928 }
929
930 weight = dgamma_sum*dx/multi;
931
932 } else if(ityp == omega) { // omega
933
934 smin = 0.0013;
935 smax = 0.6446;
936 tmax = 0.086086;
937
938 z = (Int_t)((smax-smin)/dx);
939
940 for(Int_t i=1; i <= z; i++) {
941 mx = smin+dx*i;
942
943 t1 = ((2*alpha_em)/(3*TMath::Pi()*mx))*
944 sqrt(1-(4*pow(mass_electron,2))/
945 (pow(mx,2)))*(1+(2*pow(mass_electron,2))/(pow(mx,2)));
946
947 t2 = sqrt(pow((1+(pow(mx,2))/(pow(vacmass_omega,2)-
948 pow(vacmass_pi0,2))),2)
949 -pow(((2*vacmass_omega*mx)/(pow(vacmass_omega,2)-pow(vacmass_pi0,2))),2));
950
951 f1 = (pow(lambda_omega,2)*(pow(lambda_omega,2)+pow(gamma_omega,2)))/
952 (pow((pow(lambda_omega,2)-pow(mx,2)),2)+
953 (pow(lambda_omega,2)*pow(gamma_omega,2)));
954
955 dgamma = t1*pow(t2,3)*f1*gamma_photon*(tau/gev);
956 dgamma_sum = dgamma_sum + dgamma;
957 }
958
959 weight = dgamma_sum*dx/multi;
960
961 } else if(ityp == etaprime) { // eta'
962
963 smin = 0.0013;
964 smax = 0.958;
965 tmax = 0.04092;
966
967 z = (Int_t)((smax-smin)/dx);
968
969 for(Int_t i=1; i <= z; i++) {
970 mx = smin+dx*i;
971
972 t1 = ((4*alpha_em)/(3*TMath::Pi()*mx))*
973 sqrt(1-(4*pow(mass_electron,2))/
974 (pow(mx,2)))*(1+(2*pow(mass_electron,2))/(pow(mx,2)));
975
976 t2 = pow((1-(pow(mx,2))/(pow(vacmass_etaprime,2))),3);
977
978 f1 = (pow(lambda_etaprime,2)*
979 (pow(lambda_etaprime,2)+pow(gamma_etaprime,2)))/
980 (pow((pow(lambda_etaprime,2)-pow(mx,2)),2)+
981 (pow(lambda_etaprime,2)*pow(gamma_etaprime,2)));
982
983 dgamma = t1*t2*f1*br_etaprime;
984 dgamma_sum = dgamma_sum + dgamma;
985
986 }
987
988 weight=dgamma_sum*dx/multi;
989
990 } else if(ityp == eta) { // eta
991
992 smin = 0.0013;
993 smax = 0.547;
994 tmax = 0.757592;
995
996 z = (Int_t)((smax-smin)/dx);
997
998 for(Int_t i=1; i <= z;i++) {
999 mx = smin+dx*i;
1000
1001 t1 = ((4.*alpha_em)/(3.*TMath::Pi()*mx))*
1002 sqrt(1.-(4.*pow(mass_electron,2))/
1003 (pow(mx,2))) * (1.+(2.*pow(mass_electron,2))/(pow(mx,2)));
1004
1005 t2 = pow((1.-(pow(mx,2))/(pow(vacmass_eta,2))),3);
1006 f1 = pow((1./(1.-(pow(mx,2))/(pow(lambda_eta,2)))),2);
1007
1008 dgamma = t1*t2*f1*br_eta;
1009 dgamma_sum = dgamma_sum + dgamma;
1010
1011 }
1012
1013 weight = dgamma_sum*dx/multi;
1014 }
1015}
1016
1017void PHUrDilep::dgamma_sum_delta(Double_t tau,Double_t mres,Int_t multi,
1018 Double_t& weight_delta,
1019 Double_t& smin_del,Double_t& smax_del,Double_t& tmax_del) {
1020
1021 const Double_t dx = 0.000001;
1022
1023 // delta
1024
1025 Double_t mx, pf, q0, e, f, mt, ml, lambda, gamma0, dgamma, dgamma_sum=0;
1026 //Double_t q,qr;
1027 // cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1028
1029
1030 smin_del = 0.0013;
1031 tmax_del = 0.00903;
1032 smax_del = mres-mass_nucleon;
1033
1034
1035 Int_t z_del = (Int_t)((smax_del-smin_del)/dx);
1036
1037 for(Int_t i=1; i <= z_del; i++) {
1038
1039 mx = smin_del+dx*i;
1040
1041 pf = sqrt((pow(mres,2)-pow((mass_nucleon+mx),2))*
1042 (pow(mres,2)-pow((mass_nucleon-mx),2)))/2./mres;
1043
1044 q0 = sqrt(mx*mx+pf*pf);
1045
1046/* q = sqrt((pow(mres,2)-pow((mass_nucleon+vacmass_pi0),2))*
1047 (pow(mres,2)-pow((mass_nucleon-vacmass_pi0),2))/(4.*pow(mres,2)));
1048
1049 qr = sqrt((pow(vacmass_delta,2)-pow((mass_nucleon+vacmass_pi0),2))*
1050 (pow(vacmass_delta,2)-pow((mass_nucleon-vacmass_pi0),2))/
1051 (4.*pow(vacmass_delta,2)));
1052*/
1053 f = (-3./2.)*(mres+mass_nucleon)/
1054 (mass_nucleon*(pow((mres+mass_nucleon),2)-mx*mx));
1055
1056 e = sqrt(4.*TMath::Pi()*alpha_em);
1057
1058 mt = pow((e*f*g),2)*pow(mres,2)/9./mass_nucleon*
1059 (pow(q0,2)*(5.*mres-3.*(q0+mass_nucleon))-
1060 pow(mx,2)*(mres+mass_nucleon+q0));
1061
1062 ml = pow((e*f*g),2)*pow(mres,2)/9./mass_nucleon*mx*mx*4.*
1063 (mres-mass_nucleon-q0);
1064
1065 lambda = pow(mx,4)+pow(mass_nucleon,4)+pow(mres,4)-2.*
1066 (mx*mx*pow(mass_nucleon,2)+mx*mx*pow(mres,2)+pow(mass_nucleon,2)*mres*mres);
1067
1068
1069 if(fabs(lambda) <= 1e-15) {
1070 lambda = 0;
1071 }
1072
1073 if(mx <= q0) {
1074 gamma0 = sqrt(lambda)/(16.*TMath::Pi()*mres*mres)*mass_nucleon*(2.*mt+ml);
1075 dgamma = (2.*alpha_em)/(3.*TMath::Pi()*mx)*gamma0*(tau/gev);
1076 dgamma_sum = dgamma_sum + dgamma;
1077 } else {
1078 dgamma_sum = 0; // added
1079 }
1080 }
1081
1082 weight_delta = dgamma_sum*dx/multi;
1083}
1084
1085void PHUrDilep::gamma_star(Double_t smin, Double_t smax, Double_t tmax,
1086 Double_t& s, Double_t& t,
1087 Double_t& xgstar, Double_t& ygstar, Double_t& zgstar, Double_t& tgstar,
1088 Double_t& xparticle, Double_t& yparticle, Double_t& zparticle, Double_t& tparticle) {
1089
1090 Double_t phi_gstar = 2.*TMath::Pi()*rndfunc();
1091 Double_t y = rndfunc();
1092 Double_t sin_theta_gstar = (y-0.5)*2.;
1093 Double_t theta_gstar = asin(sin_theta_gstar)+TMath::PiOver2();
1094
1095 Double_t phi_particle = phi_gstar + TMath::Pi();
1096 if (phi_particle > 2*TMath::Pi()) {
1097 phi_particle = phi_particle-2*TMath::Pi();
1098 }
1099 Double_t theta_particle = 0;
1100 if (theta_gstar <= TMath::Pi()) {
1101 theta_particle = TMath::Pi()-theta_gstar;
1102 } else {
1103 theta_particle = 0;
1104 }
1105
1106 xgstar = sin(theta_gstar)*cos(phi_gstar);
1107 ygstar = sin(theta_gstar)*sin(phi_gstar);
1108 zgstar = cos(theta_gstar);
1109 tgstar = sin(theta_gstar);
1110
1111 xparticle = sin(theta_particle)*cos(phi_particle);
1112 yparticle = sin(theta_particle)*sin(phi_particle);
1113 zparticle = cos(theta_particle);
1114 tparticle = sin(theta_particle);
1115
1116
1117 t = rndfunc()*tmax;
1118 s = smin + rndfunc()*(smax-smin);
1119
1120}
1121
1122void PHUrDilep::t_delta(Double_t tau, Double_t mres, Int_t, Double_t &tmax_del) {
1123 const Double_t smin_del = 0.0013;
1124
1125 // cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1126
1127 Double_t mx = smin_del;
1128
1129 Double_t pf = sqrt((pow(mres,2)-pow((mass_nucleon+mx),2))*
1130 (pow(mres,2)-pow((mass_nucleon-mx),2)))/2./mres;
1131
1132 Double_t q0 = sqrt(mx*mx+pf*pf);
1133
1134
1135 Double_t f = (-3./2.)*(mres+mass_nucleon)/
1136 (mass_nucleon*(pow((mres+mass_nucleon),2)-pow(mx,2)));
1137
1138 Double_t e = sqrt(4.*TMath::Pi()*alpha_em);
1139
1140 Double_t mt = pow(e*f*g,2)*pow(mres,2)/9./mass_nucleon*
1141 (q0*q0*(5.*mres-3.*(q0+mass_nucleon))-
1142 pow(mx,2)*(mres+mass_nucleon+q0));
1143
1144 Double_t ml = pow(e*f*g,2)*pow(mres,2)/9./mass_nucleon*pow(mx,2)*4.*
1145 (mres-mass_nucleon-q0);
1146
1147 Double_t lambda = pow(mx,4)+pow(mass_nucleon,4)+pow(mres,4)-2.*
1148 (pow(mx,2)*pow(mass_nucleon,2)+pow(mx,2)*pow(mres,2)+pow(mass_nucleon,2)*pow(mres,2));
1149
1150 if(fabs(lambda) <= 1e-15) {
1151 lambda = 0;
1152 }
1153
1154 if(mx <= q0) {
1155 Double_t gamma0 = sqrt(lambda)/(16.*TMath::Pi()*pow(mres,2))*mass_nucleon*(2.*mt+ml);
1156 Double_t dgamma = (2.*alpha_em)/(3.*TMath::Pi()*mx)*gamma0*(tau/gev);
1157 tmax_del = dgamma;
1158 }
1159}
1160
1161void PHUrDilep::t_omega(Double_t tau,Int_t,
1162 Double_t& smin_omega,Double_t& smax_omega,Double_t& tmax_omega) {
1163
1164 // cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
1165
1166
1167 smin_omega = 0.0013;
1168 smax_omega = 0.6446;
1169
1170
1171 Double_t mx_omega = smin_omega;
1172 Double_t t1_omega = ((2.*alpha_em)/(3.*TMath::Pi()*mx_omega))*
1173 sqrt(1-(4.*pow(mass_electron,2))/
1174 (pow(mx_omega,2))) * (1+(2.*pow(mass_electron,2))/(pow(mx_omega,2)));
1175
1176 Double_t t2_omega = sqrt(pow((1+pow(mx_omega,2)/(pow(vacmass_omega,2)-pow(vacmass_pi0,2))),2)
1177 -pow(((2.*vacmass_omega*mx_omega)/(pow(vacmass_omega,2)-pow(vacmass_pi0,2))),2));
1178
1179 Double_t f1_omega = (pow(lambda_omega,2)*(pow(lambda_omega,2)+pow(gamma_omega,2)))/
1180 (pow((pow(lambda_omega,2)-pow(mx_omega,2)),2)+(pow(lambda_omega,2)*pow(gamma_omega,2)));
1181
1182 Double_t dgamma_omega = t1_omega*pow(t2_omega,3)*f1_omega*gamma_photon*(tau/gev);
1183
1184 tmax_omega = dgamma_omega;
1185
1186}
1187
1188void PHUrDilep::Output(Int_t ityp, Double_t weight, Double_t mres, Int_t, Double_t tau) {
1189
1190 /* . ,
1191 write(*,'(I5,1X,13(E16.8E3,1X),I4,1X,2(E16.8E3,1X),I9,1X,F12.8,1X,2(E16.8,1x))')
1192 ityp,weight,mres,
1193 p0_el_lab,px_el_lab,py_el_lab,pz_el_lab,
1194 p0_po_lab,px_po_lab,py_po_lab,pz_po_lab,
1195 tau,time_start,time_end,
1196 flag,
1197 dens_cre,dens_abs,
1198 noe,bev,acce,accp
1199 */
1200 //# ityp weight mres p0_el_lab px_el_lab py_el_lab pz_el_lab p0_po_lab px_po_lab py_po_lab pz_po_lab tau time_start time_end flag dens_cre dens_abs noe bev acce accp
1201 // 17 9.377040e-06 1.940621e-03 3.300701e-01 -1.145014e-01 -1.798632e-01 2.519616e-01 1.025126e-01 -3.523874e-02 -5.656054e-02 7.789554e-02 3.753789e-01 3.410117e+00 3.792198e+00 2 3.075000e-01 3.360000e-01 1 3.941000e+00 1.000000e+00 1.000000e+00
1202
1203 if(out) {
1204 fprintf(out, "%5i %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e %15.8e ",
1205 ityp, weight, mres,
1206 em.E(), em.Px(), em.Py(), em.Pz(),
1207 ep.E(), ep.Px(), ep.Py(), ep.Pz(),
1208 tau, particleOut->t, particleIn->t
1209 );
1210 fprintf(out, "%4i %15.8e %15.8e %9i %15.8e %15.8e %15.8e \n",
1211 collisionIn->n_in, collisionOut->baryon_density, collisionIn->baryon_density,
1212 evtheader->evtnum, evtheader->impact_par, 1., 1.
1213 );
1214 }
1215
1216 OutputROOT(ityp, weight);
1217}
1218
1219void PHUrDilep::OutputROOT(Int_t, Double_t weight) {
1220
1221 Int_t nRequest = 1;
1222 if(outputLeptons) nRequest = 3;
1223
1224 Int_t nfree = reader->GetNFreeSlots();
1225 if(nfree < nRequest) {
1226
1227 Error("output()", "event reached max stacksize ! Will skip ....");
1228 return;
1229 }
1230 Int_t index = 0;
1231 PParticle *pMoth = reader->CreateParticle(index);
1232 PHUrAddon *pMothAdd = reader->CreateAddon(index);
1233
1234 Int_t pdg = particleIn->pdg;
1235
1236 pMoth->Reset(pdg, particleIn->Px(), particleIn->Py(), particleIn->Pz(), particleIn->E(), weight);
1237
1238 Int_t stable = 0;
1239 if(particleOut->t == particleIn ->t) stable =1;
1240
1241 pMothAdd->t_cre = particleOut ->t;
1242 pMothAdd->t_abs = particleIn ->t;
1243 pMothAdd->dens_cre = collisionOut->baryon_density;
1244 pMothAdd->dens_abs = collisionIn ->baryon_density;
1245 pMothAdd->fr_cre = particleOut ->fr;
1246 pMothAdd->fr_abs = particleIn ->fr;
1247 pMothAdd->stable = stable;
1248 pMothAdd->n_in = collisionIn->n_in;
1249 pMothAdd->n_out = collisionIn->n_out;
1250 pMothAdd->process_id = collisionIn->process_id;
1251 pMothAdd->n_collision = collisionIn->n_collision;
1252 pMothAdd->instance = particleIn->instance;
1253 pMothAdd->first = particleIn->first;
1254
1255 if(outputLeptons) {
1256 Int_t indexlep1 = 0;
1257 Int_t indexlep2 = 0;
1258 PParticle *pep = reader->CreateParticle(indexlep1);
1259 PParticle *pem = reader->CreateParticle(indexlep2);
1260 PHUrAddon *pepAdd = reader->CreateAddon(indexlep1);
1261 PHUrAddon *pemAdd = reader->CreateAddon(indexlep2);
1262
1263 pep->Reset(-11, ep.Px(), ep.Py(), ep.Pz(), ep.E(), weight); // ep
1264 pem->Reset( 11, em.Px(), em.Py(), em.Pz(), em.E(), weight); // em
1265 pem->SetParentId(pdg);
1266 pem->SetParentIndex(index);
1267 pep->SetParentId(pdg);
1268 pep->SetParentIndex(index);
1269
1270 *pepAdd = *pMothAdd;
1271 *pemAdd = *pMothAdd;
1272 }
1273}
1274
1275Double_t PHUrDilep::rndfunc() {
1276
1277 static long int init = 0;
1278
1279 if(init == 0) {
1280 srand(1);
1281 }
1282
1283 // 0 to RAND_MAX
1284 double val = ((double)rand())/RAND_MAX;
1285
1286 init++;
1287 return val;
1288}
ClassImp
ClassImp(PHUrDilep) PHUrDilep
Definition PHUrDilep.cc:16
f1
PFermiMomentumGA * f1
Definition analyze_g12C.C:104
f
TFile * f
Definition analyze_g12C.C:8
PHUrAddon::fr_abs
TVector3 fr_abs
Definition PHUrAddon.h:29
PHUrAddon::n_in
Int_t n_in
Definition PHUrAddon.h:20
PHUrAddon::process_id
Int_t process_id
Definition PHUrAddon.h:22
PHUrAddon::stable
Int_t stable
Definition PHUrAddon.h:30
PHUrAddon::instance
Int_t instance
Definition PHUrAddon.h:31
PHUrAddon::dens_abs
Double_t dens_abs
Definition PHUrAddon.h:27
PHUrAddon::t_abs
Double_t t_abs
Definition PHUrAddon.h:25
PHUrAddon::fr_cre
TVector3 fr_cre
Definition PHUrAddon.h:28
PHUrAddon::dens_cre
Double_t dens_cre
Definition PHUrAddon.h:26
PHUrAddon::first
Int_t first
Definition PHUrAddon.h:32
PHUrAddon::t_cre
Double_t t_cre
Definition PHUrAddon.h:24
PHUrAddon::n_out
Int_t n_out
Definition PHUrAddon.h:21
PHUrAddon::n_collision
Int_t n_collision
Definition PHUrAddon.h:23
PHUrDilep::omega
Int_t omega
Definition PHUrDilep.h:50
PHUrDilep::collisionIn
PHUrCollisionHeader * collisionIn
Definition PHUrDilep.h:43
PHUrDilep::em
TLorentzVector em
Definition PHUrDilep.h:45
PHUrDilep::ep
TLorentzVector ep
Definition PHUrDilep.h:45
PHUrDilep::particleOut
PHUrParticle * particleOut
Definition PHUrDilep.h:42
PHUrDilep::rho
Int_t rho
Definition PHUrDilep.h:51
PHUrDilep::evtheader
PHUrEventHeader * evtheader
Definition PHUrDilep.h:40
PHUrDilep::particleIn
PHUrParticle * particleIn
Definition PHUrDilep.h:41
PHUrDilep::Dilep
void Dilep()
Definition PHUrDilep.cc:77
PHUrDilep::out
FILE * out
Definition PHUrDilep.h:34
PHUrDilep::phi
Int_t phi
Definition PHUrDilep.h:52
PHUrDilep::omegadal
Int_t omegadal
Definition PHUrDilep.h:58
PHUrDilep::pion
Int_t pion
Definition PHUrDilep.h:53
PHUrDilep::reader
PHUrReader * reader
Definition PHUrDilep.h:29
PHUrDilep::collisionOut
PHUrCollisionHeader * collisionOut
Definition PHUrDilep.h:44
PHUrDilep::delta
Int_t delta
Definition PHUrDilep.h:56
PHUrDilep::eta
Int_t eta
Definition PHUrDilep.h:55
PHUrDilep::etaprime
Int_t etaprime
Definition PHUrDilep.h:54
PHUrDilep::outputLeptons
Bool_t outputLeptons
Definition PHUrDilep.h:65
PHUrDilep::omegadir
Int_t omegadir
Definition PHUrDilep.h:57
PParticle::SetParentIndex
void SetParentIndex(Int_t pInd)
Definition PParticle.h:234
PParticle::Reset
void Reset(const Int_t id, const TLorentzVector &v, const Double_t w=1.)
Definition PParticle.h:125
PParticle::SetParentId
void SetParentId(Int_t pId)
Definition PParticle.h:227
eta
PParticle eta("eta")
delta
PParticle * delta
Definition dp_2delta_dalitz.C:14
s
PParticle * s
Definition dp_mult.C:25
n
Int_t n
Definition etap+omeg.C:21
pdg
PPDGPlugin * pdg
Definition PDistributionManager.mod.cc:2
g
PParticle * g
Definition pp_delta_dalitz.C:18
rho
PChannelModel * rho
Definition rho_mass_shape.C:27