source: huonglan/MYSIMULATION/Saved/main_25_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.6 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"
7#include "Functions.C"
8#include "TRandom.h"
9#include "TH1F.h"
10#include "TF1.h"
11#include "TFile.h"
12#include "TTree.h"
13#include "LHAPDF/LHAPDF.h"
14
15using namespace std;
16
17TH1F *ctet_dis = new TH1F("ctet_dis","",100,-1.,1.);
18TH1F *ctet_dis_0 = new TH1F("ctet_dis_0","",100,-1.,1.);
19TH1F *ctet_dis_1 = new TH1F("ctet_dis_1","",100,-1.,1.);
20TF1 *f0 = new TF1("dis0","[0]*(1.- x*x)",-1,1); 
21TF1 *f1 = new TF1("dis0","[0]*(1.- x)*(1.- x)",-1,1); 
22TF1 *f = new TF1("f","(3./2.)*([0]*(1.-x)*(1.-x)/4. + [1]*(1.-x*x)/2.)");
23
24//histograms for top prime
25TH1F *pT_t4 = new TH1F("pT_t4","",100,0,2000);
26TH1F *ctet_t4 = new TH1F("cos_theta of t4","",100,-1,1);
27TH1F *eta_t4 = new TH1F("eta of t4","",100,0,10);
28TH1F *phi_t4 = new TH1F("phi of t4","",100,-M_PI,M_PI);
29
30//histogram for W
31TH1F *pT_W = new TH1F("pT_W","",100,0,1000);
32
33//histograms for charged lepton in W-decay
34TH1F *pT_lep = new TH1F("pT_lep","",100,0,1000);
35TH1F *eta_lep = new TH1F("eta_lep","",100,0,10);
36TH1F *phi_lep = new TH1F("phi_lep","",100,-M_PI,M_PI);
37
38//histograms for down quark in W-decay
39TH1F *pT_quark = new TH1F("pT_quark","",100,0,1000);
40TH1F *eta_quark = new TH1F("eta_quark","",100,0,10);
41TH1F *phi_quark = new TH1F("phi_quark","",100,-M_PI,M_PI);
42
43//histograms for the particles produced in W-decay
44TH1F *deltaR_ln = new TH1F("deltaR_ln","",100,-100,100);
45TH1F *deltaR_qq = new TH1F("deltaR_qq","",100,-100,100);
46
47//Reconstruct W and top
48TH1F *hW1 = new TH1F("hW1","",100,0,200);
49TH1F *hW2 = new TH1F("hW2","",100,0,200);
50TH1F *ht1 = new TH1F("ht1","",100,0,1000);
51TH1F *ht2 = new TH1F("ht2","",100,0,1000);
52
53
54/*
55    //######## FILL HISTO OF W ########
56    MyVector4 v4_W = a.vector4();
57    double WpT = v4_W.Pt();
58
59    pT_W -> Fill(WpT);
60
61    //#######################################
62
63
64    //######  FILL histograms for charged lepton (down quark) in W-decay ######
65    double leppT = vector4_data_b.Pt();
66    double lepeta = vector4_data_b.eta();
67    double lepphi = vector4_data_b.phi();
68
69    pT_lep -> Fill(leppT);
70    eta_lep -> Fill(lepeta);
71    phi_lep -> Fill(lepphi);
72   
73    //######  FILL histograms deltaR for particle produced in W-decay #######
74    double dR_W = vector4_data_b.deltaR(vector4_data_c);
75    cout << "Distance between the productions of W :" << dR_W << endl;
76
77    if (R < 0.216){
78      deltaR_ln -> Fill(dR_W);
79    }
80    else {
81       deltaR_qq -> Fill(dR_W);
82    }
83*/
84
85
86//############################################
87//************** MAIN PROGRAM ****************
88//############################################
89
90#define m_t4 400
91#define Ep 3500
92#define is_pp true
93
94int main()
95{
96  cout << endl;
97  double s = 4*Ep*Ep;
98  cout << endl;
99
100  // Init LHAPDF
101  const int SUBSET  = 0;
102  const std::string NAME = "MRST2007lomod";
103  LHAPDF::initPDFSet(NAME, LHAPDF::LHGRID, SUBSET);
104  const double Q = m_t4;
105  LHAPDF::initPDF(1);
106
107  std::vector<Particle> Event;
108
109  //##### Cross section #####
110  double N = 1000000;
111 
112  double tau_min = 4*m_t4*m_t4/s;
113  double y_min = log(sqrt(tau_min));
114  double y_max = -log(sqrt(tau_min));
115
116  double sigma_gg = 0.,
117    sigma_dd = 0., 
118    sigma_uu = 0.,
119    sigma_ss = 0.,
120    sigma_cc = 0.,
121    sigma_bb = 0.;
122 
123  double t_min = 1e30;
124  double t_max = - t_min;
125
126  TH1F *tau_gg = new TH1F("tau_gg","tau_gg", 1000, 0, 1);
127  TH1F *tau_qq = new TH1F("tau_qq","tau_qq", 1000, 0, 1);
128  TH1F *tau_unweigg = new TH1F("tau_unweigg","tau_unweigg", 1000, 0, 1);
129  TH1F *tau_unweiqq = new TH1F("tau_unweiqq","tau_unweiqq", 1000, 0, 1);
130
131  float m_tau, m_y, m_ctet_top, m_x, m_xf[11], m_wgg, m_wqq; 
132  TFile *resF = new TFile("MC_topprime.root","recreate");
133  TTree *resT = new TTree("T3","T3");
134  resT->Branch("tau",  &m_tau,  "tau/F");
135  resT->Branch("y",    &m_y,    "y/F");
136  resT->Branch("ctet_top", &m_ctet_top, "ctet_top/F");
137  resT->Branch("x",    &m_x,    "x/F");
138  resT->Branch("xf",   &m_xf,   "xf[11]/F");
139  resT->Branch("wgg",  &m_wgg,  "wgg/F");
140  resT->Branch("wqq",  &m_wqq,  "wqq/F");
141
142  double weight_ggmax = 1e-10;
143  double weight_qqmax = 1e-10;
144  double x1 = 0., x2 = 0.;
145
146
147  //############################################## 
148  // FIND THE WEIGHT_GGMAX AND WEIGHT_QQMAX VALUES
149  //############################################## 
150 
151  for (int i = 0; i < N; i++){
152
153    std::vector<double> vec = genTauYCos(Ep,m_t4);
154    double tau = vec[0];
155    double h_tau = vec[1];
156    double y = vec[2];
157    double h_y = vec[3];
158    double cos_theta = vec[4];
159    double h_cos = vec[5];
160    x1 = sqrt(tau)*exp(y);
161    x2 = sqrt(tau)*exp(-y);
162
163    double xf1_g = LHAPDF::xfx(x1,Q,0);
164    if (xf1_g < 0) xf1_g = 0;
165    double xf2_g = LHAPDF::xfx(x2,Q,0);
166    if (xf2_g < 0) xf2_g = 0;
167   
168    double xf1_d    = LHAPDF::xfx(x1,Q,1);
169    if (xf1_d < 0) xf1_d = 0;
170    double xf1_dbar = LHAPDF::xfx(x1,Q,-1);
171    if (xf1_dbar < 0) xf1_dbar = 0;
172    double xf2_d    = LHAPDF::xfx(x2,Q,1);
173    if (xf2_d < 0) xf2_d = 0;
174    double xf2_dbar = LHAPDF::xfx(x2,Q,-1);
175    if (xf2_dbar < 0) xf2_dbar = 0;
176   
177    double xf1_u    = LHAPDF::xfx(x1,Q,2);
178    if (xf1_u < 0) xf1_u = 0;
179    double xf1_ubar = LHAPDF::xfx(x1,Q,-2);
180    if (xf1_ubar < 0) xf1_ubar = 0;
181    double xf2_u    = LHAPDF::xfx(x2,Q,2);
182    if (xf2_u < 0) xf2_u = 0;
183    double xf2_ubar = LHAPDF::xfx(x2,Q,-2);
184    if (xf2_ubar < 0) xf2_ubar = 0;
185
186    double xf1_s    = LHAPDF::xfx(x1,Q,3);
187    if (xf1_s < 0) xf1_s = 0;
188    double xf2_sbar = LHAPDF::xfx(x2,Q,-3);
189    if (xf2_sbar < 0) xf2_sbar = 0;
190
191    double xf1_c    = LHAPDF::xfx(x1,Q,4);
192    if (xf1_c < 0) xf1_c = 0;
193    double xf2_cbar = LHAPDF::xfx(x2,Q,-4);
194    if (xf2_cbar < 0) xf2_cbar = 0;
195
196    double xf1_b    = LHAPDF::xfx(x1,Q,5);
197    if (xf1_b < 0) xf1_b = 0;
198    double xf2_bbar = LHAPDF::xfx(x2,Q,-5);
199    if (xf2_bbar < 0) xf2_bbar = 0;
200
201    double s_hat = x1*x2*s;
202    double beta = sqrt(1 - 4*m_t4*m_t4/s_hat);   
203    double u_hat = m_t4*m_t4 - s_hat/2.- beta*s_hat*cos_theta/2.;
204    double t_hat = m_t4*m_t4 - s_hat/2.+ beta*s_hat*cos_theta/2.;
205
206    //######## DSIGMAs & WEIGHTs ########
207    double dsigma_gg = dsigmaGG(s_hat,t_hat,u_hat,m_t4);
208    double dsigma_qq = dsigmaQQ(s_hat,t_hat,u_hat,m_t4);
209
210    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);
211 
212    //p-pbar case
213    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);
214    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);
215
216    //pp case
217    if (is_pp) 
218      {
219        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);
220        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);
221      }
222
223    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);
224    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);
225    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);
226    double weight_qq = weight_dd + weight_uu + weight_ss + weight_cc + weight_bb;
227
228    double r_gg = gRandom -> Uniform();
229    double r_qq = gRandom -> Uniform();
230
231    if (weight_gg > weight_ggmax)
232      weight_ggmax = weight_gg;
233   
234    if (weight_qq > weight_qqmax)
235      weight_qqmax = weight_qq;
236
237    sigma_gg += weight_gg;
238   
239    sigma_dd += weight_dd;
240    sigma_uu += weight_uu;
241    sigma_ss += weight_ss;
242    sigma_cc += weight_cc;
243    sigma_bb += weight_bb;
244
245  }
246
247  /*cout << "*********************************" << endl; 
248  cout << "weight_ggmax : " << weight_ggmax << endl;
249  cout << "weight_qqmax : " << weight_qqmax << endl;
250  cout << "*********************************" << endl; 
251
252  cout << "t_min : " << t_min << "   t_max : " << t_max << endl;*/
253 
254  double dtau = 1 - tau_min;
255  double dy = y_max - y_min;
256  double dcost = 2.;
257
258  // Here we dont have to multiply sigma
259  // *** with the phase-space volume ***
260  sigma_gg = M_PI*sigma_gg/N/s;
261  sigma_dd = M_PI*sigma_dd/N/s;
262  sigma_uu = M_PI*sigma_uu/N/s;
263  sigma_ss = M_PI*2*sigma_ss/N/s;
264  sigma_cc = M_PI*2*sigma_cc/N/s;
265  sigma_bb = M_PI*2*sigma_bb/N/s;
266
267  double sigma_qq = sigma_dd + sigma_uu + sigma_ss + sigma_cc + sigma_bb;
268
269  //### CONVERT THE UNIT OF CROSS SECTION TO PB
270
271  double sigma_convert1 = sigma_gg*0.389*1e9;
272  double sigma_convert2 = sigma_dd*0.389*1e9;
273  double sigma_convert3 = sigma_uu*0.389*1e9;
274  double sigma_convert4 = sigma_ss*0.389*1e9;
275  double sigma_convert5 = sigma_cc*0.389*1e9;
276  double sigma_convert6 = sigma_bb*0.389*1e9;
277  double sigma_convert7 = sigma_qq*0.389*1e9;
278
279  /*cout << "After convert the unit of cross section: " << endl;
280
281  cout << "gg->t't'bar      : " << sigma_convert1 << endl << endl;
282  cout << "d-dbar->t't'bar  : " << sigma_convert2 << endl;
283  cout << "u-ubar->t't'bar  : " << sigma_convert3 << endl;
284  cout << "s-sbar->t't'bar  : " << sigma_convert4 << endl;
285  cout << "c-cbar->t't'bar  : " << sigma_convert5 << endl;
286  cout << "b-bbar->t't'bar  : " << sigma_convert6 << endl << endl;
287  cout << "Total of the processes q-qbar->t't'bar  : " << sigma_convert7 << endl << endl;
288
289  cout << "TOTAL : " << sigma_convert1 + sigma_convert7 << endl;   
290  */
291
292
293  //######################################
294  //######## GENERATE REAL EVENTS ########
295  //######################################
296
297  double N1 = 1000000;
298  double N_gg = 0.;
299  double N_qq = 0.;
300
301  for (int i = 0; i < N1; i++){
302    if (i%100 == 0)
303      cout << "generating event " << i << endl;
304
305    //##################################
306    std::vector<double> vec = genTauYCos(3500,400);
307    double tau = vec[0];
308    double h_tau = vec[1];
309    double y = vec[2];
310    double h_y = vec[3];
311    double cos_theta = vec[4];
312    double sin_theta = sqrt(1-cos_theta*cos_theta);
313    double h_cos = vec[5];
314
315    //########## X1 & X2 ##########
316    x1 = sqrt(tau)*exp(y);
317    x2 = sqrt(tau)*exp(-y);
318
319    double s_hat = x1*x2*s;
320
321    double beta = sqrt(1 - 4*m_t4*m_t4/s_hat);   
322    double u_hat = m_t4*m_t4 - s_hat/2.- beta*s_hat*cos_theta/2.;
323    double t_hat = m_t4*m_t4 - s_hat/2.+ beta*s_hat*cos_theta/2.;
324   
325    //fill the tree
326    m_tau   = tau;
327    m_y     = y;
328    m_ctet_top = cos_theta;
329    m_x     = x1;
330
331    std::vector<Particle> Event;
332    MyVector4 vp1,vp2;
333       
334    double r_gg = gRandom -> Uniform();
335    double k_gg = sigma_convert1/(sigma_convert1 + sigma_convert7);
336           
337    if (r_gg < k_gg){
338      //cout << "gg -> QQbar" << endl;
339      double xf1_g = LHAPDF::xfx(x1,Q,0);
340      if (xf1_g < 0) xf1_g = 0;
341      double xf2_g = LHAPDF::xfx(x2,Q,0);
342      if (xf2_g < 0) xf2_g = 0;
343
344      double dsigma_gg = dsigmaGG(s_hat,t_hat,u_hat,m_t4);     
345      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);
346
347      //fill the tree & histos
348      m_xf[0] = xf1_g;
349      m_xf[1] = xf2_g;
350      m_wgg = weight_gg;
351      tau_gg -> Fill(tau,weight_gg);     
352
353      double r = gRandom -> Uniform();
354   
355      if (weight_gg/k_gg < r*(weight_ggmax/k_gg + weight_qqmax/(1.- k_gg))) continue;
356      else {
357        N_gg ++;
358        tau_unweigg -> Fill(tau);
359       
360        MyVector4 vp11(MyVector3(0,0,x1*Ep),x1*Ep);
361        vp1 = vp11;
362        Particle p1(vp1,21,3,0);
363        p1.barcode(1);
364        p1.motherbarcode(0);
365        MyVector4 vp22(MyVector3(0,0,- x2*Ep),x2*Ep);
366        vp2 = vp22;
367        Particle p2(vp2,21,3,0);
368        p2.barcode(2);
369        p2.motherbarcode(0);
370       
371        Event.push_back(p1);
372        Event.push_back(p2);
373      }
374    }//END of <if (r_gg < k_gg)>
375
376    else {
377      //cout << "qqbar -> QQbar" << endl;
378      double xf1_d    = LHAPDF::xfx(x1,Q,1);
379      if (xf1_d < 0) xf1_d = 0;
380      double xf1_dbar = LHAPDF::xfx(x1,Q,-1);
381      if (xf1_dbar < 0) xf1_dbar = 0;
382      double xf2_d    = LHAPDF::xfx(x2,Q,1);
383      if (xf2_d < 0) xf2_d = 0;
384      double xf2_dbar = LHAPDF::xfx(x2,Q,-1);
385      if (xf2_dbar < 0) xf2_dbar = 0;
386     
387      double xf1_u    = LHAPDF::xfx(x1,Q,2);
388      if (xf1_u < 0) xf1_u = 0;
389      double xf1_ubar = LHAPDF::xfx(x1,Q,-2);
390      if (xf1_ubar < 0) xf1_ubar = 0;
391      double xf2_u    = LHAPDF::xfx(x2,Q,2);
392      if (xf2_u < 0) xf2_u = 0;
393      double xf2_ubar = LHAPDF::xfx(x2,Q,-2);
394      if (xf2_ubar < 0) xf2_ubar = 0;
395     
396      double xf1_s    = LHAPDF::xfx(x1,Q,3);
397      if (xf1_s < 0) xf1_s = 0;
398      double xf2_sbar = LHAPDF::xfx(x2,Q,-3);
399      if (xf2_sbar < 0) xf2_sbar = 0;
400     
401      double xf1_c    = LHAPDF::xfx(x1,Q,4);
402      if (xf1_c < 0) xf1_c = 0;
403      double xf2_cbar = LHAPDF::xfx(x2,Q,-4);
404      if (xf2_cbar < 0) xf2_cbar = 0;
405     
406      double xf1_b    = LHAPDF::xfx(x1,Q,5);
407      if (xf1_b < 0) xf1_b = 0;
408      double xf2_bbar = LHAPDF::xfx(x2,Q,-5);
409      if (xf2_bbar < 0) xf2_bbar = 0;
410
411      double dsigma_qq = dsigmaQQ(s_hat,t_hat,u_hat,m_t4);
412
413      //p-pbar case
414      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);
415      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);
416
417      //pp case
418      if (is_pp) {
419        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);
420        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);
421      }
422
423      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);
424      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);
425      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);
426      double weight_qq = weight_dd + weight_uu + weight_ss + weight_cc + weight_bb;
427
428      //Fill tree and histo
429      m_xf[2] = xf1_d;
430      m_xf[3] = xf2_dbar;
431      m_xf[4] = xf1_u;
432      m_xf[5] = xf2_ubar;
433      m_xf[6] = xf1_s;
434      m_xf[7] = xf2_sbar;
435      m_xf[8] = xf1_c;
436      m_xf[9] = xf2_cbar;
437      m_xf[10] = xf1_b;
438      m_xf[11] = xf2_bbar;
439      m_wqq = weight_qq;   
440      tau_qq -> Fill(tau,weight_qq);//histo
441
442      double r = gRandom -> Uniform();
443
444      if (weight_qq/(1.- k_gg) < r*(weight_ggmax/k_gg + weight_qqmax/(1.- k_gg))) continue;
445      else {
446        N_qq ++;
447        tau_unweiqq -> Fill(tau);
448       
449        MyVector4 vp11(MyVector3(0,0,x1*Ep),x1*Ep);
450        vp1 = vp11;
451        Particle p1(vp1,2,1,0);
452        p1.barcode(1);
453        p1.motherbarcode(0);
454        MyVector4 vp22(MyVector3(0,0,- x2*Ep),x2*Ep);
455        vp2 = vp22;
456        Particle p2(vp2,-2,1,0);
457        p2.barcode(2);
458        p2.motherbarcode(0);
459       
460        Event.push_back(p1);
461        Event.push_back(p2);
462      }
463
464    }//END of else of <if (r_gg < k_gg)>
465
466    MyVector4 v_total = vp1 + vp2;
467    double z = v_total.vector3().Z();
468    double t = v_total.time();
469       
470    MyVector3 beta_v(0,0,z/t);
471
472    double E_CM = sqrt(s_hat);
473       
474    //for top production:
475    double E_t_CM = E_CM/2;
476    double p_t_CM = sqrt(E_t_CM*E_t_CM - m_t4*m_t4);
477       
478    double phi = gRandom->Uniform(-M_PI,M_PI);
479       
480    MyVector3 vt31(p_t_CM*sin_theta*cos(phi),p_t_CM*sin_theta*sin(phi),p_t_CM*cos_theta);
481    MyVector3 vec0(0,0,0);
482    MyVector3 vt32 = vec0 - vt31;
483       
484    MyVector4 vt41(vt31,E_t_CM);
485    vt41.boost(beta_v);
486    MyVector4 vt42(vt32,E_t_CM);
487    vt42.boost(beta_v);
488       
489    Particle t4(vt41,8,2,0);
490    t4.barcode(3);
491    t4.motherbarcode(1);
492    Event.push_back(t4);
493    Particle at4(vt42,-8,2,0);
494    at4.barcode(4);
495    at4.motherbarcode(1);
496    Event.push_back(at4);
497
498    //######## FILL HISTOS OF TOP PRIME ########
499    double pT_top = vt41.Pt();
500    double eta_top = vt41.eta();
501    //double phi_top = vt41.phi();
502
503    pT_t4 -> Fill(pT_top);
504    ctet_t4 -> Fill(cos_theta);
505    eta_t4 -> Fill(eta_top);
506    phi_t4 -> Fill(phi);
507    //##########################################
508
509    cout << "*******************************************************" << endl;
510    for (int i = 0; i < Event.size(); i++){
511      Particle p = Event.at(i);
512      p.barcode(i+1);
513      cout << p ;
514      if (p.STAT() == 2){
515        Particle m = Event.at(Event.at(i).motherbarcode()-1);
516        std::vector<Particle> p_dec = Decay::DecayParticle(p,m);
517        for (int j = 0; j<p_dec.size(); j++){
518          Particle p_d = p_dec.at(j);
519          Event.push_back(p_d);
520        }
521      }
522    }
523
524    //#############  TEST  ################
525    cout << endl;
526    cout << "##############  TEST INITIAL STATE AND FINAL STATE  ##############" << endl << endl;
527       
528    MyVector4 v4_test(MyVector3(0,0,0),0);
529    for (int k = 2; k < Event.size(); k++){
530      Particle p_test = Event.at(k);
531      if (p_test.STAT() == 1){
532        MyVector4 v4 = p_test.vector4();
533        v4_test += v4;
534      }
535    }
536       
537    cout << "Initial state : " << v_total << endl;
538    cout << "Final state : " << v4_test << endl;
539 
540    //############ BUILT W AND TOP-PRIME ############
541    cout << " ####### BUILT W AND TOP ########" << endl;
542
543    MyVector4 vW1, vW2, vt1, vt2;
544    double mW1, mW2, mt1, mt2;
545
546    //Create the vectors including leptons and quarks
547    std::vector<Particle> lepW;
548    std::vector<Particle> quarks;
549
550    for (int j = 2; j < Event.size(); j++){
551      Particle pb = Event.at(j);
552      pb.barcode(j+1);
553 
554      if (abs(pb.PDG())>=11 && abs(pb.PDG())<=14)
555        lepW.push_back(pb);
556
557      if (abs(pb.PDG())==1||abs(pb.PDG())==2)
558        quarks.push_back(pb);
559
560    }
561
562    //Reconstruct leptonic mass of W
563    //First we check if there are leptons in the decays
564    //Then we choose the charged and neutral leptons (i: charged; j: neutral)
565    //Applying the cut for mass of W with <if (a > -0.1 && a < 0.1)>
566
567    if (lepW.size() != 0){
568      for (int ilep = 0; ilep < lepW.size()-1; ilep ++){
569        for (int jlep = lepW.size()-1 ; jlep>0 ; jlep --){
570
571          Particle lepi = lepW.at(ilep);
572          Particle lepj = lepW.at(jlep);
573                 
574          if (lepi.PDG()%2 !=0){
575            MyVector4 vi = lepi.vector4();
576         
577            if (lepj.PDG()%2 ==0){
578              MyVector4 vj = lepj.vector4();
579
580              MyVector4 vW1 = vi + vj;
581              mW1 = vW1.mass();
582              double dW = mW1 - 80.4;
583              if (dW > -0.1 && dW < 0.1){
584                double pTW = vW1.vector3().perp();
585                double pTlep = vi.vector3().perp();
586                pT_W -> Fill(pTW);
587                pT_lep -> Fill(pTlep);
588                //cout << "lepi : " << lepi ;
589                //cout << "lepj : " << lepj ;
590                hW1 -> Fill(mW1);
591             
592                //Reconstruct leptonic mass of top-prime
593                //Combine these reconstructed masses after the cut
594                //with the quarks from top decays
595                for (int iq = 0; iq < quarks.size(); iq ++){
596                  Particle qi = quarks.at(iq);
597                  if (qi.motherbarcode()!=5 && qi.motherbarcode()!=7){
598                  MyVector4 viq = qi.vector4();
599                  MyVector4 vt1 = vW1 + viq;
600                  mt1 = vt1.mass();
601                  double dt = mt1 - 400;
602                  if (dt > -0.1 && dt < 0.1){
603                    cout << "qi : " << qi ;
604                    cout << "W : " << vW1 << endl;
605                    cout << "t : " << vt1 << endl;
606                    ht1 -> Fill(mt1);
607
608                    //###Calculate cos_theta_star from vt1,vW1,vi
609                    //boost of W and boost of t <-- use to deboost the particles
610                    MyVector3 bW(-vW1.vector3().X()/vW1.time(),-vW1.vector3().Y()/vW1.time(),-vW1.vector3().Z()/vW1.time());
611                    cout << "bW : " << bW << endl;
612                    //direction of down quark in CM of W
613                    vi.boost(bW);
614                    cout << "lepi after boost : " << vi;
615                    MyVector3 vi3 = vi.vector3();
616                    double vi3_mag = vi3.mag();         
617
618                    MyVector3 bt(-vt1.vector3().X()/vt1.time(),-vt1.vector3().Y()/vt1.time(),-vt1.vector3().Z()/vt1.time());
619                    //direction of W in top' rest frame
620                    vW1.boost(bt);
621                    cout << "W after boost : " << vW1;
622                    MyVector3 vW13 = vW1.vector3();
623                    double vW13_mag = vW13.mag();
624                                   
625                    double ctet_star = (vW13*vi3)*(1./vW13_mag/vi3_mag);
626                    ctet_dis -> Fill(ctet_star);
627                    //cout << "cos_theta : " << ctet_star << endl;
628                  }}}}}}}}
629    }// END of <Reconstruct leptonic mass of W and top'>
630
631
632    //Reconstruct hadronic mass of W
633    if (quarks.size() != 0){
634      for (int iq = 0; iq < quarks.size()-1; iq ++){
635        for (int jq = quarks.size()-1 ; jq>iq ; jq --){
636          Particle qi = quarks.at(iq);
637          MyVector4 vi = qi.vector4();
638          Particle qj = quarks.at(jq);
639          MyVector4 vj = qj.vector4();
640         
641          MyVector4 vW2 = vi + vj;
642          mW2 = vW2.mass();
643          double dW = mW2 - 80.4;
644          if (dW > -0.1 && dW < 0.1){
645            double pTW = vW2.vector3().perp();
646            double pTquark1 = vi.vector3().perp();
647            double pTquark2 = vj.vector3().perp();
648            pT_W -> Fill(pTW);
649            pT_quark -> Fill(pTquark1);
650            pT_quark -> Fill(pTquark2);
651            hW2 -> Fill(mW2);
652
653            //Reconstruct hadronic mass of top-prime
654            for (int nq = 0; nq < quarks.size(); nq ++){
655              Particle qn = quarks.at(nq);
656              if (qn.motherbarcode()!=5 && qn.motherbarcode()!=7){
657                MyVector4 vqn = qn.vector4();
658                MyVector4 vt2 = vW2 + vqn;
659                mt2 = vt2.mass();
660                double dt = mt2 - 400;
661                if (dt > -0.1 && dt < 0.1){
662                    ht2 -> Fill(mt2);
663                }}}}}}
664    }// END of <Reconstruct hadronic mass of W and top'>
665
666    resT->Fill();
667   
668  }//END for (int i = 0; i < N; i++)
669 
670  resT->Write();
671  tau_gg -> Write();
672  tau_qq -> Write();
673  tau_unweigg -> Write();
674  tau_unweiqq -> Write();
675  ctet_dis -> Write();
676  ctet_dis_0 -> Write();
677  ctet_dis_1 -> Write();
678
679  //fill histos of top prime
680  pT_t4 -> Write();
681  ctet_t4 -> Write();
682  eta_t4 -> Write();
683  phi_t4 -> Write();
684
685  //fill histo of W
686  pT_W -> Write();
687
688  //fill histos of charged lepton (down quark) produced in W-decay
689  pT_lep -> Write();
690  eta_lep -> Write();
691  phi_lep -> Write();
692  pT_quark -> Write();
693  eta_quark -> Write();
694  phi_quark -> Write();
695
696  //fill histos deltaR of particles produced in W-decay
697  deltaR_ln -> Write();
698  deltaR_qq -> Write();
699
700  //fill histos for W and t reconstructed mass
701  hW1 -> Write();
702  hW2 -> Write();
703  ht1 -> Write();
704  ht2 -> Write();
705
706  resF->Close();
707
708  cout << "number of events for gg : " << N_gg << endl;
709  cout << "number of events for qq : " << N_qq << endl;
710
711  cout << "***********************************" << endl << endl;
712
713  return 0;
714}
715
716
717#include "LHAPDF/FortranWrappers.h"
718#ifdef FC_DUMMY_MAIN
719int FC_DUMMY_MAIN() {return l;}
720#endif
Note: See TracBrowser for help on using the repository browser.