ATLAS Offline Software
UserSelections.h
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // ======================================================================
6 // UserSelections
7 // James.Catmore@cern.ch
8 // This file replaces the old "user_finsel" where arbitrary selections
9 // can be made. There is a single bool method which takes
10 // (a) a string to point the code to the appropriate selection
11 // (b) a vector of doubles which can be used to specify cuts
12 // Additional code can be included, as with "BsJpsiPhiAngles.h" below
13 // ======================================================================
14 #include "CLHEP/Vector/LorentzVector.h"
15 #include "Pythia8B_i/BsJpsiPhiAngles.h"
16 
17 CLHEP::HepLorentzVector convertVector(const Pythia8::Vec4 v) {
18  CLHEP::HepLorentzVector vect;
19  vect.setPx(v.px());
20  vect.setPy(v.py());
21  vect.setPz(v.pz());
22  vect.setE(v.e());
23  return vect;
24 }
25 
26 double BsJpsiPhi_PDF(double *params, double *x, bool useHelicity) {
27 
28  // if ( !checkInput(x) ) return 0.;
29 
30  // Human-readable aliases of the parameters and observables.
31  // Human-readable aliases of the parameters and observables.
32  double A0 = std::sqrt(params[0] * (1 - params[2]));
33  double Al = std::sqrt(params[1] * (1 - params[2]));
34  double As = std::sqrt(params[2]);
35  double &GammaS = params[3];
36  double &DeltaGamma = params[4];
37  double &DeltaM = params[5];
38  double &phiS = params[6];
39  double &delta_p = params[7];
40  double &delta_l = params[8];
41  double delta_s = params[7] - params[9];
42 
43  double &time = x[0];
44  // double &timeErr = x[1]; // pico seconds
45  double &tagprob = x[4]; // 0. anti-particle, 0.5 no tag, 1.0 particle
46 
47  // Calculate A perpendicular from the other parameters
48  double Ap = 0;
49  if ( 1. - As*As - A0*A0 - Al*Al >= 0 ) Ap = std::sqrt(1. - As*As - A0*A0 - Al*Al);
50 
51  // The tabulated function is Sum of: A[i] * ( B1[i] +/- B2[i] ) * C[i]
52 
53  double A[10], B1[10], B2[10], C[10];
54 
55  // Define A cells
56  double sinphiS = std::sin(phiS);
57  double cosphiS = std::cos(phiS);
58 
59  A[0] = 0.5 * A0*A0;
60  A[1] = 0.5 * Al*Al;
61  A[2] = 0.5 * Ap*Ap;
62  A[3] = 0.5 * A0*Al * std::cos(delta_l);
63  A[4] = Al*Ap;
64  A[5] = A0*Ap;
65  A[6] = 0.5 * As*As;
66  A[7] = As*Al;
67  A[8] = 0.5 * As*Ap * std::sin(delta_p - delta_s);
68  A[9] = As*A0;
69 
70 
71 
72  // Define B cells
73  double gammaH = GammaS - DeltaGamma/2.;
74  double gammaL = GammaS + DeltaGamma/2.;
75  double tauL = 1./gammaL;
76  double tauH = 1./gammaH;
77  double norm1 = (1./(tauL*(1.+cosphiS) + tauH*(1.-cosphiS)));
78  double norm2 = (1./(tauL*(1.-cosphiS) + tauH*(1.+cosphiS)));
79  double norm3 = (1./sqrt((tauL-tauH)*(tauL-tauH)*sinphiS*sinphiS+4.*tauL*tauH));
80 
81  double ExpGSSinMT = 0.0;
82  double ExpGSCosMT = 0.0;
83 
84 
85 
86  // Calculate convolution of exp(-gammaL * time) with a gaussian. Not convoluted if sigma == 0.
87  double ExpGL = std::exp(-time * gammaL);
88 
89  // Calculate convolution of exp(-gammaH * time) with a gaussian. Not convoluted if sigma == 0.
90  double ExpGH = std::exp(-time * gammaH);
91 
92  B1[0] = ( (1. + cosphiS) * ExpGL + (1. - cosphiS) * ExpGH ) * norm1;
93  B1[1] = B1[0];
94  B1[2] = ( (1. - cosphiS) * ExpGL + (1. + cosphiS) * ExpGH ) * norm2;
95  B1[3] = B1[0];
96  B1[4] = 0.5 * std::cos(delta_p - delta_l) * sinphiS * (ExpGL - ExpGH) * norm3;
97  B1[5] = 0.5 * std::cos(delta_p) * sinphiS * (ExpGL - ExpGH) * norm3;
98  B1[6] = B1[2];
99  B1[7] = 0.5 * sinphiS * std::sin(delta_l - delta_s) * (ExpGL - ExpGH) * norm3;
100  B1[8] = B1[2];
101  B1[9] = 0.5 * sinphiS * std::sin(delta_s) * (ExpGH - ExpGL) * norm3;
102 
103 
104 
105 
106 
107  // Tagged analysis
108  if( std::abs(tagprob - 0.5) > 1e-6 ){
109 
110 
111  ExpGSSinMT = std::exp(-time * GammaS) * std::sin(DeltaM * time);
112  ExpGSCosMT = std::exp(-time * GammaS) * std::cos(DeltaM * time);
113 
114 
115  B2[0] = 2. * ExpGSSinMT * sinphiS * norm1;
116  B2[1] = B2[0];
117  B2[2] = -2. * ExpGSSinMT * sinphiS * norm2;
118  B2[3] = B2[0];
119  B2[4] = ( ExpGSCosMT * std::sin(delta_p - delta_l) - cosphiS * std::cos(delta_p - delta_l) * ExpGSSinMT ) * norm3;
120  B2[5] = ( ExpGSCosMT * std::sin(delta_p) - cosphiS * std::cos(delta_p) * ExpGSSinMT ) * norm3;
121  B2[6] = B2[2];
122  B2[7] = ( ExpGSCosMT * std::cos(delta_l - delta_s) - cosphiS * std::sin(delta_l - delta_s) * ExpGSSinMT ) * norm3;
123  B2[8] = B2[2];
124  B2[9] = ( ExpGSCosMT * std::cos(delta_s) + cosphiS * std::sin(delta_s) * ExpGSSinMT ) * norm3;
125 
126  } // End of tagged terms
127  else{
128  // Untagged analysis
129  for (int k = 0; k<10; k++) B2[k] = 0;
130 
131  }
132  // Define C cells
133 
134  if (useHelicity) {
135 
136  // Helicity Basis: Input variables are costhetal, costhetak and chi
137  double &costhetal = x[1];
138  double &costhetak = x[2];
139  double &chi = x[3];
140 
141  double sinsqthetal = 1. - (costhetal * costhetal);
142  double sinsqthetak = 1. - (costhetak * costhetak);
143  double cossqthetak = costhetak * costhetak;
144  double coschi = std::cos(chi);
145  double cossqchi = coschi * coschi;
146  double sinchi = std::sin(chi);
147  double sinsqchi = sinchi * sinchi;
148  double sin2thetak = 2. * std::sqrt(1. - costhetak * costhetak) * costhetak;
149  double sin2thetal = 2. * std::sqrt(1. - costhetal * costhetal) * costhetal;
150  double sin2chi = std::sin(2. * chi);
151  double sinthetak = std::sqrt(sinsqthetak);
152 
153  C[0] = 2. * cossqthetak * sinsqthetal; //cossqpsi * (1. - sinsqtheta * cossqphi);
154  C[1] = sinsqthetak * (1 - sinsqthetal * cossqchi); //sinsqpsi * (1. - sinsqtheta * sinsqphi);
155  C[2] = sinsqthetak * (1 - sinsqthetal * sinsqchi); //sinsqpsi * sinsqtheta;
156  C[3] = 1. / std::sqrt(2.) * sin2thetak * sin2thetal * coschi; //* sin2psi * sinsqtheta * sin2phi;
157  C[4] = -1. * sinsqthetak * sinsqthetal * sin2chi; //sinsqpsi * sin2theta * sinphi;
158  C[5] = 1. / std::sqrt(2.) * sin2thetak * sin2thetal * sinchi; //* sin2psi * sin2theta * cosphi;
159  C[6] = 2. / 3. * sinsqthetal; //(1. - sinsqtheta * cossqphi);
160  C[7] = 1. / 3. * std::sqrt(6.) * sinthetak * sin2thetal * coschi; //sinpsi * sinsqtheta * sin2phi;
161  C[8] = 1. / 3. * std::sqrt(6.) * sinthetak * sin2thetal * sinchi; //sinpsi * sin2theta * cosphi;
162  C[9] = 4. / 3. * std::sqrt(3.) * costhetak * sinsqthetal; // cospsi * (1. - sinsqtheta * cossqphi);
163 
164  }
165  else {
166 
167  // Transversity Basis: Input variables are costheta, cospsi and phi
168  double &costheta = x[1];
169  double &cospsi = x[2];
170  double &phi = x[3];
171 
172  double sinsqtheta = 1. - (costheta * costheta);
173  double sinsqpsi = 1. - (cospsi * cospsi);
174  double cossqpsi = cospsi * cospsi;
175  double sin2theta = 2. * std::sqrt(1. - costheta * costheta) * costheta;
176  double sin2psi = 2. * std::sqrt(1. - cospsi * cospsi) * cospsi;
177  double cosphi = std::cos(phi);
178  double cossqphi = cosphi * cosphi;
179  double sinphi = std::sin(phi);
180  double sinsqphi = sinphi * sinphi;
181  double sin2phi = std::sin(2. * phi);
182  double sinpsi = std::sqrt(1. - cospsi*cospsi);
183 
184  C[0] = 2. * cossqpsi * (1. - sinsqtheta * cossqphi);
185  C[1] = sinsqpsi * (1.-sinsqtheta * sinsqphi);
186  C[2] = sinsqpsi * sinsqtheta;
187  C[3] = 1. / sqrt(2.) * sin2psi * sinsqtheta * sin2phi;
188  C[4] = -1. * sinsqpsi * sin2theta * sinphi;
189  C[5] = 1. / sqrt(2.) * sin2psi * sin2theta * cosphi;
190  C[6] = 2. / 3. *(1. - sinsqtheta * cossqphi);
191  C[7] = 1. / 3. *sqrt(6.) *sinpsi * sinsqtheta * sin2phi;
192  C[8] = 1. / 3. *sqrt(6.) *sinpsi * sin2theta * cosphi;
193  C[9] = 4. / 3. *sqrt(3.) *cospsi * (1. - sinsqtheta * cossqphi);
194 
195  }
196 
197 
198 
199  // Sum the tabulated formula
200  double Wplus = 0;
201  double Wminus = 0;
202 
203 
204  for ( int i=0; i<10; ++i ) {
205 
206  Wplus += A[i] * ( B1[i] + B2[i] ) * C[i];
207  Wminus += A[i] * ( B1[i] - B2[i] ) * C[i];
208 
209  }
210 
211  double value = tagprob * Wplus + (1. - tagprob) * Wminus;
212  if (value < 0)
213  value = 0;
214  return value;
215 }
216 
217 //double transversityAnglePhi(const HepLorentzVector &jpsi,
218 // const HepLorentzVector &phi, const HepLorentzVector &kp,
219 // const HepLorentzVector &mup) {
220 //
221 // HepLorentzVector phi_jpsi = phi;
222 // HepLorentzVector kp_jpsi = kp;
223 // HepLorentzVector mup_jpsi = mup;
224 // phi_jpsi.boost(-jpsi.boostVector());
225 // kp_jpsi .boost(-jpsi.boostVector());
226 // mup_jpsi.boost(-jpsi.boostVector());
227 //
228 // Hep3Vector phi_direction = phi_jpsi.vect().unit();
229 // Hep3Vector phi_orthogonal = (kp_jpsi.vect().unit() - phi_direction
230 // * (kp_jpsi.vect().unit().dot(phi_direction))).unit();
231 //
232 // return atan2(mup_jpsi.vect().unit().dot(phi_orthogonal),
233 // mup_jpsi.vect().unit().dot(phi_direction));
234 //}
235 
236 bool Pythia8B_i::userSelection(Pythia8::Event &event, std::string userString,
237  std::vector<double> userVars) {
238 
239  using CLHEP::HepLorentzVector;
240  bool accept(false);
241 
242  // ==================================================
243  // code up your selection here; as many blocks as
244  // needed can be inserted
245  if (userString == "EXAMPLE" && event.size() > 0 && userVars.size() > 0) {
246 
247  accept = true;
248 
249  }
250 
251  // ==================================================
252  // selects decays where a B-hadron decays to a J/psi
253  // by any route
254  // Accepts event if event contains a B which (eventually)
255  // goes to a J/psi
256 
257  else if (userString == "BJPSIINCLUSIVE" && event.size() > 0) {
258  int chargeConj = 1;
259  if (userVars.size()==0) {
260  ATH_MSG_INFO("User selection BJPSIINCLUSIVE with B-state");
261  }
262  if (userVars.size()>0) {
263  if (userVars.size()>1) ATH_MSG_WARNING("User selection BJPSIINCLUSIVE with more than one argument! Check job options");
264  if (userVars.at(0)>0) ATH_MSG_DEBUG("User selection BJPSIINCLUSIVE with B-state");
265  if (userVars.at(0)<0) {ATH_MSG_DEBUG("User selection BJPSIINCLUSIVE with anti-B-state"); chargeConj = -1;}
266  }
267 
268  // B-decay codes which can go to charmonium
269  std::vector<int> bToCharmoniaCodes;
270  bToCharmoniaCodes.push_back(511*chargeConj);
271  bToCharmoniaCodes.push_back(521*chargeConj);
272  bToCharmoniaCodes.push_back(531*chargeConj);
273  bToCharmoniaCodes.push_back(541*chargeConj);
274  bToCharmoniaCodes.push_back(-5122*chargeConj);
275  bToCharmoniaCodes.push_back(-5132*chargeConj);
276  bToCharmoniaCodes.push_back(-5232*chargeConj);
277  bToCharmoniaCodes.push_back(-5332*chargeConj);
278 
279  int eventSize = event.size();
280 
281  bool isBtoJpsi(false);
282  bool containsB(false);
283  // loop over all particles in the event
284  for (int i = 0; i < eventSize; i++) {
285  int pdgID = event[i].id();
286  std::vector < Pythia8::Particle > decayMembers;
287  bool isB(false);
288  // Does event contain B which can go to J/psi?
289  for (unsigned int k = 0; k < bToCharmoniaCodes.size(); ++k) {
290  if (pdgID == bToCharmoniaCodes[k]) {
291  containsB = true;
292  isB = true;
293  break;
294  }
295  }
296  // Get decay chain of the B; see if there is a J/psi
297  if (isB) {
298  descendThroughDecay(event, decayMembers, i);
299  std::vector<int> pdgCodes = getCodes(decayMembers);
300  for (unsigned int k = 0; k < pdgCodes.size(); ++k) {
301  if (pdgCodes[k] == 443)
302  isBtoJpsi = true;
303  }
304  }
305  if (isBtoJpsi)
306  break;
307  }
308  if (containsB && isBtoJpsi)
309  accept = true;
310  if (containsB && !isBtoJpsi)
311  accept = false;
312  if (!containsB)
313  accept = false;
314 
315  }
316 
317  // ==================================================
318  // prototype for Bs->J/psi phi angular distribution
319  // (Adam Barton)
320 
321  else if ((userString == "BJPSIPHI_TRANS" || userString == "BJPSIPHI_HEL")
322  && event.size() > 0) {
323 
324  const bool debug = false;
325  bool flat;
326 
327  //Read
328  //userVars[0] 0(flat)/1(angle)
329  //[1] prob_limit
330  //[2] tag mode
331  //[3] A0
332  //[4] Al
333  //[5] As
334  //[6] GammaS
335  //[7] DeltaGamma
336  //[8] DeltaM
337  //[9] phiS
338  //[10] delta_p
339  //[11] delta_l
340  //[12]delta_s
341 
342  if (userVars.size() < 13) {
343  ATH_MSG_WARNING(
344  "REQUIRED userVARs not provided for BsJpsiPhi pdf defaulting to flat");
345  flat = true;
346  } else {
347  flat = userVars[0] == 0.0;
348 
349  }
350 
351  // for(int i=0; i<10;i++){
352  // ATH_MSG_INFO("BJPSIPHI_TRANS: " << i << " random number: " << Rdmengine->flat());
353  //
354  // }
355 
356  // std::vector<Pythia8::Particle> decayMembers;
357 
358  int i_Bs = 0;
359  int i_Muplus = 0;
360  int i_Muminus = 0;
361  int i_Kplus = 0;
362  int i_Kminus = 0;
363  int i_Phi = 0;
364  int i_Jpsi = 0;
365 
366  bool isBsJpsiPhi = false;
367  int eventSize = event.size();
368  for (int i = 0; i < eventSize; i++) {
369 
370  int pID = event[i].id();
371  if (std::abs(pID) == 531) { //NOTE THIS WILL FIND BS AND ANTIBS
372  i_Bs = i;
373  std::vector<int> daughterlist = event.daughterList(i);
374 
375  if (daughterlist.size() != 2)
376  continue;
377  bool isjpsi = false;
378  bool isphi = false;
379  if (event[daughterlist[0]].id() == 443) {
380  isjpsi = true;
381  i_Jpsi = daughterlist[0];
382  }
383  if (event[daughterlist[1]].id() == 443) {
384  isjpsi = true;
385  i_Jpsi = daughterlist[1];
386  }
387  if (event[daughterlist[0]].id() == 333) {
388  isphi = true;
389  i_Phi = daughterlist[0];
390  }
391  if (event[daughterlist[1]].id() == 333) {
392  isphi = true;
393  i_Phi = daughterlist[1];
394  }
395  if (!isphi || !isjpsi)
396  continue;
397  std::vector<int> daughterlistJpsi = event.daughterList(i_Jpsi);
398  std::vector<int> daughterlistPhi = event.daughterList(i_Phi);
399  if (daughterlistJpsi.size() != 2 || daughterlistPhi.size() != 2)
400  continue;
401  //resets values to avoid possible bug
402 
403  if (event[daughterlistJpsi[0]].id() == 13)
404  i_Muminus = daughterlistJpsi[0];
405  else if (event[daughterlistJpsi[1]].id() == 13)
406  i_Muminus = daughterlistJpsi[1];
407  else
408  i_Muminus = 0;
409 
410  if (event[daughterlistJpsi[0]].id() == -13)
411  i_Muplus = daughterlistJpsi[0];
412  else if (event[daughterlistJpsi[1]].id() == -13)
413  i_Muplus = daughterlistJpsi[1];
414  else
415  i_Muplus = 0;
416 
417  if (event[daughterlistPhi[0]].id() == 321)
418  i_Kplus = daughterlistPhi[0];
419  else if (event[daughterlistPhi[1]].id() == 321)
420  i_Kplus = daughterlistPhi[1];
421  else
422  i_Kplus = 0;
423 
424  if (event[daughterlistPhi[0]].id() == -321)
425  i_Kminus = daughterlistPhi[0];
426  else if (event[daughterlistPhi[1]].id() == -321)
427  i_Kminus = daughterlistPhi[1];
428  else
429  i_Kminus = 0;
430  if (i_Muminus == 0 || i_Muplus == 0 || i_Kminus == 0 || i_Kplus
431  == 0)
432  continue;
433  else {
434  if (debug)
435  ATH_MSG_INFO(
436  "found entire Bs->Jpsi(mumu)phi(KK) decay ");
437  isBsJpsiPhi = true;
438  break;
439  }
440  }
441  }
442 
443  if (!isBsJpsiPhi)
444  return false;
445  if (flat)
446  return true;
447 
448  Pythia8::Particle &p_Bs = event[i_Bs];
449  Pythia8::Particle &p_Muplus = event[i_Muplus];
450  Pythia8::Particle &p_Muminus = event[i_Muminus];
451  Pythia8::Particle &p_Kplus = event[i_Kplus];
452  Pythia8::Particle &p_Kminus = event[i_Kminus];
453  Pythia8::Particle &p_Phi = event[i_Phi];
454  Pythia8::Particle &p_Jpsi = event[i_Jpsi];
455 
456  // cout << "Bs " << p_Bs.id() << endl;
457  // cout << "|> " << p_Jpsi.id() << " + " << p_Phi.id() << endl;
458  // cout << "|> " << p_Muplus.id() << " + " << p_Muminus.id() << " + " << p_Kplus.id() << " + " << p_Kminus.id() << endl;
459 
460 
461  HepLorentzVector v_Bs = convertVector(p_Bs.p());
462  HepLorentzVector v_Muplus = convertVector(p_Muplus.p());
463  HepLorentzVector v_Muminus = convertVector(p_Muminus.p());
464  HepLorentzVector v_Kplus = convertVector(p_Kplus.p());
465  HepLorentzVector v_Kminus = convertVector(p_Kminus.p());
466  HepLorentzVector v_Phi = convertVector(p_Phi.p());
467  HepLorentzVector v_Jpsi = convertVector(p_Jpsi.p());
468 
469  BsJpsiPhiAngles angles(v_Kplus, v_Muplus, v_Phi, v_Jpsi, v_Bs);
470 
471  CLHEP::HepRandomEngine* Rdmengine = Pythia8B_i::p_rndmEngine;
472  const double gentau = Pythia8_i::m_pythia->particleData.tau0(531);
473  const double correctionfactor = 0.299792458;
474  const double gentauCorrect = gentau / correctionfactor;
475  if (debug) {
476  ATH_MSG_INFO("Lifetime for 531: " << gentau);
477  ATH_MSG_INFO("correct lifetime " << gentauCorrect);
478  }
479  const double Bstau = p_Bs.tau() / correctionfactor;
480  double prob1 = 1000;
481 
482  //PUT PDFS HERE
483  if (userString == "BJPSIPHI_TRANS") {
484  double x[5];
485  x[0] = Bstau;
486  x[1] = angles.thetatrfix();
487  x[2] = angles.thetaKfix();
488  x[3] = angles.phitrfix();
489  x[4] = userVars[2];
490 
491  prob1 = BsJpsiPhi_PDF(&userVars[3], x, false);
492  } else if (userString == "BJPSIPHI_HEL") {
493  double x[5];
494  x[0] = Bstau;
495  x[1] = angles.thetaLfix();
496  x[2] = angles.thetaKfix();
497  x[3] = angles.chifix();
498  x[4] = userVars[2];
499 
500  prob1 = BsJpsiPhi_PDF(&userVars[3], x, true);
501  }
502 
503  const double prob2 = exp(Bstau / gentauCorrect) * gentauCorrect;
504  if (Bstau > 20)
505  ATH_MSG_WARNING("Warning Bstau > 20 ps, this could indicate a bug");
506  const double prob_norm = userVars[1];
507  if (debug) {
508  ATH_MSG_INFO("prob limit set to " << prob_norm);
509  ATH_MSG_INFO("This Bs has lifetime " << Bstau);
510  }
511  double rand = Rdmengine->flat() * prob_norm;
512  if (prob1 * prob2 > prob_norm) {
513  ATH_MSG_WARNING(
514  "Max prob exceeded, too many of these indicates a problem, a few is fine");
515  }
516  if (debug)
517  std::cout << "totalprob " << prob1 * prob2 << std::endl;
518  if (rand < (prob1 * prob2)) {
519 
520  if (debug)
521  ATH_MSG_INFO("Passed PDF filter");
522  accept = true;
523  } else {
524  if (debug)
525  ATH_MSG_INFO("Rejected PDF filter");
526  accept = false;
527  }
528 
529  } // end of Bs->J/psiphi
530  else if(userString == "BD_BPLUS_TAUCUT"){
531  const double correctionfactor = 0.299792458;
532  bool pass = true;
533  int eventSize = event.size();
534  int foundcount=0;
535  for (int i = 0; i < eventSize; i++) {
536  int pID = event[i].id();
537  if (pID == 511) {
538  const double Bdtau = event[i].tau() / correctionfactor;
539  pass &= Bdtau > userVars[0] && Bdtau < userVars[1];
540  foundcount++;
541  }
542  if (pID == 521) {
543  const double Bptau = event[i].tau() / correctionfactor;
544  pass &= Bptau > userVars[0] && Bptau < userVars[1];
545  foundcount++;
546  }
547  }
548  accept = (foundcount > 0) & pass;
549  }
550 
551  else if ((userString == "BDJPSIKSTAR_TRANS") && event.size() > 0) {
552  const bool debug = false;
553  bool flat;
554  if (userVars.size() < 13) {
555  ATH_MSG_WARNING(
556  "REQUIRED userVARs not provided for BdJpsiKstar pdf defaulting to flat");
557  flat = true;
558  } else {
559  flat = userVars[0] == 0.0;
560  }
561 
562  int i_Bd = 0;
563  int i_Muplus = 0;
564  int i_Muminus = 0;
565  int i_Kplus = 0;
566  int i_piminus = 0;
567  int i_Kstar = 0;
568  int i_Jpsi = 0;
569 
570  bool isBsJpsiKstar = false;
571  const int eventSize = event.size();
572  for (int i = 0; i < eventSize; i++) {
573 
574  const int pID = event[i].id();
575  if (std::abs(pID) == 511) { //NOTE THIS FIND BD and Anti-Bd
576  i_Bd = i;
577  std::vector<int> daughterlist = event.daughterList(i);
578 
579  if (daughterlist.size() != 2)
580  continue;
581  bool isjpsi = false;
582  bool iskstar = false;
583  if (event[daughterlist[0]].id() == 443) {
584  isjpsi = true;
585  i_Jpsi = daughterlist[0];
586  }
587  if (event[daughterlist[1]].id() == 443) {
588  isjpsi = true;
589  i_Jpsi = daughterlist[1];
590  }
591  if (std::abs(event[daughterlist[0]].id()) == 313) { //This will find kstar or KstarBar
592  iskstar = true;
593  i_Kstar = daughterlist[0];
594  }
595  if (std::abs(event[daughterlist[1]].id()) == 313) { //This will find kstar or KstarBar
596  iskstar = true;
597  i_Kstar = daughterlist[1];
598  }
599  if (!iskstar || !isjpsi)
600  continue;
601  std::vector<int> daughterlistJpsi = event.daughterList(i_Jpsi);
602  std::vector<int> daughterlistKstar = event.daughterList(i_Kstar);
603  if (daughterlistJpsi.size() != 2 || daughterlistKstar.size() != 2)
604  continue;
605  //resets values to avoid possible bug
606 
607  if (event[daughterlistJpsi[0]].id() == 13)
608  i_Muminus = daughterlistJpsi[0];
609  else if (event[daughterlistJpsi[1]].id() == 13)
610  i_Muminus = daughterlistJpsi[1];
611  else
612  i_Muminus = 0;
613 
614  if (event[daughterlistJpsi[0]].id() == -13)
615  i_Muplus = daughterlistJpsi[0];
616  else if (event[daughterlistJpsi[1]].id() == -13)
617  i_Muplus = daughterlistJpsi[1];
618  else
619  i_Muplus = 0;
620 
621  if (std::abs(event[daughterlistKstar[0]].id()) == 321)
622  i_Kplus = daughterlistKstar[0];
623  else if (std::abs(event[daughterlistKstar[1]].id()) == 321)
624  i_Kplus = daughterlistKstar[1];
625  else
626  i_Kplus = 0;
627 
628  if (std::abs(event[daughterlistKstar[0]].id()) == 211)
629  i_piminus = daughterlistKstar[0];
630  else if (std::abs(event[daughterlistKstar[1]].id()) == 211)
631  i_piminus = daughterlistKstar[1];
632  else
633  i_piminus = 0;
634 
635  if (i_Muminus == 0 || i_Muplus == 0 || i_piminus == 0 || i_Kplus
636  == 0)
637  continue;
638 
639  if (debug)
640  ATH_MSG_INFO(
641  "found entire Bd->Jpsi(mumu)Kstar(KPi) decay ");
642  isBsJpsiKstar = true;
643  break;
644 
645  }
646  }
647 
648  if (!isBsJpsiKstar)
649  return false;
650  if (flat)
651  return true;
652 
653  Pythia8::Particle &p_Bd = event[i_Bd];
654  Pythia8::Particle &p_Muplus = event[i_Muplus];
655  Pythia8::Particle &p_Muminus = event[i_Muminus];
656  Pythia8::Particle &p_Kplus = event[i_Kplus];
657  //Pythia8::Particle &p_piminus = event[i_piminus];
658  Pythia8::Particle &p_Kstar = event[i_Kstar];
659  Pythia8::Particle &p_Jpsi = event[i_Jpsi];
660 
661  // cout << "Bs " << p_Bs.id() << endl;
662  // cout << "|> " << p_Jpsi.id() << " + " << p_Phi.id() << endl;
663  // cout << "|> " << p_Muplus.id() << " + " << p_Muminus.id() << " + " << p_Kplus.id() << " + " << p_Kminus.id() << endl;
664 
665 
666  HepLorentzVector v_Bd = convertVector(p_Bd.p());
667  HepLorentzVector v_Muplus = convertVector(p_Muplus.p());
668  HepLorentzVector v_Muminus = convertVector(p_Muminus.p());
669  HepLorentzVector v_Kplus = convertVector(p_Kplus.p());
670  //HepLorentzVector v_piminus = convertVector(p_piminus.p());
671  HepLorentzVector v_Kstar = convertVector(p_Kstar.p());
672  HepLorentzVector v_Jpsi = convertVector(p_Jpsi.p());
673 
674  BsJpsiPhiAngles angles(v_Kplus, v_Muplus, v_Kstar, v_Jpsi, v_Bd);
675 
676  CLHEP::HepRandomEngine* Rdmengine = Pythia8B_i::p_rndmEngine;
677  const double gentau = Pythia8_i::m_pythia->particleData.tau0(511);
678  const double correctionfactor = 0.299792458;
679  const double gentauCorrect = gentau / correctionfactor;
680  if (debug) {
681  ATH_MSG_INFO("Lifetime for 511: " << gentau);
682  ATH_MSG_INFO("correct lifetime " << gentauCorrect);
683  }
684  const double Bdtau = p_Bd.tau() / correctionfactor;
685  double prob1;
686 
687  //PUT PDFS HERE
688 
689  double x[5];
690  x[0] = Bdtau;
691  x[1] = angles.thetatrfix();
692  x[2] = angles.thetaKfix();
693  x[3] = angles.phitrfix();
694  x[4] = userVars[2];
695  prob1 = BsJpsiPhi_PDF(&userVars[3], x, false);
696 
697 
698  const double prob2 = exp(Bdtau / gentauCorrect) * gentauCorrect;
699  if (Bdtau > 20)
700  ATH_MSG_WARNING("Warning Bdtau > 20 ps, this could indicate a bug");
701  const double prob_norm = userVars[1];
702  if (debug) {
703  ATH_MSG_INFO("prob limit set to " << prob_norm);
704  ATH_MSG_INFO("This Bd has lifetime " << Bdtau);
705  }
706  const double rand = Rdmengine->flat() * prob_norm;
707  if (prob1 * prob2 > prob_norm) {
708  ATH_MSG_WARNING(
709  "Max prob exceeded, too many of these indicates a problem, a few is fine");
710  }
711  if (debug)
712  std::cout << "totalprob " << prob1 * prob2 << std::endl;
713  if (rand < (prob1 * prob2)) {
714  if (debug)
715  ATH_MSG_INFO("Passed PDF filter");
716  accept = true;
717  } else {
718  if (debug)
719  ATH_MSG_INFO("Rejected PDF filter");
720  accept = false;
721  }
722  }
723 
724  return (accept);
725 
726 }
727 
728 
phi
Scalar phi() const
phi method
Definition: AmgMatrixBasePlugin.h:64
ATH_MSG_INFO
#define ATH_MSG_INFO(x)
Definition: AthMsgStreamMacros.h:31
CutsMETMaker::accept
StatusCode accept(const xAOD::Muon *mu)
Definition: CutsMETMaker.cxx:18
DMTest::C
C_v1 C
Definition: C.h:26
JiveXML::Event
struct Event_t Event
Definition: ONCRPCServer.h:65
BsJpsiPhiAngles
Definition: BsJpsiPhiAngles.h:9
BsJpsiPhiAngles.h
athena.value
value
Definition: athena.py:122
convertVector
CLHEP::HepLorentzVector convertVector(const Pythia8::Vec4 v)
Definition: UserSelections.h:17
Pythia8B_i::getCodes
std::vector< int > getCodes(const std::vector< Pythia8::Particle > &) const
Definition: Pythia8B_i.cxx:632
drawFromPickle.cos
cos
Definition: drawFromPickle.py:36
drawFromPickle.exp
exp
Definition: drawFromPickle.py:36
BsJpsiPhiAngles::phitrfix
double phitrfix() const
Definition: BsJpsiPhiAngles.h:32
BsJpsiPhi_PDF
double BsJpsiPhi_PDF(double *params, double *x, bool useHelicity)
Definition: UserSelections.h:26
x
#define x
xAOD::Particle
Particle_v1 Particle
Define the latest version of the particle class.
Definition: Event/xAOD/xAODParticleEvent/xAODParticleEvent/Particle.h:17
Pythia8B_i::userSelection
bool userSelection(Pythia8::Event &, std::string, std::vector< double >)
Definition: UserSelections.h:236
dqt_zlumi_alleff_HIST.A
A
Definition: dqt_zlumi_alleff_HIST.py:110
BsJpsiPhiAngles::thetatrfix
double thetatrfix() const
Definition: BsJpsiPhiAngles.h:31
LArG4FSStartPointFilter.rand
rand
Definition: LArG4FSStartPointFilter.py:80
Pythia8B_i::descendThroughDecay
void descendThroughDecay(Pythia8::Event &, std::vector< Pythia8::Particle > &, int) const
Definition: Pythia8B_i.cxx:617
event
POOL::TEvent event(POOL::TEvent::kClassAccess)
lumiFormat.i
int i
Definition: lumiFormat.py:92
ATH_MSG_DEBUG
#define ATH_MSG_DEBUG(x)
Definition: AthMsgStreamMacros.h:29
debug
const bool debug
Definition: MakeUncertaintyPlots.cxx:53
BsJpsiPhiAngles::thetaKfix
double thetaKfix() const
Definition: BsJpsiPhiAngles.h:28
Pythia8B_i::p_rndmEngine
static CLHEP::HepRandomEngine * p_rndmEngine
Definition: Pythia8B_i.h:41
python.PyAthena.v
v
Definition: PyAthena.py:157
DiTauMassTools::MaxHistStrategyV2::e
e
Definition: PhysicsAnalysis/TauID/DiTauMassTools/DiTauMassTools/HelperFunctions.h:26
CaloSwCorrections.time
def time(flags, cells_name, *args, **kw)
Definition: CaloSwCorrections.py:242
ATH_MSG_WARNING
#define ATH_MSG_WARNING(x)
Definition: AthMsgStreamMacros.h:32
BsJpsiPhiAngles::thetaLfix
double thetaLfix() const
Definition: BsJpsiPhiAngles.h:29
drawFromPickle.sin
sin
Definition: drawFromPickle.py:36
PowhegControl_ttFCNC_NLO.params
params
Definition: PowhegControl_ttFCNC_NLO.py:226
Pythia8_i::m_pythia
std::unique_ptr< Pythia8::Pythia > m_pythia
Definition: Pythia8_i.h:88
fitman.k
k
Definition: fitman.py:528
BsJpsiPhiAngles::chifix
double chifix() const
Definition: BsJpsiPhiAngles.h:30
Generated on Tue Jul 2 2024 21:32:11 for ATLAS Offline Software by   doxygen 1.8.18