source: huonglan/MYSIMULATION/Saved/main_12_5.C @ 15

Last change on this file since 15 was 15, checked in by huonglan, 13 years ago

fast simulation of t't'bar pair production from my stage
File size: 21.0 KB
Line 
1#include <iostream>
2#include <math.h>
3#include "MyVector4.h"
4#include "MyVector3.h"
5#include "Particle.h"
6#include "TRandom.h"
7#include "TH1F.h"
8#include "TFile.h"
9#include "TTree.h"
10#include "LHAPDF/LHAPDF.h"
11
12using namespace std;
13
14std::vector<Particle> Decay(Particle a)
15{
16  //######## 2-BODY DECAYS ###########
17
18  short stat_a = a.STAT();
19  if (stat_a != 2)
20    {
21      std::vector<Particle> vec;
22      return vec;
23    }
24
25  MyVector4 vec4 = a.vector4();
26  MyVector3 vec3 = vec4.vector3();
27  double m_a = vec4.mass();
28  double t_a = vec4.time();
29  int pdg_a = a.PDG();
30  //int sign_pdg = pdg_a/abs(pdg_a);
31  int sign_pdg = pdg_a > 0 ? 1 : -1;
32  short pol_a = a.POL();
33
34  MyVector3 beta(vec3.X()/t_a,vec3.Y()/t_a,vec3.Z()/t_a);
35  double beta_longueur = beta.mag();
36
37  std::vector<Particle> par_data;
38
39  double m_b   = 0.;
40  int pdg_b    = 0;
41  short stat_b = 1;
42  short pol_b  = 0;
43
44  double m_c = 0.;
45  int pdg_c = 0;
46  short stat_c = 1;
47  short pol_c = 0;
48
49  if (abs(pdg_a) == 8){
50  // #### W-particle ####
51    m_b   = 80.4;
52    pdg_b    = 24*sign_pdg;
53    stat_b = 2;
54    pol_b  = 0;
55 
56  // ##### down-quark #####
57    m_c = 0.;
58    pdg_c = 1*sign_pdg;
59    stat_c = 1;
60    pol_c = 0;
61  }
62
63  else {
64  // #### up-particle ####
65    m_b   = 0.;
66    stat_b = 1;
67    pol_b  = 0;
68 
69  // ##### down-quark #####
70    m_c = 0.;
71    stat_c = 1;
72    pol_c = 0;
73
74    double R = gRandom->Uniform();
75    //  cout << "R = " << R << endl;
76
77    if (R < 0.108){
78      pdg_b = (-11)*sign_pdg;
79      pdg_c = 12*sign_pdg;
80    }
81    else if (R < 0.216){
82      pdg_b = (-13)*sign_pdg;
83      pdg_c = 14*sign_pdg;
84    }
85    else {
86      pdg_b = (2)*sign_pdg;
87      pdg_c = (-1)*sign_pdg;
88    }
89
90  }
91 
92  double E_b_CM = (m_a*m_a - m_c*m_c + m_b*m_b)/(2*m_a);
93  double p_b_CM = sqrt(E_b_CM*E_b_CM - m_b*m_b);
94 
95  double E_c_CM = (m_a*m_a - m_b*m_b + m_c*m_c)/(2*m_a);
96  double p_c_CM = sqrt(E_c_CM*E_c_CM - m_c*m_c);
97 
98  double phi = gRandom->Uniform(-M_PI,M_PI);
99  double cos_theta = gRandom->Uniform(-1,1);
100  cout << "cos_theta : " << cos_theta << endl;
101  double sin_theta = sqrt(1 - cos_theta*cos_theta);
102  cout << "sin_theta : " << sin_theta << endl;
103
104  MyVector3 vp_b_CM(p_b_CM*sin_theta*cos(phi),p_b_CM*sin_theta*sin(phi),p_b_CM*cos_theta);
105  MyVector4 vector4_data_b(vp_b_CM,E_b_CM);
106  vector4_data_b.boost(beta);
107  Particle par_b(vector4_data_b,pdg_b,stat_b,pol_b);
108  par_data.push_back(par_b);
109 
110  MyVector3 vec0(0,0,0);
111  MyVector3 vp_c_CM = vec0 - vp_b_CM;
112  MyVector4 vector4_data_c(vp_c_CM,E_c_CM);
113  vector4_data_c.boost(beta);
114  Particle par_c(vector4_data_c,pdg_c,stat_c,pol_c);
115  par_data.push_back(par_c);
116 
117  return par_data;
118}
119
120//##################################################
121
122//## NEW FUNCTIONS FOR CROSS SECTION CALCULATIONS ##
123
124//##################################################
125
126
127double alphas(double mu, double lamda = 0.24)
128{
129  double n_f = 5.;
130  double b0 = 11 - (2./3.)*n_f;
131  double b1 = 51 - (19./3.)*n_f;
132  double b2 = 2857 - (5033./9.)*n_f + (325./27.)*n_f*n_f;
133
134  double f1 = 4*M_PI/(b0*log(mu*mu/(lamda*lamda)));
135  double f2 = -(2*b1/(b0*b0))*(log(log(mu*mu/(lamda*lamda)))/log(mu*mu/(lamda*lamda)));
136  double f3 = (4*b1*b1/(b0*b0*b0*b0*pow(log(mu*mu/(lamda*lamda)),2)));
137  double f4 = pow((log(log(mu*mu/(lamda*lamda))) - 1./2.),2) + b2*b0/(8*b1*b1) - 5./4.;
138
139  double alphas = f1*(1 + f2 + f3*f4);
140
141  return alphas;
142}
143
144double dsigmaGG(double s, double t, double u, double m_Q)
145{
146  double alphaS = alphas(m_Q);
147  double s2 = s*s;
148  double t2 = t*t;
149  double u2 = u*u;
150  double m2 = m_Q*m_Q;
151
152  double coef = M_PI*alphaS*alphaS/(8*s2);
153  double f1 = 6*(t - m2)*(u - m2)/s2;
154
155  double f2 = (4./3.)*((t - m2)*(u - m2) - 2*m2*(t + m2))/((t - m2)*(t - m2)) + 3*((t - m2)*(u - m2) + m2*(u - t))/(s*(t - m2));
156
157  double f3 = (4./3.)*((u - m2)*(t - m2) - 2*m2*(u + m2))/((u - m2)*(u - m2)) + 3*((u - m2)*(t - m2) + m2*(t - u))/(s*(u - m2));
158
159  double f4 = - m2*(s - 4*m2)/(3*(t - m2)*(u - m2));
160
161  double dsigma = coef*(f1 + f2 + f3 + f4);
162
163  return dsigma;
164}
165
166double dsigmaQQ(double s_hat, double t_hat, double u_hat, double m_Q)
167{
168  double alphaS = alphas(m_Q);
169
170  double coef = 4*M_PI*alphaS*alphaS/(9*s_hat*s_hat);
171
172  double n = (t_hat - m_Q*m_Q)*(t_hat - m_Q*m_Q) + (u_hat - m_Q*m_Q)*(u_hat - m_Q*m_Q) + 2*m_Q*m_Q*s_hat;
173  double d = s_hat*s_hat;
174
175  double dsigma = coef*n/d;
176
177}
178
179//##################################
180
181//## NEW METHOD FOR CROSS SECTION ##
182
183//##################################
184
185std::vector<double> genTauYCos(double E, double m)
186{
187  double Ep = E;
188  double m_t4 = m;
189  double s = 4*Ep*Ep;
190
191  //CALCULATE TAU
192  double tau_min = 4*m_t4*m_t4/s;
193  double C1 = 0.5, C2 = 0.5;
194  double I1 = -log(tau_min);
195  double I2 = 1./tau_min - 1;
196  double tau = tau_min; //initial value
197  double r = gRandom -> Uniform();
198  if (r < C1){
199    double u = gRandom -> Uniform(log(tau_min),0);
200    tau = exp(u);
201  }
202  else {
203    double u = gRandom -> Uniform(1,1/tau_min);
204    tau = 1/u;
205  } 
206
207  //CALCULATE Y
208  double C3 = 0.4;
209  double C4 = 0.6; 
210  double y_min = log(sqrt(tau));
211  double y_max = -log(sqrt(tau));
212  double y = y_max; //initial value
213  double I3 = y_max - y_min;
214  double I4 = 2*(atan(exp(y_max)) - atan(exp(y_min)));
215  double ry = gRandom -> Uniform();
216  if (ry < C3){
217    y = gRandom -> Uniform(y_min,y_max);
218  }
219 
220  else {
221    double v = gRandom -> Uniform(2*atan(exp(y_min)),2*atan(exp(y_max)));
222    y = log(tan(v/2.));
223  }
224
225  //H_TAU AND H_Y
226  double h_tau = (C1/I1)*(1/tau) + (C2/I2)*(1/(tau*tau));
227  double h_y = C3/I3 + C4/I4/cosh(y);
228
229  //COS_THETA
230  double cos_theta = gRandom -> Uniform(-1,1);
231  double h_cos = 1/2.;
232
233  std::vector<double> vect; 
234  vect.reserve(6);
235  vect.push_back(tau);
236  vect.push_back(h_tau);
237  vect.push_back(y);
238  vect.push_back(h_y);
239  vect.push_back(cos_theta); 
240  vect.push_back(h_cos);
241
242  return vect;
243 
244}
245
246//#####################################
247//#####       MAIN PROGRAM        #####
248//#####################################
249
250#define m_t4 400
251#define Ep 3500
252#define is_pp true
253
254int main()
255{
256  cout << endl;
257  double s = 4*Ep*Ep;
258  cout << endl;
259
260  // Init LHAPDF
261  const int SUBSET  = 0;
262  const std::string NAME = "MRST2007lomod";
263  LHAPDF::initPDFSet(NAME, LHAPDF::LHGRID, SUBSET);
264  const double Q = m_t4;
265  LHAPDF::initPDF(1);
266
267  std::vector<Particle> Event;
268
269  //##### Cross section #####
270  double N = 1000000;
271 
272  double tau_min = 4*m_t4*m_t4/s;
273  double y_min = log(sqrt(tau_min));
274  double y_max = -log(sqrt(tau_min));
275
276  double sigma_gg = 0.,
277    sigma_dd = 0., 
278    sigma_uu = 0.,
279    sigma_ss = 0.,
280    sigma_cc = 0.,
281    sigma_bb = 0.;
282 
283  double N_gg = 0.;
284  double N_qq = 0.;
285
286  double t_min = 1e30;
287  double t_max = - t_min;
288
289  TH1F *tau_gg = new TH1F("tau_gg","tau_gg", 1000, 0, 1);
290  TH1F *tau_qq = new TH1F("tau_qq","tau_qq", 1000, 0, 1);
291  TH1F *tau_unweigg = new TH1F("tau_unweigg","tau_unweigg", 1000, 0, 1);
292  TH1F *tau_unweiqq = new TH1F("tau_unweiqq","tau_unweiqq", 1000, 0, 1);
293
294  float m_tau, m_y, m_ctet, m_x, m_xf[11], m_wgg, m_wqq; 
295  TFile *resF = new TFile("new_res_MC.root","recreate");
296  TTree *resT = new TTree("T3","T3");
297  resT->Branch("tau",  &m_tau,  "tau/F");
298  resT->Branch("y",    &m_y,    "y/F");
299  resT->Branch("ctet", &m_ctet, "ctet/F");
300  resT->Branch("x",    &m_x,    "x/F");
301  resT->Branch("xf",   &m_xf,   "xf[11]/F");
302  resT->Branch("wgg",  &m_wgg,  "wgg/F");
303  resT->Branch("wqq",  &m_wqq,  "wqq/F");
304
305  double weight_ggmax = 1e-10;
306  double weight_qqmax = 1e-10;
307  double x1 = 0., x2 = 0.;
308 
309  // FIND THE WEIGHT_GGMAX AND WEIGHT_QQMAX VALUES
310 
311  for (int i = 0; i < N; i++){
312
313    std::vector<double> vec = genTauYCos(Ep,m_t4);
314    double tau = vec[0];
315    double h_tau = vec[1];
316    double y = vec[2];
317    double h_y = vec[3];
318    double cos_theta = vec[4];
319    double h_cos = vec[5];
320    x1 = sqrt(tau)*exp(y);
321    x2 = sqrt(tau)*exp(-y);
322
323    double xf1_g = LHAPDF::xfx(x1,Q,0);
324    if (xf1_g < 0) xf1_g = 0;
325    double xf2_g = LHAPDF::xfx(x2,Q,0);
326    if (xf2_g < 0) xf2_g = 0;
327   
328    double xf1_d    = LHAPDF::xfx(x1,Q,1);
329    if (xf1_d < 0) xf1_d = 0;
330    double xf1_dbar = LHAPDF::xfx(x1,Q,-1);
331    if (xf1_dbar < 0) xf1_dbar = 0;
332    double xf2_d    = LHAPDF::xfx(x2,Q,1);
333    if (xf2_d < 0) xf2_d = 0;
334    double xf2_dbar = LHAPDF::xfx(x2,Q,-1);
335    if (xf2_dbar < 0) xf2_dbar = 0;
336   
337    double xf1_u    = LHAPDF::xfx(x1,Q,2);
338    if (xf1_u < 0) xf1_u = 0;
339    double xf1_ubar = LHAPDF::xfx(x1,Q,-2);
340    if (xf1_ubar < 0) xf1_ubar = 0;
341    double xf2_u    = LHAPDF::xfx(x2,Q,2);
342    if (xf2_u < 0) xf2_u = 0;
343    double xf2_ubar = LHAPDF::xfx(x2,Q,-2);
344    if (xf2_ubar < 0) xf2_ubar = 0;
345
346    double xf1_s    = LHAPDF::xfx(x1,Q,3);
347    if (xf1_s < 0) xf1_s = 0;
348    double xf2_sbar = LHAPDF::xfx(x2,Q,-3);
349    if (xf2_sbar < 0) xf2_sbar = 0;
350
351    double xf1_c    = LHAPDF::xfx(x1,Q,4);
352    if (xf1_c < 0) xf1_c = 0;
353    double xf2_cbar = LHAPDF::xfx(x2,Q,-4);
354    if (xf2_cbar < 0) xf2_cbar = 0;
355
356    double xf1_b    = LHAPDF::xfx(x1,Q,5);
357    if (xf1_b < 0) xf1_b = 0;
358    double xf2_bbar = LHAPDF::xfx(x2,Q,-5);
359    if (xf2_bbar < 0) xf2_bbar = 0;
360
361    double s_hat = x1*x2*s;
362    double beta = sqrt(1 - 4*m_t4*m_t4/s_hat);   
363    double u_hat = m_t4*m_t4 - s_hat/2.- beta*s_hat*cos_theta/2.;
364    double t_hat = m_t4*m_t4 - s_hat/2.+ beta*s_hat*cos_theta/2.;
365
366    //######## DSIGMAs & WEIGHTs ########
367    double dsigma_gg = dsigmaGG(s_hat,t_hat,u_hat,m_t4);
368    double dsigma_qq = dsigmaQQ(s_hat,t_hat,u_hat,m_t4);
369
370    double weight_gg = xf1_g * xf2_g * dsigma_gg * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
371 
372    //p-pbar case
373    double weight_dd = (xf1_d*xf2_d + xf1_dbar*xf2_dbar) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
374    double weight_uu = (xf1_u*xf2_u + xf1_ubar*xf2_ubar) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
375
376    //pp case
377    if (is_pp) 
378      {
379        weight_dd = (xf1_d * xf2_dbar + xf1_dbar * xf2_d) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
380        weight_uu = (xf1_u * xf2_ubar + xf1_ubar * xf2_u) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
381      }
382
383    double weight_ss = xf1_s * xf2_sbar * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
384    double weight_cc = xf1_c * xf2_cbar * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
385    double weight_bb = xf1_b * xf2_bbar * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
386    double weight_qq = weight_dd + weight_uu + weight_ss + weight_cc + weight_bb;
387
388    double r_gg = gRandom -> Uniform();
389    double r_qq = gRandom -> Uniform();
390
391    if (weight_gg > weight_ggmax)
392      weight_ggmax = weight_gg;
393   
394    if (weight_qq > weight_qqmax)
395      weight_qqmax = weight_qq;
396
397    sigma_gg += weight_gg;
398   
399    sigma_dd += weight_dd;
400    sigma_uu += weight_uu;
401    sigma_ss += weight_ss;
402    sigma_cc += weight_cc;
403    sigma_bb += weight_bb;
404
405  }
406  cout << "*********************************" << endl; 
407  cout << "weight_ggmax : " << weight_ggmax << endl;
408  cout << "weight_qqmax : " << weight_qqmax << endl; 
409  cout << "*********************************" << endl; 
410 
411  //######################################
412  //######## GENERATE REAL EVENTS ########
413  //######################################
414
415  for (int i = 0; i < N; i++){
416   
417    if (i%10000 == 0)
418      cout << "generating event " << i << endl;
419   
420    //##################################
421    std::vector<double> vec = genTauYCos(3500,400);
422    double tau = vec[0];
423    double h_tau = vec[1];
424    double y = vec[2];
425    double h_y = vec[3];
426    double cos_theta = vec[4];
427    double h_cos = vec[5];
428   
429    //########## X1 & X2 ##########
430    x1 = sqrt(tau)*exp(y);
431    x2 = sqrt(tau)*exp(-y);
432
433    //######## USING PDFS #########
434    double xf1_g = LHAPDF::xfx(x1,Q,0);
435    if (xf1_g < 0) xf1_g = 0;
436    double xf2_g = LHAPDF::xfx(x2,Q,0);
437    if (xf2_g < 0) xf2_g = 0;
438   
439    double xf1_d    = LHAPDF::xfx(x1,Q,1);
440    if (xf1_d < 0) xf1_d = 0;
441    double xf1_dbar = LHAPDF::xfx(x1,Q,-1);
442    if (xf1_dbar < 0) xf1_dbar = 0;
443    double xf2_d    = LHAPDF::xfx(x2,Q,1);
444    if (xf2_d < 0) xf2_d = 0;
445    double xf2_dbar = LHAPDF::xfx(x2,Q,-1);
446    if (xf2_dbar < 0) xf2_dbar = 0;
447   
448    double xf1_u    = LHAPDF::xfx(x1,Q,2);
449    if (xf1_u < 0) xf1_u = 0;
450    double xf1_ubar = LHAPDF::xfx(x1,Q,-2);
451    if (xf1_ubar < 0) xf1_ubar = 0;
452    double xf2_u    = LHAPDF::xfx(x2,Q,2);
453    if (xf2_u < 0) xf2_u = 0;
454    double xf2_ubar = LHAPDF::xfx(x2,Q,-2);
455    if (xf2_ubar < 0) xf2_ubar = 0;
456
457    double xf1_s    = LHAPDF::xfx(x1,Q,3);
458    if (xf1_s < 0) xf1_s = 0;
459    double xf2_sbar = LHAPDF::xfx(x2,Q,-3);
460    if (xf2_sbar < 0) xf2_sbar = 0;
461
462    double xf1_c    = LHAPDF::xfx(x1,Q,4);
463    if (xf1_c < 0) xf1_c = 0;
464    double xf2_cbar = LHAPDF::xfx(x2,Q,-4);
465    if (xf2_cbar < 0) xf2_cbar = 0;
466
467    double xf1_b    = LHAPDF::xfx(x1,Q,5);
468    if (xf1_b < 0) xf1_b = 0;
469    double xf2_bbar = LHAPDF::xfx(x2,Q,-5);
470    if (xf2_bbar < 0) xf2_bbar = 0;
471
472    //######## FILL THE TREE ########
473    m_tau   = tau;
474    m_y     = y;
475    m_ctet = cos_theta;
476    m_x     = x1;
477    m_xf[0] = xf1_g;
478    m_xf[1] = xf2_g;
479    m_xf[2] = xf1_d;
480    m_xf[3] = xf2_dbar;
481    m_xf[4] = xf1_u;
482    m_xf[5] = xf2_ubar;
483    m_xf[6] = xf1_s;
484    m_xf[7] = xf2_sbar;
485    m_xf[8] = xf1_c;
486    m_xf[9] = xf2_cbar;
487    m_xf[10] = xf1_b;
488    m_xf[11] = xf2_bbar;
489
490    MyVector4 vg1(MyVector3(0,0,x1*Ep),x1*Ep);
491    Particle g1(vg1,21,3,0);
492    MyVector4 vg2(MyVector3(0,0,- x2*Ep),x2*Ep);
493    Particle g2(vg2,21,3,0);
494   
495    Event.push_back(g1);
496    Event.push_back(g2);
497 
498    double s_hat = x1*x2*s;
499    double beta = sqrt(1 - 4*m_t4*m_t4/s_hat);   
500    double u_hat = m_t4*m_t4 - s_hat/2.- beta*s_hat*cos_theta/2.;
501    double t_hat = m_t4*m_t4 - s_hat/2.+ beta*s_hat*cos_theta/2.;
502
503    if (t_hat > t_max)
504      t_max = t_hat;
505
506    if (t_hat < t_min)
507      t_min = t_hat;
508
509
510    //######## DSIGMAs & WEIGHTs ########
511    double dsigma_gg = dsigmaGG(s_hat,t_hat,u_hat,m_t4);
512    double dsigma_qq = dsigmaQQ(s_hat,t_hat,u_hat,m_t4);
513
514    double weight_gg = xf1_g * xf2_g * dsigma_gg * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
515 
516    //p-pbar case
517    double weight_dd = (xf1_d*xf2_d + xf1_dbar*xf2_dbar) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
518    double weight_uu = (xf1_u*xf2_u + xf1_ubar*xf2_ubar) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
519
520    //pp case
521    if (is_pp) 
522      {
523        weight_dd = (xf1_d * xf2_dbar + xf1_dbar * xf2_d) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
524        weight_uu = (xf1_u * xf2_ubar + xf1_ubar * xf2_u) * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
525      }
526
527    double weight_ss = xf1_s * xf2_sbar * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
528    double weight_cc = xf1_c * xf2_cbar * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
529    double weight_bb = xf1_b * xf2_bbar * dsigma_qq * s_hat * s_hat * beta /(2*tau*tau * h_tau * h_y * h_cos * M_PI);
530    double weight_qq = weight_dd + weight_uu + weight_ss + weight_cc + weight_bb;
531    m_wgg = weight_gg;
532    m_wqq = weight_qq;   
533
534    //######## FIL HISTOS ###########
535
536    tau_gg -> Fill(tau,weight_gg);
537    tau_qq -> Fill(tau,weight_qq);
538
539    //###############################
540
541    double r_gg = gRandom -> Uniform();
542    double r_qq = gRandom -> Uniform();
543   
544    if (weight_gg > r_gg*weight_ggmax){
545      N_gg ++;
546      tau_unweigg -> Fill(tau);
547    }
548   
549    if (weight_qq > r_qq*weight_qqmax){
550      N_qq ++;
551      tau_unweiqq -> Fill(tau);
552    }
553   
554    resT->Fill();
555   
556  }//END for (int i = 0; i < N; i++)
557 
558  resT->Write();
559  tau_gg -> Write();
560  tau_qq -> Write();
561  tau_unweigg -> Write();
562  tau_unweiqq -> Write();
563  resF->Close();
564
565  cout << "t_min : " << t_min << "   t_max : " << t_max << endl;
566 
567  double dtau = 1 - tau_min;
568  double dy = y_max - y_min;
569  double dcost = 2.;
570
571  // Here we dont have to multiply sigma
572  // *** with the phase-space volume ***
573  sigma_gg = M_PI*sigma_gg/N/s;
574  sigma_dd = M_PI*sigma_dd/N/s;
575  sigma_uu = M_PI*sigma_uu/N/s;
576  sigma_ss = M_PI*2*sigma_ss/N/s;
577  sigma_cc = M_PI*2*sigma_cc/N/s;
578  sigma_bb = M_PI*2*sigma_bb/N/s;
579 
580  cout << "sigma_gg = " << sigma_gg << endl << endl;
581  cout << "sigma_dd = " << sigma_dd << endl;
582  cout << "sigma_uu = " << sigma_uu << endl;
583  cout << "sigma_ss = " << sigma_ss << endl;
584  cout << "sigma_cc = " << sigma_cc << endl;
585  cout << "sigma_bb = " << sigma_bb << endl << endl;
586
587  double sigma_qq = sigma_dd + sigma_uu + sigma_ss + sigma_cc + sigma_bb;
588 
589  //### CONVERT THE UNIT OF CROSS SECTION TO PB
590
591  double sigma_convert1 = sigma_gg*0.389*1e9;
592  double sigma_convert2 = sigma_dd*0.389*1e9;
593  double sigma_convert3 = sigma_uu*0.389*1e9;
594  double sigma_convert4 = sigma_ss*0.389*1e9;
595  double sigma_convert5 = sigma_cc*0.389*1e9;
596  double sigma_convert6 = sigma_bb*0.389*1e9;
597  double sigma_convert7 = sigma_qq*0.389*1e9;
598
599  cout << "After convert the unit of cross section: " << endl;
600
601  cout << "gg->t't'bar      : " << sigma_convert1 << endl << endl;
602  cout << "d-dbar->t't'bar  : " << sigma_convert2 << endl;
603  cout << "u-ubar->t't'bar  : " << sigma_convert3 << endl;
604  cout << "s-sbar->t't'bar  : " << sigma_convert4 << endl;
605  cout << "c-cbar->t't'bar  : " << sigma_convert5 << endl;
606  cout << "b-bbar->t't'bar  : " << sigma_convert6 << endl << endl;
607  cout << "Total of the processes q-qbar->t't'bar  : " << sigma_convert7 << endl << endl;
608
609  cout << "TOTAL : " << sigma_convert1 + sigma_convert7 << endl;
610
611  cout << "number of events for gg : " << N_gg << endl;
612  cout << "number of events for qq : " << N_qq << endl;
613
614  cout << "***********************************" << endl << endl;
615 
616  //############ PRODUCTION ############
617  double r = gRandom -> Uniform();
618  double r_gg = sigma_convert1/(sigma_convert1 + sigma_convert7);
619
620  cout << "r = " << r << ";    r_gg = " << r_gg << endl;
621
622  if (r < r_gg){
623
624    std::vector<Particle> Event;
625
626    MyVector4 vg1(MyVector3(0,0,x1*Ep),x1*Ep);
627    Particle g1(vg1,21,3,0);
628    MyVector4 vg2(MyVector3(0,0,- x2*Ep),x2*Ep);
629    Particle g2(vg2,21,3,0);
630
631    Event.push_back(g1);
632    Event.push_back(g2);
633
634    MyVector4 g_total = vg1 + vg2;
635    //    cout << "g1 + g2 : " << g_diff << endl;
636    double z = g_total.vector3().Z();
637    double t = g_total.time();
638
639    MyVector3 beta_g(0,0,z/t);
640
641    double s_hat = x1*x2*s;
642    double E_CM = sqrt(s_hat);
643
644    //for top production:
645    double E_t_CM = E_CM/2;
646    double p_t_CM = sqrt(E_t_CM*E_t_CM - m_t4*m_t4);
647
648    double phi_s = gRandom->Uniform(-M_PI,M_PI);
649    double cos_theta_s = gRandom->Uniform(-1,1);
650    double sin_theta_s = sqrt(1 - cos_theta_s*cos_theta_s);
651
652    MyVector3 vt31(p_t_CM*sin_theta_s*cos(phi_s),p_t_CM*sin_theta_s*sin(phi_s),p_t_CM*cos_theta_s);
653    MyVector3 vec0(0,0,0);
654    MyVector3 vt32 = vec0 - vt31;
655
656    MyVector4 vt41(vt31,E_t_CM);
657    vt41.boost(beta_g);
658    MyVector4 vt42(vt32,E_t_CM);
659    vt42.boost(beta_g);
660
661    Particle t4(vt41,8,2,0);
662    Event.push_back(t4);
663    Particle at4(vt42,-8,2,0);
664    Event.push_back(at4);
665
666    cout << "****************************************************************** "  << endl;
667
668    for (int i = 0; i < Event.size(); i++){
669      Particle p = Event.at(i);
670      cout << "p " << i+1 << " : "  << p << endl;
671      cout << "****************************************************************** "  << endl;
672      if (p.STAT() == 2){
673        std::vector<Particle> p_dec = Decay(p);
674        for (int j = 0; j<p_dec.size(); j++)
675          Event.push_back(p_dec.at(j));
676      }
677    }
678  }//END of  if (r < r_gg)
679 
680  else {
681    std::vector<Particle> Event;
682
683    MyVector4 vq(MyVector3(0,0,x1*Ep),x1*Ep);
684    Particle q(vq,2,1,0);
685    MyVector4 vqbar(MyVector3(0,0,- x2*Ep),x2*Ep);
686    Particle qbar(vqbar,-2,1,0);
687
688    Event.push_back(q);
689    Event.push_back(qbar);
690
691    MyVector4 q_total = vq + vqbar;
692    //    cout << "g1 + g2 : " << g_diff << endl;
693    double z = q_total.vector3().Z();
694    double t = q_total.time();
695
696    MyVector3 beta_g(0,0,z/t);
697
698    double s_hat = x1*x2*s;
699    double E_CM = sqrt(s_hat);
700
701    //for top production:
702    double E_t_CM = E_CM/2;
703    double p_t_CM = sqrt(E_t_CM*E_t_CM - m_t4*m_t4);
704   
705    double phi_s = gRandom->Uniform(-M_PI,M_PI);
706    double cos_theta_s = gRandom->Uniform(-1,1);
707    double sin_theta_s = sqrt(1 - cos_theta_s*cos_theta_s);
708   
709    MyVector3 vt31(p_t_CM*sin_theta_s*cos(phi_s),p_t_CM*sin_theta_s*sin(phi_s),p_t_CM*cos_theta_s);
710    MyVector3 vec0(0,0,0);
711    MyVector3 vt32 = vec0 - vt31;
712
713    MyVector4 vt41(vt31,E_t_CM);
714    vt41.boost(beta_g);
715    MyVector4 vt42(vt32,E_t_CM);
716    vt42.boost(beta_g);
717
718    Particle t4(vt41,8,2,0);
719    Event.push_back(t4);
720    Particle at4(vt42,-8,2,0);
721    Event.push_back(at4);
722
723    cout << "****************************************************************** "  << endl;
724
725    for (int i = 0; i < Event.size(); i++){
726      Particle p = Event.at(i);
727      cout << "p " << i+1 << " : "  << p << endl;
728      cout << "****************************************************************** "  << endl;
729      if (p.STAT() == 2){
730        std::vector<Particle> p_dec = Decay(p);
731        for (int j = 0; j<p_dec.size(); j++)
732          Event.push_back(p_dec.at(j));
733      }
734    }
735   
736  }//END of else after  if (r < r_gg)
737       
738  return 0;
739}
740
741
742#include "LHAPDF/FortranWrappers.h"
743#ifdef FC_DUMMY_MAIN
744int FC_DUMMY_MAIN() {return l;}
745#endif
Note: See TracBrowser for help on using the repository browser.