source: huonglan/MYSIMULATION/main.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.9 KB
Line 
1#include <iostream>
2#include <math.h>
3#include "MyVector4.h"
4#include "MyVector3.h"
5#include "Particle.h"
6//#include "Decay.h"//For top_prime
7#include "Decaybottom.h"//For bottom prime with 2 same sign leptons
8#include "Functions.C"
9#include "TRandom.h"
10#include "TH1F.h"
11#include "TF1.h"
12#include "TFile.h"
13#include "TTree.h"
14#include "LHAPDF/LHAPDF.h"
15
16using namespace std;
17
18MyVector4 smear(Particle p){
19
20  int    pdg = p.PDG();
21  double E   = p.vector4().E();
22  double pT  = p.vector4().Pt();
23  double eta = p.vector4().eta();
24
25  //****** case of electrons ******
26  if (abs(pdg) == 11){
27    double sigma_Ee;
28    if (fabs(eta) < 1.4){
29      sigma_Ee = E*0.095/sqrt(E);
30    }
31    if (fabs(eta) > 1.4){
32      sigma_Ee = E*0.195/sqrt(E);
33    }
34   
35    double E_smeared = gRandom->Gaus(E,sigma_Ee);
36    double k = E_smeared/E;
37    MyVector3 v3 = p.vector4().vector3();
38    MyVector4 v4(k*v3,E_smeared);
39
40    return v4;
41  }
42
43  //****** case of jets ******
44  if (abs(pdg) <= 5){
45    double sigma_Eq;
46    if (fabs(eta) < 0.5){
47      sigma_Eq = E*0.65/sqrt(E);
48    }
49    if (fabs(eta) > 0.5 && fabs(eta) < 1.5){
50      sigma_Eq = sqrt(E)*(0.35*fabs(eta) + 0.475);
51    }
52    if (fabs(eta) > 1.5 && fabs(eta) < 5){
53      sigma_Eq = E*1.0/sqrt(E);
54    }
55       
56    double E_smeared = gRandom->Gaus(E,sigma_Eq);
57    double k = E_smeared/E;
58    MyVector3 v3 = p.vector4().vector3();
59    MyVector4 v4(k*v3,E_smeared);
60
61    return v4;
62  }
63
64  //****** case of muons ******
65  if (abs(pdg) == 13){
66    double delta_pTm;
67    if (fabs(eta) < 1.0){
68      delta_pTm = 0.02*pT;
69    }
70    if (fabs(eta) > 1.0 && fabs(eta) < 1.8){
71      delta_pTm = pT*(0.225*fabs(eta) - 0.025);
72    }
73    if (fabs(eta) > 1.8 && fabs(eta) < 2.0){
74      delta_pTm = pT*(-0.4*fabs(eta) + 1.1);
75    }
76    if (fabs(eta) > 2.0){
77      delta_pTm = 0.03*pT;
78    }
79   
80    double pTm_smeared = gRandom->Gaus(pT,delta_pTm);
81    double k = pTm_smeared/pT;
82    MyVector3 v3 = p.vector4().vector3();
83    MyVector4 v4(k*v3,k*E);
84
85    return v4;
86  }
87
88  else {
89    cout << "Unknown particle type!" << "  pdg : " << pdg << endl;
90    return p.vector4();
91  }
92}
93
94
95//############################################
96//************** MAIN PROGRAM ****************
97//############################################
98
99#define m_t4 400
100#define Ep 3500
101#define is_pp true
102
103int main()
104{
105  cout << endl;
106  double s = 4*Ep*Ep;
107  cout << endl;
108
109  // Init LHAPDF
110  const int SUBSET  = 0;
111  const std::string NAME = "MRST2007lomod";
112  LHAPDF::initPDFSet(NAME, LHAPDF::LHGRID, SUBSET);
113  const double Q = m_t4;
114  LHAPDF::initPDF(1);
115
116  std::vector<Particle> Event;
117
118  //##### Cross section #####
119  double N = 1000000;
120 
121  double tau_min = 4*m_t4*m_t4/s;
122  double y_min = log(sqrt(tau_min));
123  double y_max = -log(sqrt(tau_min));
124
125  double sigma_gg = 0.,
126    sigma_dd = 0., 
127    sigma_uu = 0.,
128    sigma_ss = 0.,
129    sigma_cc = 0.,
130    sigma_bb = 0.;
131 
132  double t_min = 1e30;
133  double t_max = - t_min;
134
135  TFile *resF = new TFile("MC_bottom_prime.root","recreate");
136  TTree *resT = new TTree("T3","T3");
137  TTree *resT_ana = new TTree("T_ana","T_ana");
138
139  float m_tau, m_y, m_ctet_top, m_x, m_xf[11], m_wgg, m_wqq;
140
141  resT->Branch("tau",  &m_tau,  "tau/F");
142  resT->Branch("y",    &m_y,    "y/F");
143  resT->Branch("ctet_top", &m_ctet_top, "ctet_top/F");
144  resT->Branch("x",    &m_x,    "x/F");
145  resT->Branch("xf",   &m_xf,   "xf[11]/F");
146  resT->Branch("wgg",  &m_wgg,  "wgg/F");
147  resT->Branch("wqq",  &m_wqq,  "wqq/F");
148
149  // An EventNumber
150  int m_evt;
151  resT_ana->Branch("evt",  &m_evt,  "evt/I");
152
153  // Truth information
154  static const int m_nparmax = 30;
155  int m_npar, m_pdgId[m_nparmax], m_stat[m_nparmax], m_barc[m_nparmax];
156  float m_etaP[m_nparmax], m_phiP[m_nparmax], m_pTP[m_nparmax], m_EP[m_nparmax];
157  resT_ana->Branch("npar",  &m_npar,  "npar/I");
158  resT_ana->Branch("pdgId", &m_pdgId, "pdgId[npar]/I");
159  resT_ana->Branch("stat",  &m_stat,  "stat[npar]/I");
160  resT_ana->Branch("barc",  &m_barc,  "barc[npar]/I");
161  resT_ana->Branch("etaP",  &m_etaP,  "etaP[npar]/F");
162  resT_ana->Branch("phiP",  &m_phiP,  "phiP[npar]/F");
163  resT_ana->Branch("pTP",   &m_pTP,   "pTP[npar]/F");
164  resT_ana->Branch("EP",    &m_EP,    "EP[npar]/F");
165
166  // Electrons
167  static const int m_nelemax = 10;
168  int m_nele, m_pdgIdE[m_nelemax];
169  float m_etaE[m_nelemax], m_phiE[m_nelemax], m_pTE[m_nelemax], m_EE[m_nelemax];
170
171  resT_ana->Branch("nele",  &m_nele,  "nele/I");
172  resT_ana->Branch("pdgIdE", &m_pdgIdE, "pdgIdE[nele]/I");
173  resT_ana->Branch("etaE",  &m_etaE,  "etaE[nele]/F");
174  resT_ana->Branch("phiE",  &m_phiE,  "phiE[nele]/F");
175  resT_ana->Branch("pTE",   &m_pTE,   "pTE[nele]/F");
176  resT_ana->Branch("EE",    &m_EE,    "EE[nele]/F");
177
178  // Jets
179  static const int m_nqmax = 20;
180  int m_nq, m_pdgIdQ[m_nqmax];
181  float m_etaQ[m_nqmax], m_phiQ[m_nqmax], m_pTQ[m_nqmax], m_EQ[m_nqmax];
182
183  resT_ana->Branch("nq",  &m_nq,  "nq/I");
184  resT_ana->Branch("pdgIdQ", &m_pdgIdQ, "pdgIdQ[nq]/I");
185  resT_ana->Branch("etaQ",  &m_etaQ,  "etaQ[nq]/F");
186  resT_ana->Branch("phiQ",  &m_phiQ,  "phiQ[nq]/F");
187  resT_ana->Branch("pTQ",   &m_pTQ,   "pTQ[nq]/F");
188  resT_ana->Branch("EQ",    &m_EQ,    "EQ[nq]/F");
189
190  // Muons
191  static const int m_nmumax = 10;
192  int m_nmu, m_pdgIdM[m_nmumax];
193  float m_etaM[m_nmumax], m_phiM[m_nmumax], m_pTM[m_nmumax], m_EM[m_nmumax];
194
195  resT_ana->Branch("nmu",  &m_nmu,  "nmu/I");
196  resT_ana->Branch("pdgIdM", &m_pdgIdM, "pdgIdM[nmu]/I");
197  resT_ana->Branch("etaM",  &m_etaM,  "etaM[nmu]/F");
198  resT_ana->Branch("phiM",  &m_phiM,  "phiM[nmu]/F");
199  resT_ana->Branch("pTM",   &m_pTM,   "pTM[nmu]/F");
200  resT_ana->Branch("EM",    &m_EM,    "EM[nmu]/F");
201
202  //Missing Energy
203  float m_ETmisX, m_ETmisY, m_ETmisZ, m_SumET;
204  resT_ana->Branch("ETmisX",  &m_ETmisX,  "ETmisX/F");
205  resT_ana->Branch("ETmisY",  &m_ETmisY,  "ETmisY/F");
206  resT_ana->Branch("ETmisZ",  &m_ETmisZ,  "ETmisZ/F");
207  resT_ana->Branch("SumET",   &m_SumET,   "SumET/F");
208
209  //Measured Missing Energy
210  float m_ETmisXmea, m_ETmisYmea;
211  resT_ana->Branch("ETmisXmea",  &m_ETmisXmea,  "ETmisXmea/F");
212  resT_ana->Branch("ETmisYmea",  &m_ETmisYmea,  "ETmisYmea/F");
213 
214
215  //############################################## 
216  // FIND THE WEIGHT_GGMAX AND WEIGHT_QQMAX VALUES
217  //############################################## 
218
219  double weight_ggmax = 1e-10;
220  double weight_qqmax = 1e-10;
221  double x1 = 0., x2 = 0.;
222 
223  for (int i = 0; i < N; i++){
224
225    std::vector<double> vec = genTauYCos(Ep,m_t4);
226    double tau = vec[0];
227    double h_tau = vec[1];
228    double y = vec[2];
229    double h_y = vec[3];
230    double cos_theta = vec[4];
231    double h_cos = vec[5];
232    x1 = sqrt(tau)*exp(y);
233    x2 = sqrt(tau)*exp(-y);
234
235    double xf1_g = LHAPDF::xfx(x1,Q,0);
236    if (xf1_g < 0) xf1_g = 0;
237    double xf2_g = LHAPDF::xfx(x2,Q,0);
238    if (xf2_g < 0) xf2_g = 0;
239   
240    double xf1_d    = LHAPDF::xfx(x1,Q,1);
241    if (xf1_d < 0) xf1_d = 0;
242    double xf1_dbar = LHAPDF::xfx(x1,Q,-1);
243    if (xf1_dbar < 0) xf1_dbar = 0;
244    double xf2_d    = LHAPDF::xfx(x2,Q,1);
245    if (xf2_d < 0) xf2_d = 0;
246    double xf2_dbar = LHAPDF::xfx(x2,Q,-1);
247    if (xf2_dbar < 0) xf2_dbar = 0;
248   
249    double xf1_u    = LHAPDF::xfx(x1,Q,2);
250    if (xf1_u < 0) xf1_u = 0;
251    double xf1_ubar = LHAPDF::xfx(x1,Q,-2);
252    if (xf1_ubar < 0) xf1_ubar = 0;
253    double xf2_u    = LHAPDF::xfx(x2,Q,2);
254    if (xf2_u < 0) xf2_u = 0;
255    double xf2_ubar = LHAPDF::xfx(x2,Q,-2);
256    if (xf2_ubar < 0) xf2_ubar = 0;
257
258    double xf1_s    = LHAPDF::xfx(x1,Q,3);
259    if (xf1_s < 0) xf1_s = 0;
260    double xf2_sbar = LHAPDF::xfx(x2,Q,-3);
261    if (xf2_sbar < 0) xf2_sbar = 0;
262
263    double xf1_c    = LHAPDF::xfx(x1,Q,4);
264    if (xf1_c < 0) xf1_c = 0;
265    double xf2_cbar = LHAPDF::xfx(x2,Q,-4);
266    if (xf2_cbar < 0) xf2_cbar = 0;
267
268    double xf1_b    = LHAPDF::xfx(x1,Q,5);
269    if (xf1_b < 0) xf1_b = 0;
270    double xf2_bbar = LHAPDF::xfx(x2,Q,-5);
271    if (xf2_bbar < 0) xf2_bbar = 0;
272
273    double s_hat = x1*x2*s;
274    double beta = sqrt(1 - 4*m_t4*m_t4/s_hat);   
275    double u_hat = m_t4*m_t4 - s_hat/2.- beta*s_hat*cos_theta/2.;
276    double t_hat = m_t4*m_t4 - s_hat/2.+ beta*s_hat*cos_theta/2.;
277
278    //######## DSIGMAs & WEIGHTs ########
279    double dsigma_gg = dsigmaGG(s_hat,t_hat,u_hat,m_t4);
280    double dsigma_qq = dsigmaQQ(s_hat,t_hat,u_hat,m_t4);
281
282    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);
283 
284    //p-pbar case
285    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);
286    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);
287
288    //pp case
289    if (is_pp) 
290      {
291        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);
292        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);
293      }
294
295    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);
296    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);
297    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);
298    double weight_qq = weight_dd + weight_uu + weight_ss + weight_cc + weight_bb;
299
300    double r_gg = gRandom -> Uniform();
301    double r_qq = gRandom -> Uniform();
302
303    if (weight_gg > weight_ggmax)
304      weight_ggmax = weight_gg;
305   
306    if (weight_qq > weight_qqmax)
307      weight_qqmax = weight_qq;
308
309    sigma_gg += weight_gg;
310   
311    sigma_dd += weight_dd;
312    sigma_uu += weight_uu;
313    sigma_ss += weight_ss;
314    sigma_cc += weight_cc;
315    sigma_bb += weight_bb;
316
317  }
318 
319  double dtau = 1 - tau_min;
320  double dy = y_max - y_min;
321  double dcost = 2.;
322
323  // Here we dont have to multiply sigma
324  // *** with the phase-space volume ***
325  sigma_gg = M_PI*sigma_gg/N/s;
326  sigma_dd = M_PI*sigma_dd/N/s;
327  sigma_uu = M_PI*sigma_uu/N/s;
328  sigma_ss = M_PI*2*sigma_ss/N/s;
329  sigma_cc = M_PI*2*sigma_cc/N/s;
330  sigma_bb = M_PI*2*sigma_bb/N/s;
331
332  double sigma_qq = sigma_dd + sigma_uu + sigma_ss + sigma_cc + sigma_bb;
333
334  //### CONVERT THE UNIT OF CROSS SECTION TO PB
335
336  double sigma_convert1 = sigma_gg*0.389*1e9;
337  double sigma_convert2 = sigma_dd*0.389*1e9;
338  double sigma_convert3 = sigma_uu*0.389*1e9;
339  double sigma_convert4 = sigma_ss*0.389*1e9;
340  double sigma_convert5 = sigma_cc*0.389*1e9;
341  double sigma_convert6 = sigma_bb*0.389*1e9;
342  double sigma_convert7 = sigma_qq*0.389*1e9;
343
344
345  //######################################################################
346
347  //####################### GENERATE REAL EVENTS #########################
348
349  //######################################################################
350
351  double N1 = 1000000;
352  double N_gg = 0.;
353  double N_qq = 0.;
354
355  for (int i = 0; i < N1; i++){
356    if (i%10000 == 0)
357      cout << "generating event " << i << endl;
358
359    //##################################
360    std::vector<double> vec = genTauYCos(3500,400);
361    double tau = vec[0];
362    double h_tau = vec[1];
363    double y = vec[2];
364    double h_y = vec[3];
365    double cos_theta = vec[4];
366    double sin_theta = sqrt(1-cos_theta*cos_theta);
367    double h_cos = vec[5];
368
369    //########## X1 & X2 ##########
370    x1 = sqrt(tau)*exp(y);
371    x2 = sqrt(tau)*exp(-y);
372
373    double s_hat = x1*x2*s;
374
375    double beta = sqrt(1 - 4*m_t4*m_t4/s_hat);   
376    double u_hat = m_t4*m_t4 - s_hat/2.- beta*s_hat*cos_theta/2.;
377    double t_hat = m_t4*m_t4 - s_hat/2.+ beta*s_hat*cos_theta/2.;
378   
379    //### fill tree T3 ###
380    m_tau   = tau;
381    m_y     = y;
382    m_ctet_top = cos_theta;
383    m_x     = x1;
384
385    std::vector<Particle> Event;
386    MyVector4 vp1,vp2;
387       
388    double r_gg = gRandom -> Uniform();
389    double k_gg = sigma_convert1/(sigma_convert1 + sigma_convert7);
390           
391    if (r_gg < k_gg){
392      double xf1_g = LHAPDF::xfx(x1,Q,0);
393      if (xf1_g < 0) xf1_g = 0;
394      double xf2_g = LHAPDF::xfx(x2,Q,0);
395      if (xf2_g < 0) xf2_g = 0;
396
397      double dsigma_gg = dsigmaGG(s_hat,t_hat,u_hat,m_t4);     
398      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);
399
400      //### fill tree T3 ###
401      m_xf[0] = xf1_g;
402      m_xf[1] = xf2_g;
403      m_wgg = weight_gg;
404 
405      double r = gRandom -> Uniform();
406   
407      if (weight_gg/k_gg < r*(weight_ggmax/k_gg + weight_qqmax/(1.- k_gg))) continue;
408      else {
409        N_gg ++;
410       
411        MyVector4 vp11(MyVector3(0,0,x1*Ep),x1*Ep);
412        vp1 = vp11;
413        Particle p1(vp1,21,3,0);
414        p1.barcode(1);
415        p1.motherbarcode(0);
416        MyVector4 vp22(MyVector3(0,0,- x2*Ep),x2*Ep);
417        vp2 = vp22;
418        Particle p2(vp2,21,3,0);
419        p2.barcode(2);
420        p2.motherbarcode(0);
421       
422        Event.push_back(p1);
423        Event.push_back(p2);
424      }
425    }//END of <if (r_gg < k_gg)>
426
427    else {
428      //cout << "qqbar -> QQbar" << endl;
429      double xf1_d    = LHAPDF::xfx(x1,Q,1);
430      if (xf1_d < 0) xf1_d = 0;
431      double xf1_dbar = LHAPDF::xfx(x1,Q,-1);
432      if (xf1_dbar < 0) xf1_dbar = 0;
433      double xf2_d    = LHAPDF::xfx(x2,Q,1);
434      if (xf2_d < 0) xf2_d = 0;
435      double xf2_dbar = LHAPDF::xfx(x2,Q,-1);
436      if (xf2_dbar < 0) xf2_dbar = 0;
437     
438      double xf1_u    = LHAPDF::xfx(x1,Q,2);
439      if (xf1_u < 0) xf1_u = 0;
440      double xf1_ubar = LHAPDF::xfx(x1,Q,-2);
441      if (xf1_ubar < 0) xf1_ubar = 0;
442      double xf2_u    = LHAPDF::xfx(x2,Q,2);
443      if (xf2_u < 0) xf2_u = 0;
444      double xf2_ubar = LHAPDF::xfx(x2,Q,-2);
445      if (xf2_ubar < 0) xf2_ubar = 0;
446     
447      double xf1_s    = LHAPDF::xfx(x1,Q,3);
448      if (xf1_s < 0) xf1_s = 0;
449      double xf2_sbar = LHAPDF::xfx(x2,Q,-3);
450      if (xf2_sbar < 0) xf2_sbar = 0;
451     
452      double xf1_c    = LHAPDF::xfx(x1,Q,4);
453      if (xf1_c < 0) xf1_c = 0;
454      double xf2_cbar = LHAPDF::xfx(x2,Q,-4);
455      if (xf2_cbar < 0) xf2_cbar = 0;
456     
457      double xf1_b    = LHAPDF::xfx(x1,Q,5);
458      if (xf1_b < 0) xf1_b = 0;
459      double xf2_bbar = LHAPDF::xfx(x2,Q,-5);
460      if (xf2_bbar < 0) xf2_bbar = 0;
461
462      double dsigma_qq = dsigmaQQ(s_hat,t_hat,u_hat,m_t4);
463
464      //p-pbar case
465      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);
466      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);
467
468      //pp case
469      if (is_pp) {
470        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);
471        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);
472      }
473
474      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);
475      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);
476      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);
477      double weight_qq = weight_dd + weight_uu + weight_ss + weight_cc + weight_bb;
478
479      //### Fill tree T3 ###
480      m_xf[2] = xf1_d;
481      m_xf[3] = xf2_dbar;
482      m_xf[4] = xf1_u;
483      m_xf[5] = xf2_ubar;
484      m_xf[6] = xf1_s;
485      m_xf[7] = xf2_sbar;
486      m_xf[8] = xf1_c;
487      m_xf[9] = xf2_cbar;
488      m_xf[10] = xf1_b;
489      m_xf[11] = xf2_bbar;
490      m_wqq = weight_qq;   
491
492      double r = gRandom -> Uniform();
493
494      if (weight_qq/(1.- k_gg) < r*(weight_ggmax/k_gg + weight_qqmax/(1.- k_gg))) continue;
495      else {
496        N_qq ++;
497       
498        MyVector4 vp11(MyVector3(0,0,x1*Ep),x1*Ep);
499        vp1 = vp11;
500        Particle p1(vp1,2,1,0);
501        p1.barcode(1);
502        p1.motherbarcode(0);
503        MyVector4 vp22(MyVector3(0,0,- x2*Ep),x2*Ep);
504        vp2 = vp22;
505        Particle p2(vp2,-2,1,0);
506        p2.barcode(2);
507        p2.motherbarcode(0);
508       
509        Event.push_back(p1);
510        Event.push_back(p2);
511      }
512
513    }//END of else of <if (r_gg < k_gg)>
514
515    MyVector4 v_total = vp1 + vp2;
516    double z = v_total.vector3().Z();
517    double t = v_total.time();
518       
519    MyVector3 beta_v(0,0,z/t);
520
521    double E_CM = sqrt(s_hat);
522       
523    //##### BOTTOM PRODUCTION #####
524    double E_t_CM = E_CM/2;
525    double p_t_CM = sqrt(E_t_CM*E_t_CM - m_t4*m_t4);
526       
527    double phi = gRandom->Uniform(-M_PI,M_PI);
528       
529    MyVector3 vt31(p_t_CM*sin_theta*cos(phi),p_t_CM*sin_theta*sin(phi),p_t_CM*cos_theta);
530    MyVector3 vec0(0,0,0);
531    MyVector3 vt32 = vec0 - vt31;
532       
533    MyVector4 vt41(vt31,E_t_CM);
534    vt41.boost(beta_v);
535    MyVector4 vt42(vt32,E_t_CM);
536    vt42.boost(beta_v);
537       
538    Particle t4(vt41,7,2,0);//Change the PDG from 8 to 7 (t' to b')
539    t4.barcode(3);
540    t4.motherbarcode(1);
541    Event.push_back(t4);
542    Particle at4(vt42,-7,2,0);
543    at4.barcode(4);
544    at4.motherbarcode(1);
545    Event.push_back(at4);
546
547    //##### Generate full event and fill tree T_ana for Particle #####
548    if (i < 10)
549      cout << "*******************************************************" << endl;
550    for (int ip = 0; ip < Event.size(); ip++){
551      Particle p = Event.at(ip);
552      p.barcode(ip+1);
553
554      if (ip < m_nparmax){
555        m_pdgId[ip] = p.PDG();
556        m_stat[ip] = p.STAT();
557        m_barc[ip] = p.barcode();
558        m_phiP[ip] = p.vector4().phi();     
559        m_pTP[ip] = p.vector4().vector3().perp();
560        m_EP[ip] = p.vector4().E();
561        m_etaP[ip] = p.vector4().eta();
562      }
563      else
564        cout << "Too many particles !!! " << ip << endl;
565
566      if (i < 10)
567        cout << p ;
568      if (p.STAT() == 2){
569        Particle m = Event.at(Event.at(ip).motherbarcode()-1);
570        std::vector<Particle> p_dec = Decaybottom::DecayParticle(p,m);
571        cout << "this particle decays into : " << endl;
572        for (int j = 0; j<p_dec.size(); j++){
573          Particle p_d = p_dec.at(j);
574          Event.push_back(p_d);
575          cout << j+1 << " " << p_d;
576        }
577      }
578    }
579
580    if ( Event.size() <= m_nparmax)
581      m_npar = Event.size();
582    else {
583      cout << "Too many particles! " << Event.size() << endl;
584      m_npar = m_nparmax;
585    }
586
587    m_etaP[0] = 10.;
588    m_etaP[1] = -10.;
589
590    //######## Create vectors of electrons, jets and muons ########
591    int nele = 0, nmuo = 0, njet = 0;
592    double SumET = 0., ETmisX = 0., ETmisY = 0.;
593    double ETmisXmea = 0., ETmisYmea = 0.;
594    MyVector3 pMET(0,0,0);
595    MyVector3 pMETmea(0,0,0);
596    std::vector<Particle> electrons;
597    std::vector<Particle> jets;
598    std::vector<Particle> muons;
599
600    for (int j = 4; j < Event.size(); j++){
601      Particle pb = Event.at(j);
602
603      //vector of electrons
604      if (pb.STAT() == 1 && abs(pb.PDG()) == 11) {
605        electrons.push_back(pb);
606        jets.push_back(pb);   
607      }
608
609      //vector of jets
610      if (pb.STAT() == 1 && abs(pb.PDG()) <= 5)
611        jets.push_back(pb);
612
613      //vector of muons
614      if (pb.STAT() == 1 && abs(pb.PDG()) == 13)
615        muons.push_back(pb);
616
617     //neutrinos
618      if (abs(pb.PDG())==12 || abs(pb.PDG())==14){
619        MyVector3 v3pb = pb.vector4().vector3();
620        pMET += v3pb;
621      }
622
623      if (pb.STAT() == 1 && abs(pb.PDG())!=12 && abs(pb.PDG())!=14){
624        MyVector3 v3 = pb.vector4().vector3();
625        pMETmea -= v3;
626      }
627
628    }//END of <for (int j = 4; j < Event.size(); j++)>
629
630    //True missing energy calculated by using neutrinos
631    ETmisX = pMET.X();
632    ETmisY = pMET.Y();
633    ETmisY = pMET.Z();
634   
635    //Missing energy calculated by using all jets and muons
636    ETmisXmea = pMETmea.X();
637    ETmisYmea = pMETmea.Y();
638
639
640    //****** Fill tree T_ana : ELECTRONS ******
641    if (electrons.size()!=0){
642      for (int ie = 0; ie < electrons.size(); ie ++){
643        Particle pe = electrons.at(ie);
644        double eta = pe.vector4().eta();
645       
646        MyVector4 hlve = pe.vector4();
647        double E_e = hlve.E();
648       
649        if (fabs(hlve.eta()) < 5.) {
650          double ET = hlve.E()*sin(hlve.vector3().theta());
651          SumET += ET;
652        }
653
654        //Smear the energy of electrons
655        MyVector4 v4s = smear(pe);
656
657        // ... and fill with condition of eta
658        if (fabs(v4s.eta()) < 2.5){
659          if (ie < m_nelemax) {
660            m_pdgIdE[ie] = pe.PDG();
661            m_etaE[ie]   = v4s.eta();
662            m_phiE[ie]   = v4s.phi();
663            m_pTE[ie]    = v4s.Pt();
664            m_EE[ie]     = v4s.E();
665            nele++;
666          } else
667            cout << "Too many electrons ! Skipping..." << endl;
668        } 
669      }//END of <for (int ie = 0; ie < jets.size(); ie ++)>
670    }//END of <if (electrons.size()!=0)>
671
672    m_nele = nele;
673   
674    //******** Fill tree T_ana : JETS ********
675    if (jets.size()!=0){
676      for (int iq = 0; iq < jets.size(); iq ++){
677        Particle pq = jets.at(iq);
678        double eta = pq.vector4().eta();
679       
680        MyVector4 hlvq = pq.vector4();
681        double E_q = hlvq.E();
682       
683        if (fabs(hlvq.eta()) < 5.) {
684          double ET = hlvq.E()*sin(hlvq.vector3().theta());
685          SumET += ET;
686        }
687
688        //smear the energy of jets
689        MyVector4 v4s = smear(pq);
690
691        // ... and fill with condition of eta
692        if (fabs(v4s.eta()) < 5.0){
693          if (iq < m_nqmax){
694            m_pdgIdQ[iq] = pq.PDG();
695            m_etaQ[iq]   = v4s.eta();
696            m_phiQ[iq]   = v4s.phi();
697            m_pTQ[iq]    = v4s.Pt();
698            m_EQ[iq]     = v4s.E();
699            njet++;
700          } else
701            cout << "Too many jets ! Skipping..." << endl;
702        }
703      }//END of <for (int iq = 0; iq < jets.size(); iq ++)>
704    }//END of <if (jets.size()!=0)>
705
706    m_nq = njet;
707   
708    //******** Fill tree T_ana : MUONS ********
709    if (muons.size()!=0){
710      for (int im = 0; im < muons.size(); im ++){
711        Particle pm = muons.at(im);
712        double eta = pm.vector4().eta();
713       
714        MyVector4 hlvm = pm.vector4();
715        double pT_m = hlvm.Pt();
716        double E_m = hlvm.E();
717       
718        //smear the transverse momentum of muons
719        MyVector4 v4s = smear(pm);
720
721        // ... and fill with condition of eta
722        if (fabs(v4s.eta()) < 2.8){
723          if (im < m_nmumax){
724            m_pdgIdM[im] = pm.PDG();
725            m_etaM[im]   = v4s.eta();
726            m_phiM[im]   = v4s.phi();
727            m_pTM[im]    = v4s.Pt();
728            m_EM[im]     = v4s.E();
729            nmuo++;
730          } else
731            cout << "Too many muons ! Skipping..." << endl;
732        }
733      }//END of <for (int im = 0; im < muons.size(); im ++)>
734    }//END of <if (muons.size()!=0)>
735
736    m_nmu = nmuo;
737
738    //*** Fill tree T_ana : Missing Energy after smearing ***
739    double sigma_MET = 0.45*sqrt(SumET);
740    m_ETmisX = gRandom->Gaus(ETmisX,sigma_MET);
741    m_ETmisY = gRandom->Gaus(ETmisY,sigma_MET);
742    m_ETmisXmea = gRandom->Gaus(ETmisXmea,sigma_MET);
743    m_ETmisYmea = gRandom->Gaus(ETmisYmea,sigma_MET);
744    m_SumET  = gRandom->Gaus(SumET,sigma_MET);
745
746    m_evt = i;
747    resT->Fill();
748    resT_ana->Fill();
749   
750  }//END for (int i = 0; i < N; i++)
751 
752  resT->Write();
753  resT_ana->Write();
754  resF->Close();
755
756  return 0;
757}
758
759#include "LHAPDF/FortranWrappers.h"
760#ifdef FC_DUMMY_MAIN
761int FC_DUMMY_MAIN() {return l;}
762#endif
Note: See TracBrowser for help on using the repository browser.