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