ATLAS Offline Software
Loading...
Searching...
No Matches
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
17CLHEP::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
26double 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
236bool 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: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(double *params, double *x, bool useHelicity)
Definition UserSelections.h:26
convertVector
CLHEP::HepLorentzVector convertVector(const Pythia8::Vec4 v)
Definition UserSelections.h:17
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 &, std::string, std::vector< double >)
Definition UserSelections.h:236
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.14.0