evtgen is hosted by Hepforge, IPPP Durham
EvtGen  2.0.0
Monte Carlo generator of particle decays, in particular the weak decays of heavy flavour particles such as B mesons.
EvtLambdaB2LambdaV.cpp
Go to the documentation of this file.
1 
2 /***********************************************************************
3 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
4 * *
5 * This file is part of EvtGen. *
6 * *
7 * EvtGen is free software: you can redistribute it and/or modify *
8 * it under the terms of the GNU General Public License as published by *
9 * the Free Software Foundation, either version 3 of the License, or *
10 * (at your option) any later version. *
11 * *
12 * EvtGen is distributed in the hope that it will be useful, *
13 * but WITHOUT ANY WARRANTY; without even the implied warranty of *
14 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
15 * GNU General Public License for more details. *
16 * *
17 * You should have received a copy of the GNU General Public License *
18 * along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
19 ***********************************************************************/
20 
21 #include "EvtGenModels/EvtLambdaB2LambdaV.hh"
22 
23 #include "EvtGenBase/EvtGenKine.hh"
24 #include "EvtGenBase/EvtPDL.hh"
25 #include "EvtGenBase/EvtParticle.hh"
26 #include "EvtGenBase/EvtPatches.hh"
27 #include "EvtGenBase/EvtRandom.hh"
28 #include "EvtGenBase/EvtReport.hh"
29 
30 #include <fstream>
31 #include <stdio.h>
32 #include <stdlib.h>
33 #include <string>
34 
35 using std::fstream;
36 //************************************************************************
37 //* *
38 //* Class EvtLambdaB2LambdaV *
39 //* *
40 //************************************************************************
41 //DECLARE_COMPONENT( EvtLambdaB2LambdaV );
42 
43 EvtLambdaB2LambdaV::EvtLambdaB2LambdaV()
44 {
45  //set facility name
46  fname = "EvtGen.EvtLambdaB2LambdaV";
47 }
48 
49 //------------------------------------------------------------------------
50 // Method 'getName'
51 //------------------------------------------------------------------------
52 std::string EvtLambdaB2LambdaV::getName()
53 {
54  return "LAMBDAB2LAMBDAV";
55 }
56 
57 //------------------------------------------------------------------------
58 // Method 'clone'
59 //------------------------------------------------------------------------
60 EvtDecayBase* EvtLambdaB2LambdaV::clone()
61 {
62  return new EvtLambdaB2LambdaV;
63 }
64 
65 //------------------------------------------------------------------------
66 // Method 'initProbMax'
67 //------------------------------------------------------------------------
68 void EvtLambdaB2LambdaV::initProbMax()
69 {
70  //maximum (case where C=0)
71  double Max = 1 + fabs( A * B );
72  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
73  << " PDF max value : " << Max << std::endl;
74  setProbMax( Max );
75 }
76 
77 //------------------------------------------------------------------------
78 // Method 'init'
79 //------------------------------------------------------------------------
80 void EvtLambdaB2LambdaV::init()
81 {
82  bool antiparticle = false;
83 
84  //introduction
85  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
86  << "*************************************************" << std::endl;
87  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
88  << "* Event Model Class : EvtLambdaB2LambdaV *" << std::endl;
89  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
90  << "*************************************************" << std::endl;
91 
92  //check the number of arguments
93  checkNArg( 2 );
94 
95  //check the number of daughters
96  checkNDaug( 2 );
97 
98  const EvtId Id_mother = getParentId();
99  const EvtId Id_daug1 = getDaug( 0 );
100  const EvtId Id_daug2 = getDaug( 1 );
101 
102  //lambdab identification
103  if ( Id_mother == EvtPDL::getId( "Lambda_b0" ) ) {
104  antiparticle = false;
105  } else if ( Id_mother == EvtPDL::getId( "anti-Lambda_b0" ) ) {
106  antiparticle = true;
107  } else {
108  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
109  << " Mother is not a Lambda_b0 or an anti-Lambda_b0, but a "
110  << EvtPDL::name( Id_mother ) << std::endl;
111  abort();
112  }
113 
114  //lambda
115  if ( !( Id_daug1 == EvtPDL::getId( "Lambda0" ) && !antiparticle ) &&
116  !( Id_daug1 == EvtPDL::getId( "anti-Lambda0" ) && antiparticle ) ) {
117  if ( !antiparticle ) {
118  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
119  << " Daughter1 is not a Lambda0, but a "
120  << EvtPDL::name( Id_daug1 ) << std::endl;
121  } else {
122  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
123  << " Daughter1 is not an anti-Lambda0, but a "
124  << EvtPDL::name( Id_daug1 ) << std::endl;
125  }
126  abort();
127  }
128 
129  //identification meson V
130  if ( getArg( 1 ) == 1 )
131  Vtype = VID::JPSI;
132  else if ( getArg( 1 ) == 2 )
133  Vtype = VID::RHO;
134  else if ( getArg( 1 ) == 3 )
135  Vtype = VID::OMEGA;
136  else if ( getArg( 1 ) == 4 )
137  Vtype = VID::RHO_OMEGA_MIXING;
138  else {
139  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
140  << " Vtype " << getArg( 1 ) << " is unknown" << std::endl;
141  abort();
142  }
143 
144  //check vector meson V
145  if ( Id_daug2 == EvtPDL::getId( "J/psi" ) && Vtype == VID::JPSI ) {
146  if ( !antiparticle )
147  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
148  << " Decay mode successfully initialized : Lambda_b0 -> Lambda J/psi"
149  << std::endl;
150  else
151  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
152  << " Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda J/psi"
153  << std::endl;
154  } else if ( Id_daug2 == EvtPDL::getId( "rho0" ) && Vtype == VID::RHO ) {
155  if ( !antiparticle )
156  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
157  << " Decay mode successfully initialized : Lambda_b0 -> Lambda rho0"
158  << std::endl;
159  else
160  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
161  << " Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda rho0"
162  << std::endl;
163  } else if ( Id_daug2 == EvtPDL::getId( "omega" ) && Vtype == VID::OMEGA ) {
164  if ( !antiparticle )
165  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
166  << " Decay mode successfully initialized : Lambda_b0 -> Lambda omega"
167  << std::endl;
168  else
169  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
170  << " Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda omega"
171  << std::endl;
172  } else if ( ( Id_daug2 == EvtPDL::getId( "omega" ) ||
173  Id_daug2 == EvtPDL::getId( "rho0" ) ) &&
174  Vtype == VID::RHO_OMEGA_MIXING ) {
175  if ( !antiparticle )
176  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
177  << " Decay mode successfully initialized : "
178  << "Lambda_b0 -> Lambda rho-omega-mixing" << std::endl;
179  else
180  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
181  << " Decay mode successfully initialized : "
182  << "anti-Lambda_b0 -> anti-Lambda rho-omega-mixing" << std::endl;
183  }
184 
185  else {
186  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
187  << " Daughter2 is not a J/psi, phi or rho0 but a "
188  << EvtPDL::name( Id_daug2 ) << std::endl;
189  abort();
190  }
191 
192  //fix dynamics parameters
193  B = (double)getArg( 0 );
194  C = EvtComplex( ( sqrt( 2. ) / 2. ), ( sqrt( 2. ) / 2. ) );
195  switch ( Vtype ) {
196  case VID::JPSI:
197  A = 0.490;
198  break;
199  case VID::RHO:
200  case VID::OMEGA:
201  case VID::RHO_OMEGA_MIXING:
202  A = 0.194;
203  break;
204  default:
205  A = 0;
206  break;
207  }
208  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
209  << " LambdaB decay parameters : " << std::endl;
210  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
211  << " - lambda asymmetry A = " << A << std::endl;
212  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
213  << " - lambdab polarisation B = " << B << std::endl;
214  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
215  << " - lambdab density matrix rho+- C = " << C << std::endl;
216 }
217 
218 //------------------------------------------------------------------------
219 // Method 'decay'
220 //------------------------------------------------------------------------
221 void EvtLambdaB2LambdaV::decay( EvtParticle* lambdab )
222 {
223  lambdab->makeDaughters( getNDaug(), getDaugs() );
224 
225  //get lambda and V particles
226  EvtParticle* lambda = lambdab->getDaug( 0 );
227  EvtParticle* V = lambdab->getDaug( 1 );
228 
229  //get resonance masses
230  // - LAMBDAB -> mass given by the preceding class
231  // - LAMBDA -> nominal mass
232  // - V -> getVmass method
233  double MASS_LAMBDAB = lambdab->mass();
234  double MASS_LAMBDA = EvtPDL::getMeanMass( EvtPDL::getId( "Lambda0" ) );
235  double MASS_V = getVMass( MASS_LAMBDAB, MASS_LAMBDA );
236 
237  //generate random angles
238  double phi = EvtRandom::Flat( 0, 2 * EvtConst::pi );
239  double theta = acos( EvtRandom::Flat( -1, +1 ) );
240  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
241  << " Angular angles : theta = " << theta << " ; phi = " << phi
242  << std::endl;
243  //computate resonance quadrivectors
244  double E_lambda = ( MASS_LAMBDAB * MASS_LAMBDAB +
245  MASS_LAMBDA * MASS_LAMBDA - MASS_V * MASS_V ) /
246  ( 2 * MASS_LAMBDAB );
247  double E_V = ( MASS_LAMBDAB * MASS_LAMBDAB + MASS_V * MASS_V -
248  MASS_LAMBDA * MASS_LAMBDA ) /
249  ( 2 * MASS_LAMBDAB );
250  double P = sqrt( E_lambda * E_lambda - lambda->mass() * lambda->mass() );
251 
252  EvtVector4R P_lambdab = lambdab->getP4();
253 
254  double px = P_lambdab.get( 1 );
255  double py = P_lambdab.get( 2 );
256  double pz = P_lambdab.get( 3 );
257  double E = P_lambdab.get( 0 );
258  EvtGenReport( EVTGEN_INFO, fname.c_str() )
259  << "E of lambdab: " << P_lambdab.get( 0 ) << std::endl;
260  EvtGenReport( EVTGEN_INFO, fname.c_str() )
261  << "E of lambdab: " << E << std::endl;
262 
263  EvtVector4R q_lambdab2( E,
264  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
265  ( ( px * ( px ) ) + ( py * ( py ) ) ) ),
266  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
267  ( -( ( py ) * ( px ) ) + ( px * ( py ) ) ) ),
268  ( pz ) );
269 
270  EvtVector4R q_lambdab3( E, q_lambdab2.get( 3 ), q_lambdab2.get( 1 ),
271  q_lambdab2.get( 2 ) );
272 
273  EvtVector4R q_lambda0( E_lambda, P * sin( theta ) * cos( phi ),
274  P * sin( theta ) * sin( phi ), P * cos( theta ) );
275 
276  EvtVector4R q_V0( E_V, -P * sin( theta ) * cos( phi ),
277  -P * sin( theta ) * sin( phi ), -P * cos( theta ) );
278 
279  EvtVector4R q_lambda1( E_lambda, q_lambda0.get( 2 ), q_lambda0.get( 3 ),
280  q_lambda0.get( 1 ) );
281 
282  EvtVector4R q_V1( E_V, q_V0.get( 2 ), q_V0.get( 3 ), q_V0.get( 1 ) );
283 
284  EvtVector4R q_lambda(
285  E_lambda,
286  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
287  ( ( px * ( q_lambda1.get( 1 ) ) ) - ( py * ( q_lambda1.get( 2 ) ) ) ) ),
288  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
289  ( ( py * ( q_lambda1.get( 1 ) ) ) + ( px * ( q_lambda1.get( 2 ) ) ) ) ),
290  ( q_lambda1.get( 3 ) ) );
291 
292  EvtVector4R q_V(
293  E_V,
294  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
295  ( ( px * ( q_V1.get( 1 ) ) ) - ( py * ( q_V1.get( 2 ) ) ) ) ),
296  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
297  ( ( py * ( q_V1.get( 1 ) ) ) + ( px * ( q_V1.get( 2 ) ) ) ) ),
298  ( q_V1.get( 3 ) ) );
299 
300  lambda->getP4();
301  V->getP4();
302  EvtGenReport( EVTGEN_INFO, fname.c_str() )
303  << " LambdaB px: " << px << std::endl;
304  EvtGenReport( EVTGEN_INFO, fname.c_str() )
305  << " LambdaB py: " << py << std::endl;
306  EvtGenReport( EVTGEN_INFO, fname.c_str() )
307  << " LambdaB pz: " << pz << std::endl;
308  EvtGenReport( EVTGEN_INFO, fname.c_str() )
309  << " LambdaB E: " << E << std::endl;
310 
311  EvtGenReport( EVTGEN_INFO, fname.c_str() )
312  << " Lambdab3 E: " << q_lambdab3.get( 0 ) << std::endl;
313  EvtGenReport( EVTGEN_INFO, fname.c_str() )
314  << " Lambda 0 px: " << q_lambda0.get( 1 ) << std::endl;
315  EvtGenReport( EVTGEN_INFO, fname.c_str() )
316  << " Lambda 0 py: " << q_lambda0.get( 2 ) << std::endl;
317  EvtGenReport( EVTGEN_INFO, fname.c_str() )
318  << " Lambda 0 pz: " << q_lambda0.get( 3 ) << std::endl;
319  EvtGenReport( EVTGEN_INFO, fname.c_str() )
320  << " Lambda 0 E: " << q_lambda0.get( 0 ) << std::endl;
321  EvtGenReport( EVTGEN_INFO, fname.c_str() )
322  << " Lambda 1 px: " << q_lambda1.get( 1 ) << std::endl;
323  EvtGenReport( EVTGEN_INFO, fname.c_str() )
324  << " Lambda 1 py: " << q_lambda1.get( 2 ) << std::endl;
325  EvtGenReport( EVTGEN_INFO, fname.c_str() )
326  << " Lambda 1 pz: " << q_lambda1.get( 3 ) << std::endl;
327  EvtGenReport( EVTGEN_INFO, fname.c_str() )
328  << " Lambda 1 E: " << q_lambda1.get( 0 ) << std::endl;
329  EvtGenReport( EVTGEN_INFO, fname.c_str() )
330  << " Lambda px: " << q_lambda.get( 1 ) << std::endl;
331  EvtGenReport( EVTGEN_INFO, fname.c_str() )
332  << " Lambda py: " << q_lambda.get( 2 ) << std::endl;
333  EvtGenReport( EVTGEN_INFO, fname.c_str() )
334  << " Lambda pz: " << q_lambda.get( 3 ) << std::endl;
335  EvtGenReport( EVTGEN_INFO, fname.c_str() )
336  << " Lambda E: " << q_lambda0.get( 3 ) << std::endl;
337  EvtGenReport( EVTGEN_INFO, fname.c_str() )
338  << " V 0 px: " << q_V0.get( 1 ) << std::endl;
339  EvtGenReport( EVTGEN_INFO, fname.c_str() )
340  << " V 0 py: " << q_V0.get( 2 ) << std::endl;
341  EvtGenReport( EVTGEN_INFO, fname.c_str() )
342  << " V 0 pz: " << q_V0.get( 3 ) << std::endl;
343  EvtGenReport( EVTGEN_INFO, fname.c_str() )
344  << " V 0 E: " << q_V0.get( 0 ) << std::endl;
345  EvtGenReport( EVTGEN_INFO, fname.c_str() )
346  << " V 1 px: " << q_V1.get( 1 ) << std::endl;
347  EvtGenReport( EVTGEN_INFO, fname.c_str() )
348  << " V 1 py: " << q_V1.get( 2 ) << std::endl;
349  EvtGenReport( EVTGEN_INFO, fname.c_str() )
350  << " V 1 pz: " << q_V1.get( 3 ) << std::endl;
351  EvtGenReport( EVTGEN_INFO, fname.c_str() )
352  << " V 1 E: " << q_V1.get( 0 ) << std::endl;
353  EvtGenReport( EVTGEN_INFO, fname.c_str() )
354  << " V px: " << q_V.get( 1 ) << std::endl;
355  EvtGenReport( EVTGEN_INFO, fname.c_str() )
356  << " V py: " << q_V.get( 2 ) << std::endl;
357  EvtGenReport( EVTGEN_INFO, fname.c_str() )
358  << " V pz: " << q_V.get( 3 ) << std::endl;
359  EvtGenReport( EVTGEN_INFO, fname.c_str() )
360  << " V E: " << q_V0.get( 3 ) << std::endl;
361  //set quadrivectors to particles
362  lambda->init( getDaugs()[0], q_lambda );
363  V->init( getDaugs()[1], q_V );
364 
365  //computate pdf
366  double pdf = 1 + A * B * cos( theta ) +
367  2 * A * real( C * EvtComplex( cos( phi ), sin( phi ) ) ) *
368  sin( theta );
369 
370  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
371  << " LambdaB decay pdf value : " << pdf << std::endl;
372  //set probability
373  setProb( pdf );
374 
375  return;
376 }
377 
378 //------------------------------------------------------------------------
379 // Method 'BreitWignerRelPDF'
380 //------------------------------------------------------------------------
381 double EvtLambdaB2LambdaV::BreitWignerRelPDF( double m, double _m0, double _g0 )
382 {
383  double a = ( _m0 * _g0 ) * ( _m0 * _g0 );
384  double b = ( m * m - _m0 * _m0 ) * ( m * m - _m0 * _m0 );
385  return a / ( b + a );
386 }
387 
388 //------------------------------------------------------------------------
389 // Method 'RhoOmegaMixingPDF'
390 //------------------------------------------------------------------------
391 double EvtLambdaB2LambdaV::RhoOmegaMixingPDF( double m, double _mr, double _gr,
392  double _mo, double _go )
393 {
394  double a = m * m - _mr * _mr;
395  double b = m * m - _mo * _mo;
396  double c = _gr * _mr;
397  double d = _go * _mo;
398  double re_pi = -3500e-6; //GeV^2
399  double im_pi = -300e-6; //GeV^2
400  double va_pi = re_pi + im_pi;
401 
402  //computate pdf value
403  double f = 1 / ( a * a + c * c ) *
404  ( 1 + ( va_pi * va_pi + 2 * b * re_pi + 2 * d * im_pi ) /
405  ( b * b + d * d ) );
406 
407  //computate pdf max value
408  a = 0;
409  b = _mr * _mr - _mo * _mo;
410 
411  double maxi = 1 / ( a * a + c * c ) *
412  ( 1 + ( va_pi * va_pi + 2 * b * re_pi + 2 * d * im_pi ) /
413  ( b * b + d * d ) );
414 
415  return f / maxi;
416 }
417 
418 //------------------------------------------------------------------------
419 // Method 'getVMass'
420 //------------------------------------------------------------------------
421 double EvtLambdaB2LambdaV::getVMass( double MASS_LAMBDAB, double MASS_LAMBDA )
422 {
423  //JPSI case
424  if ( Vtype == VID::JPSI ) {
425  return EvtPDL::getMass( EvtPDL::getId( "J/psi" ) );
426  }
427 
428  //RHO,OMEGA,MIXING case
429  else {
430  //parameters
431  double MASS_RHO = EvtPDL::getMeanMass( EvtPDL::getId( "rho0" ) );
432  double MASS_OMEGA = EvtPDL::getMeanMass( EvtPDL::getId( "omega" ) );
433  double WIDTH_RHO = EvtPDL::getWidth( EvtPDL::getId( "rho0" ) );
434  double WIDTH_OMEGA = EvtPDL::getWidth( EvtPDL::getId( "omega" ) );
435  double MASS_PION = EvtPDL::getMeanMass( EvtPDL::getId( "pi-" ) );
436  double _max = MASS_LAMBDAB - MASS_LAMBDA;
437  double _min = 2 * MASS_PION;
438 
439  double mass = 0;
440  bool test = false;
441  int ntimes = 10000;
442  do {
443  double y = EvtRandom::Flat( 0, 1 );
444 
445  //generate mass
446  mass = ( _max - _min ) * EvtRandom::Flat( 0, 1 ) + _min;
447 
448  //pdf calculate
449  double f = 0;
450  switch ( Vtype ) {
451  case VID::RHO:
452  f = BreitWignerRelPDF( mass, MASS_RHO, WIDTH_RHO );
453  break;
454  case VID::OMEGA:
455  f = BreitWignerRelPDF( mass, MASS_OMEGA, WIDTH_OMEGA );
456  break;
457  case VID::RHO_OMEGA_MIXING:
458  f = RhoOmegaMixingPDF( mass, MASS_RHO, WIDTH_RHO,
459  MASS_OMEGA, WIDTH_OMEGA );
460  break;
461  default:
462  f = 1;
463  break;
464  }
465 
466  ntimes--;
467  if ( y < f )
468  test = true;
469  } while ( ntimes && !test );
470 
471  //looping 10000 times
472  if ( ntimes == 0 ) {
473  EvtGenReport( EVTGEN_INFO, fname.c_str() )
474  << "Tried accept/reject:10000"
475  << " times, and rejected all the times!" << std::endl;
476  EvtGenReport( EVTGEN_INFO, fname.c_str() )
477  << "Is therefore accepting the last event!" << std::endl;
478  }
479 
480  //return V particle mass
481  return mass;
482  }
483 }
484 
485 //************************************************************************
486 //* *
487 //* Class EvtLambda2PPiForLambdaB2LambdaV *
488 //* *
489 //************************************************************************
490 
491 //------------------------------------------------------------------------
492 // Constructor
493 //------------------------------------------------------------------------
494 EvtLambda2PPiForLambdaB2LambdaV::EvtLambda2PPiForLambdaB2LambdaV()
495 {
496  //set facility name
497  fname = "EvtGen.EvtLambda2PPiForLambdaB2LambdaV";
498 }
499 
500 //------------------------------------------------------------------------
501 // Method 'getName'
502 //------------------------------------------------------------------------
503 std::string EvtLambda2PPiForLambdaB2LambdaV::getName()
504 {
505  return "LAMBDA2PPIFORLAMBDAB2LAMBDAV";
506 }
507 
508 //------------------------------------------------------------------------
509 // Method 'clone'
510 //------------------------------------------------------------------------
511 EvtDecayBase* EvtLambda2PPiForLambdaB2LambdaV::clone()
512 {
513  return new EvtLambda2PPiForLambdaB2LambdaV;
514 }
515 
516 //------------------------------------------------------------------------
517 // Method 'initProbMax'
518 //------------------------------------------------------------------------
519 void EvtLambda2PPiForLambdaB2LambdaV::initProbMax()
520 {
521  //maximum (case where D is real)
522  double Max = 0;
523  if ( A == 0 )
524  Max = 1;
525  else if ( C == 0 || real( D ) == 0 )
526  Max = 1 + fabs( A * B );
527  else if ( B == 0 )
528  Max = 1 + EvtConst::pi / 2.0 * fabs( C * A * real( D ) );
529  else {
530  double theta_max = atan( -EvtConst::pi / 2.0 * C * real( D ) / B );
531  double max1 = 1 + fabs( A * B * cos( theta_max ) -
532  EvtConst::pi / 2.0 * C * A * real( D ) *
533  sin( theta_max ) );
534  double max2 = 1 + fabs( A * B );
535  if ( max1 > max2 )
536  Max = max1;
537  else
538  Max = max2;
539  }
540  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
541  << " PDF max value : " << Max << std::endl;
542  setProbMax( Max );
543 }
544 
545 //------------------------------------------------------------------------
546 // Method 'init'
547 //------------------------------------------------------------------------
548 void EvtLambda2PPiForLambdaB2LambdaV::init()
549 {
550  bool antiparticle = false;
551 
552  //introduction
553  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
554  << " ***********************************************************"
555  << std::endl;
556  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
557  << " * Event Model Class : EvtLambda2PPiForLambdaB2LambdaV *"
558  << std::endl;
559  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
560  << " ***********************************************************"
561  << std::endl;
562 
563  //check the number of arguments
564  checkNArg( 2 );
565 
566  //check the number of daughters
567  checkNDaug( 2 );
568 
569  const EvtId Id_mother = getParentId();
570  const EvtId Id_daug1 = getDaug( 0 );
571  const EvtId Id_daug2 = getDaug( 1 );
572 
573  //lambda identification
574  if ( Id_mother == EvtPDL::getId( "Lambda0" ) ) {
575  antiparticle = false;
576  } else if ( Id_mother == EvtPDL::getId( "anti-Lambda0" ) ) {
577  antiparticle = true;
578  } else {
579  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
580  << " Mother is not a Lambda0 or an anti-Lambda0, but a "
581  << EvtPDL::name( Id_mother ) << std::endl;
582  abort();
583  }
584 
585  //proton identification
586  if ( !( Id_daug1 == EvtPDL::getId( "p+" ) && !antiparticle ) &&
587  !( Id_daug1 == EvtPDL::getId( "anti-p-" ) && antiparticle ) ) {
588  if ( !antiparticle ) {
589  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
590  << " Daughter1 is not a p+, but a " << EvtPDL::name( Id_daug1 )
591  << std::endl;
592  } else {
593  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
594  << " Daughter1 is not an anti-p-, but a "
595  << EvtPDL::name( Id_daug1 ) << std::endl;
596  }
597  abort();
598  }
599 
600  //pion identification
601  if ( !( Id_daug2 == EvtPDL::getId( "pi-" ) && !antiparticle ) &&
602  !( Id_daug2 == EvtPDL::getId( "pi+" ) && antiparticle ) ) {
603  if ( !antiparticle ) {
604  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
605  << " Daughter2 is not a p-, but a " << EvtPDL::name( Id_daug1 )
606  << std::endl;
607  } else {
608  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
609  << " Daughter2 is not an p+, but a " << EvtPDL::name( Id_daug1 )
610  << std::endl;
611  }
612  abort();
613  }
614  if ( !antiparticle )
615  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
616  << " Decay mode successfully initialized : Lambda0 -> p+ pi-"
617  << std::endl;
618  else
619  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
620  << " Decay mode successfully initialized : Anti-Lambda0 -> anti-p- pi+"
621  << std::endl;
622 
623  //identification meson V
624  if ( getArg( 1 ) == 1 ) {
625  Vtype = VID::JPSI;
626  if ( !antiparticle )
627  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
628  << " From : Lambda_b0 -> Lambda J/psi" << std::endl;
629  else
630  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
631  << " From : anti-Lambda_b0 -> anti-Lambda J/psi" << std::endl;
632  } else if ( getArg( 1 ) == 2 ) {
633  Vtype = VID::RHO;
634  if ( !antiparticle )
635  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
636  << " From : Lambda_b0 -> Lambda rho0" << std::endl;
637  else
638  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
639  << " From : anti-Lambda_b0 -> anti-Lambda rho0" << std::endl;
640  } else if ( getArg( 1 ) == 3 ) {
641  Vtype = VID::OMEGA;
642  if ( !antiparticle )
643  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
644  << " From : Lambda_b0 -> Lambda omega" << std::endl;
645  else
646  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
647  << " From : anti-Lambda_b0 -> anti-Lambda omega" << std::endl;
648  } else if ( getArg( 1 ) == 4 ) {
649  Vtype = VID::RHO_OMEGA_MIXING;
650  } else {
651  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
652  << " Vtype " << getArg( 1 ) << " is unknown" << std::endl;
653  if ( !antiparticle ) {
654  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
655  << " From : Lambda_b0 -> Lambda rho-omega-mixing" << std::endl;
656  } else {
657  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
658  << " From : anti-Lambda_b0 -> anti-Lambda rho-omega-mixing"
659  << std::endl;
660  }
661  abort();
662  }
663 
664  //constants
665  A = 0.642;
666  C = (double)getArg( 0 );
667  switch ( Vtype ) {
668  case VID::JPSI:
669  B = -0.167;
670  D = EvtComplex( 0.25, 0.0 );
671  break;
672  case VID::RHO:
673  case VID::OMEGA:
674  case VID::RHO_OMEGA_MIXING:
675  B = -0.21;
676  D = EvtComplex( 0.31, 0.0 );
677  break;
678  default:
679  B = 0;
680  D = EvtComplex( 0, 0 );
681  break;
682  }
683 
684  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
685  << " Lambda decay parameters : " << std::endl;
686  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
687  << " - proton asymmetry A = " << A << std::endl;
688  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
689  << " - lambda polarisation B = " << B << std::endl;
690  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
691  << " - lambdaB polarisation C = " << C << std::endl;
692  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
693  << " - lambda density matrix rho+- D = " << D << std::endl;
694 }
695 
696 //------------------------------------------------------------------------
697 // Method 'decay'
698 //------------------------------------------------------------------------
699 void EvtLambda2PPiForLambdaB2LambdaV::decay( EvtParticle* lambda )
700 {
701  lambda->makeDaughters( getNDaug(), getDaugs() );
702 
703  //get proton and pion particles
704  EvtParticle* proton = lambda->getDaug( 0 );
705  EvtParticle* pion = lambda->getDaug( 1 );
706 
707  //get resonance masses
708  // - LAMBDA -> mass given by EvtLambdaB2LambdaV class
709  // - PROTON & PION -> nominal mass
710  double MASS_LAMBDA = lambda->mass();
711  double MASS_PROTON = EvtPDL::getMeanMass( EvtPDL::getId( "p+" ) );
712  double MASS_PION = EvtPDL::getMeanMass( EvtPDL::getId( "pi-" ) );
713 
714  //generate random angles
715  double phi = EvtRandom::Flat( 0, 2 * EvtConst::pi );
716  double theta = acos( EvtRandom::Flat( -1, +1 ) );
717  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
718  << " Angular angles : theta = " << theta << " ; phi = " << phi
719  << std::endl;
720 
721  //computate resonance quadrivectors
722  double E_proton = ( MASS_LAMBDA * MASS_LAMBDA + MASS_PROTON * MASS_PROTON -
723  MASS_PION * MASS_PION ) /
724  ( 2 * MASS_LAMBDA );
725  double E_pion = ( MASS_LAMBDA * MASS_LAMBDA + MASS_PION * MASS_PION -
726  MASS_PROTON * MASS_PROTON ) /
727  ( 2 * MASS_LAMBDA );
728  double P = sqrt( E_proton * E_proton - proton->mass() * proton->mass() );
729 
730  EvtVector4R P_lambda = lambda->getP4();
731  EvtParticle* Mother_lambda = lambda->getParent();
732  EvtVector4R lambdab = Mother_lambda->getP4();
733 
734  double E_lambda = P_lambda.get( 0 );
735  double px_M = lambdab.get( 1 );
736  double py_M = lambdab.get( 2 );
737  double pz_M = lambdab.get( 3 );
738  double E_M = lambdab.get( 0 );
739 
740  EvtVector4R q_lambdab2(
741  E_M,
742  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
743  ( ( px_M * ( px_M ) ) + ( py_M * ( py_M ) ) ) ),
744  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
745  ( -( ( py_M ) * ( px_M ) ) + ( px_M * ( py_M ) ) ) ),
746  ( pz_M ) );
747 
748  EvtVector4R q_lambdab3( E_M, q_lambdab2.get( 3 ), q_lambdab2.get( 1 ),
749  q_lambdab2.get( 2 ) );
750 
751  EvtVector4R q_lambda1( E_lambda,
752  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
753  ( ( px_M * ( P_lambda.get( 1 ) ) ) +
754  ( py_M * ( P_lambda.get( 2 ) ) ) ) ),
755  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
756  ( -( py_M * ( P_lambda.get( 1 ) ) ) +
757  ( px_M * ( P_lambda.get( 2 ) ) ) ) ),
758  P_lambda.get( 3 ) );
759 
760  EvtVector4R q_lambda2( E_lambda, q_lambda1.get( 3 ), q_lambda1.get( 1 ),
761  q_lambda1.get( 2 ) );
762 
763  double px = q_lambda2.get( 1 );
764  double py = q_lambda2.get( 2 );
765  double pz = q_lambda2.get( 3 );
766 
767  EvtVector4R q_lambda4(
768  q_lambda2.get( 0 ),
769  ( ( 1 / ( sqrt( pow( q_lambda2.get( 1 ), 2 ) + pow( q_lambda2.get( 2 ), 2 ) +
770  pow( q_lambda2.get( 3 ), 2 ) ) ) ) *
771  ( 1 / ( sqrt( pow( q_lambda2.get( 1 ), 2 ) +
772  pow( q_lambda2.get( 2 ), 2 ) ) ) ) *
773  ( ( q_lambda2.get( 1 ) ) * ( q_lambda2.get( 1 ) ) *
774  ( q_lambda2.get( 3 ) ) +
775  ( ( q_lambda2.get( 2 ) ) * ( q_lambda2.get( 2 ) ) *
776  ( q_lambda2.get( 3 ) ) ) -
777  ( ( q_lambda2.get( 3 ) ) * ( pow( q_lambda2.get( 1 ), 2 ) +
778  pow( q_lambda2.get( 2 ), 2 ) ) ) ) ),
779  ( ( ( ( q_lambda2.get( 2 ) ) * ( q_lambda2.get( 1 ) ) ) -
780  ( ( q_lambda2.get( 1 ) ) * ( q_lambda2.get( 2 ) ) ) ) /
781  ( sqrt( pow( q_lambda2.get( 1 ), 2 ) + pow( q_lambda2.get( 2 ), 2 ) ) ) ),
782  ( ( ( 1 / sqrt( pow( q_lambda2.get( 1 ), 2 ) + pow( q_lambda2.get( 2 ), 2 ) +
783  pow( q_lambda2.get( 3 ), 2 ) ) ) *
784  ( ( ( q_lambda2.get( 1 ) ) * ( q_lambda2.get( 1 ) ) ) +
785  ( ( q_lambda2.get( 2 ) ) * ( q_lambda2.get( 2 ) ) ) +
786  ( ( q_lambda2.get( 3 ) ) * ( q_lambda2.get( 3 ) ) ) ) ) ) );
787 
788  EvtVector4R q_proton1( E_proton, P * sin( theta ) * cos( phi ),
789  P * sin( theta ) * sin( phi ), P * cos( theta ) );
790  EvtVector4R q_pion1( E_pion, -P * sin( theta ) * cos( phi ),
791  -P * sin( theta ) * sin( phi ), -P * cos( theta ) );
792 
793  EvtVector4R q_proton3(
794  E_proton,
795  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
796  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
797  ( ( q_proton1.get( 1 ) ) * ( px ) * ( pz ) -
798  ( ( q_proton1.get( 2 ) ) * ( py ) *
799  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) +
800  ( ( ( q_proton1.get( 3 ) ) ) *
801  ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) * ( px ) ) ) ),
802  ( ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) *
803  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
804  ( ( ( q_proton1.get( 1 ) ) ) * ( py ) * ( pz ) +
805  ( ( q_proton1.get( 2 ) ) * ( px ) *
806  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) +
807  ( ( ( q_proton1.get( 3 ) ) ) *
808  ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) * ( py ) ) ) ),
809  ( ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) *
810  ( ( -( q_proton1.get( 1 ) ) ) *
811  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) +
812  ( ( q_proton1.get( 3 ) ) * ( pz ) ) ) ) );
813 
814  EvtVector4R q_pion3(
815  E_pion,
816  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
817  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
818  ( ( q_pion1.get( 1 ) ) * ( px ) * ( pz ) -
819  ( ( q_pion1.get( 2 ) ) * ( py ) *
820  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) +
821  ( ( ( q_pion1.get( 3 ) ) ) *
822  ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) * ( px ) ) ) ),
823  ( ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) *
824  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
825  ( ( q_pion1.get( 1 ) ) * ( py ) * ( pz ) +
826  ( ( q_pion1.get( 2 ) ) * ( px ) *
827  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) +
828  ( ( ( q_pion1.get( 3 ) ) ) *
829  ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) * ( py ) ) ) ),
830  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
831  ( ( -( q_pion1.get( 1 ) ) ) * ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) +
832  ( ( q_pion1.get( 3 ) ) * ( pz ) ) ) ) );
833 
834  EvtVector4R q_proton5( q_proton3.get( 0 ), ( q_proton3.get( 2 ) ),
835  ( q_proton3.get( 3 ) ), ( q_proton3.get( 1 ) ) );
836 
837  EvtVector4R q_pion5( q_pion3.get( 0 ), ( q_pion3.get( 2 ) ),
838  ( q_pion3.get( 3 ) ), ( q_pion3.get( 1 ) ) );
839 
840  EvtVector4R q_proton( q_proton5.get( 0 ),
841  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
842  ( ( px_M * ( q_proton5.get( 1 ) ) ) -
843  ( py_M * ( q_proton5.get( 2 ) ) ) ) ),
844  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
845  ( ( py_M * ( q_proton5.get( 1 ) ) ) +
846  ( px_M * ( q_proton5.get( 2 ) ) ) ) ),
847  ( q_proton5.get( 3 ) ) );
848 
849  EvtVector4R q_pion(
850  q_pion5.get( 0 ),
851  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
852  ( ( px_M * ( q_pion5.get( 1 ) ) ) - ( py_M * ( q_pion5.get( 2 ) ) ) ) ),
853  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
854  ( ( py_M * ( q_pion5.get( 1 ) ) ) + ( px_M * ( q_pion5.get( 2 ) ) ) ) ),
855  ( q_pion5.get( 3 ) ) );
856 
857  EvtGenReport( EVTGEN_INFO, fname.c_str() )
858  << " Lambdab px: " << px_M << std::endl;
859  EvtGenReport( EVTGEN_INFO, fname.c_str() )
860  << " Lambdab py: " << py_M << std::endl;
861  EvtGenReport( EVTGEN_INFO, fname.c_str() )
862  << " Lambdab pz: " << pz_M << std::endl;
863  EvtGenReport( EVTGEN_INFO, fname.c_str() )
864  << " Lambdab E: " << E_M << std::endl;
865  EvtGenReport( EVTGEN_INFO, fname.c_str() )
866  << " Lambdab2 px: " << q_lambdab2.get( 1 ) << std::endl;
867  EvtGenReport( EVTGEN_INFO, fname.c_str() )
868  << " Lambdab2 py: " << q_lambdab2.get( 2 ) << std::endl;
869  EvtGenReport( EVTGEN_INFO, fname.c_str() )
870  << " Lambdab2 pz: " << q_lambdab2.get( 3 ) << std::endl;
871  EvtGenReport( EVTGEN_INFO, fname.c_str() )
872  << " Lambdab2 E: " << q_lambdab2.get( 0 ) << std::endl;
873  EvtGenReport( EVTGEN_INFO, fname.c_str() )
874  << " Lambdab3 px: " << q_lambdab3.get( 1 ) << std::endl;
875  EvtGenReport( EVTGEN_INFO, fname.c_str() )
876  << " Lambdab3 py: " << q_lambdab3.get( 2 ) << std::endl;
877  EvtGenReport( EVTGEN_INFO, fname.c_str() )
878  << " Lambdab3 pz: " << q_lambdab3.get( 3 ) << std::endl;
879  EvtGenReport( EVTGEN_INFO, fname.c_str() )
880  << " Lambdab3 E: " << q_lambdab3.get( 0 ) << std::endl;
881  EvtGenReport( EVTGEN_INFO, fname.c_str() )
882  << " Lambda 0 px: " << P_lambda.get( 1 ) << std::endl;
883  EvtGenReport( EVTGEN_INFO, fname.c_str() )
884  << " Lambda 0 py: " << P_lambda.get( 2 ) << std::endl;
885  EvtGenReport( EVTGEN_INFO, fname.c_str() )
886  << " Lambda 0 pz: " << P_lambda.get( 3 ) << std::endl;
887  EvtGenReport( EVTGEN_INFO, fname.c_str() )
888  << " Lambda 0 E: " << P_lambda.get( 0 ) << std::endl;
889  EvtGenReport( EVTGEN_INFO, fname.c_str() )
890  << " Lambda 1 px: " << q_lambda1.get( 1 ) << std::endl;
891  EvtGenReport( EVTGEN_INFO, fname.c_str() )
892  << " Lambda 1 py: " << q_lambda1.get( 2 ) << std::endl;
893  EvtGenReport( EVTGEN_INFO, fname.c_str() )
894  << " Lambda 1 pz: " << q_lambda1.get( 3 ) << std::endl;
895  EvtGenReport( EVTGEN_INFO, fname.c_str() )
896  << " Lambda 1 E: " << q_lambda1.get( 0 ) << std::endl;
897  EvtGenReport( EVTGEN_INFO, fname.c_str() )
898  << " Lambda 2 px: " << q_lambda2.get( 1 ) << std::endl;
899  EvtGenReport( EVTGEN_INFO, fname.c_str() )
900  << " Lambda 2 py: " << q_lambda2.get( 2 ) << std::endl;
901  EvtGenReport( EVTGEN_INFO, fname.c_str() )
902  << " Lambda 2 pz: " << q_lambda2.get( 3 ) << std::endl;
903  EvtGenReport( EVTGEN_INFO, fname.c_str() )
904  << " Lambda 2 E: " << q_lambda2.get( 0 ) << std::endl;
905 
906  EvtGenReport( EVTGEN_INFO, fname.c_str() )
907  << " Lambda px: " << px << std::endl;
908  EvtGenReport( EVTGEN_INFO, fname.c_str() )
909  << " Lambda py: " << py << std::endl;
910  EvtGenReport( EVTGEN_INFO, fname.c_str() )
911  << " Lambda pz: " << pz << std::endl;
912 
913  EvtGenReport( EVTGEN_INFO, fname.c_str() )
914  << " pion 1 px: " << q_pion1.get( 1 ) << std::endl;
915  EvtGenReport( EVTGEN_INFO, fname.c_str() )
916  << " pion 1 py: " << q_pion1.get( 2 ) << std::endl;
917  EvtGenReport( EVTGEN_INFO, fname.c_str() )
918  << " pion 1 pz: " << q_pion1.get( 3 ) << std::endl;
919  EvtGenReport( EVTGEN_INFO, fname.c_str() )
920  << " pion 1 E: " << q_pion1.get( 0 ) << std::endl;
921 
922  EvtGenReport( EVTGEN_INFO, fname.c_str() )
923  << " pion 3 px: " << q_pion3.get( 1 ) << std::endl;
924  EvtGenReport( EVTGEN_INFO, fname.c_str() )
925  << " pion 3 px: " << q_pion3.get( 1 ) << std::endl;
926  EvtGenReport( EVTGEN_INFO, fname.c_str() )
927  << " pion 3 py: " << q_pion3.get( 2 ) << std::endl;
928  EvtGenReport( EVTGEN_INFO, fname.c_str() )
929  << " pion 3 pz: " << q_pion3.get( 3 ) << std::endl;
930  EvtGenReport( EVTGEN_INFO, fname.c_str() )
931  << " pion 3 E: " << q_pion3.get( 0 ) << std::endl;
932 
933  EvtGenReport( EVTGEN_INFO, fname.c_str() )
934  << " pion 5 px: " << q_pion5.get( 1 ) << std::endl;
935  EvtGenReport( EVTGEN_INFO, fname.c_str() )
936  << " pion 5 py: " << q_pion5.get( 2 ) << std::endl;
937  EvtGenReport( EVTGEN_INFO, fname.c_str() )
938  << " pion 5 pz: " << q_pion5.get( 3 ) << std::endl;
939  EvtGenReport( EVTGEN_INFO, fname.c_str() )
940  << " pion 5 E: " << q_pion5.get( 0 ) << std::endl;
941 
942  EvtGenReport( EVTGEN_INFO, fname.c_str() )
943  << " proton 1 px: " << q_proton1.get( 1 ) << std::endl;
944  EvtGenReport( EVTGEN_INFO, fname.c_str() )
945  << " proton 1 py: " << q_proton1.get( 2 ) << std::endl;
946  EvtGenReport( EVTGEN_INFO, fname.c_str() )
947  << " proton 1 pz: " << q_proton1.get( 3 ) << std::endl;
948  EvtGenReport( EVTGEN_INFO, fname.c_str() )
949  << " proton 1 E: " << q_proton1.get( 0 ) << std::endl;
950 
951  EvtGenReport( EVTGEN_INFO, fname.c_str() )
952  << " proton 3 px: " << q_proton3.get( 1 ) << std::endl;
953  EvtGenReport( EVTGEN_INFO, fname.c_str() )
954  << " proton 3 py: " << q_proton3.get( 2 ) << std::endl;
955  EvtGenReport( EVTGEN_INFO, fname.c_str() )
956  << " proton 3 pz: " << q_proton3.get( 3 ) << std::endl;
957  EvtGenReport( EVTGEN_INFO, fname.c_str() )
958  << " proton 3 E: " << q_proton3.get( 0 ) << std::endl;
959 
960  EvtGenReport( EVTGEN_INFO, fname.c_str() )
961  << " proton 5 px: " << q_proton5.get( 1 ) << std::endl;
962  EvtGenReport( EVTGEN_INFO, fname.c_str() )
963  << " proton 5 py: " << q_proton5.get( 2 ) << std::endl;
964  EvtGenReport( EVTGEN_INFO, fname.c_str() )
965  << " proton 5 pz: " << q_proton5.get( 3 ) << std::endl;
966  EvtGenReport( EVTGEN_INFO, fname.c_str() )
967  << " proton 5 E: " << q_proton5.get( 0 ) << std::endl;
968 
969  EvtGenReport( EVTGEN_INFO, fname.c_str() )
970  << " proton px: " << q_proton.get( 1 ) << std::endl;
971  EvtGenReport( EVTGEN_INFO, fname.c_str() )
972  << " proton py: " << q_proton.get( 2 ) << std::endl;
973  EvtGenReport( EVTGEN_INFO, fname.c_str() )
974  << "proton pz: " << q_proton.get( 3 ) << std::endl;
975  EvtGenReport( EVTGEN_INFO, fname.c_str() )
976  << " pion px: " << q_pion.get( 1 ) << std::endl;
977  EvtGenReport( EVTGEN_INFO, fname.c_str() )
978  << " pion py: " << q_pion.get( 2 ) << std::endl;
979  EvtGenReport( EVTGEN_INFO, fname.c_str() )
980  << " pion pz: " << q_pion.get( 3 ) << std::endl;
981 
982  ;
983 
985 
986  //set quadrivectors to particles
987  proton->init( getDaugs()[0], q_proton );
988  pion->init( getDaugs()[1], q_pion );
989 
990  //computate pdf
991  //double pdf = 1 + A*B*cos(theta) - EvtConst::pi/2.0*C*A*real(D*EvtComplex(cos(phi),sin(phi)))*sin(theta);
992  double pdf = 1 + A * B * cos( theta ) +
993  2 * A * real( D * EvtComplex( cos( phi ), sin( phi ) ) ) *
994  sin( theta );
995  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
996  << " Lambda decay pdf value : " << pdf << std::endl;
997  //set probability
998  setProb( pdf );
999 
1000  return;
1001 }
1002 
1003 //************************************************************************
1004 //* *
1005 //* Class EvtLambda2PPiForLambdaB2LambdaV *
1006 //* *
1007 //************************************************************************
1008 
1009 //------------------------------------------------------------------------
1010 // Constructor
1011 //------------------------------------------------------------------------
1012 EvtV2VpVmForLambdaB2LambdaV::EvtV2VpVmForLambdaB2LambdaV()
1013 {
1014  //set facility name
1015  fname = "EvtGen.EvtV2VpVmForLambdaB2LambdaV";
1016 }
1017 
1018 //------------------------------------------------------------------------
1019 // Method 'getName'
1020 //------------------------------------------------------------------------
1021 std::string EvtV2VpVmForLambdaB2LambdaV::getName()
1022 {
1023  return "V2VPVMFORLAMBDAB2LAMBDAV";
1024 }
1025 
1026 //------------------------------------------------------------------------
1027 // Method 'clone'
1028 //------------------------------------------------------------------------
1029 EvtDecayBase* EvtV2VpVmForLambdaB2LambdaV::clone()
1030 {
1031  return new EvtV2VpVmForLambdaB2LambdaV;
1032 }
1033 
1034 //------------------------------------------------------------------------
1035 // Method 'initProbMax'
1036 //------------------------------------------------------------------------
1037 void EvtV2VpVmForLambdaB2LambdaV::initProbMax()
1038 {
1039  //maximum
1040  double Max = 0;
1041  if ( Vtype == VID::JPSI ) {
1042  if ( ( 1 - 3 * A ) > 0 )
1043  Max = 2 * ( 1 - A );
1044  else
1045  Max = 1 + A;
1046  } else {
1047  if ( ( 3 * A - 1 ) >= 0 )
1048  Max = 2 * A;
1049  else
1050  Max = 1 - A;
1051  }
1052 
1053  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1054  << " PDF max value : " << Max << std::endl;
1055  setProbMax( Max );
1056 }
1057 
1058 //------------------------------------------------------------------------
1059 // Method 'init'
1060 //------------------------------------------------------------------------
1061 void EvtV2VpVmForLambdaB2LambdaV::init()
1062 {
1063  //introduction
1064  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1065  << " ***********************************************************"
1066  << std::endl;
1067  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1068  << " * Event Model Class : EvtV2VpVmForLambdaB2LambdaV *"
1069  << std::endl;
1070  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1071  << " ***********************************************************"
1072  << std::endl;
1073 
1074  //check the number of arguments
1075  checkNArg( 2 );
1076 
1077  //check the number of daughters
1078  checkNDaug( 2 );
1079 
1080  const EvtId Id_mother = getParentId();
1081  const EvtId Id_daug1 = getDaug( 0 );
1082  const EvtId Id_daug2 = getDaug( 1 );
1083 
1084  //identification meson V
1085  if ( getArg( 1 ) == 1 )
1086  Vtype = VID::JPSI;
1087  else if ( getArg( 1 ) == 2 )
1088  Vtype = VID::RHO;
1089  else if ( getArg( 1 ) == 3 )
1090  Vtype = VID::OMEGA;
1091  else if ( getArg( 1 ) == 4 )
1092  Vtype = VID::RHO_OMEGA_MIXING;
1093  else {
1094  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1095  << " Vtype " << getArg( 1 ) << " is unknown" << std::endl;
1096  abort();
1097  }
1098 
1099  //vector meson V
1100  if ( Id_mother == EvtPDL::getId( "J/psi" ) && Vtype == VID::JPSI ) {
1101  } else if ( Id_mother == EvtPDL::getId( "omega" ) && Vtype == VID::OMEGA ) {
1102  } else if ( Id_mother == EvtPDL::getId( "rho0" ) && Vtype == VID::RHO ) {
1103  } else if ( ( Id_mother == EvtPDL::getId( "rho0" ) ||
1104  Id_mother == EvtPDL::getId( "omega" ) ) &&
1105  Vtype == VID::RHO_OMEGA_MIXING ) {
1106  } else {
1107  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1108  << " Mother is not a J/psi, phi or rho0 but a "
1109  << EvtPDL::name( Id_mother ) << std::endl;
1110  abort();
1111  }
1112 
1113  //daughters for each V possibility
1114  switch ( Vtype ) {
1115  case VID::JPSI:
1116  if ( Id_daug1 != EvtPDL::getId( "mu+" ) ) {
1117  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1118  << " Daughter1 is not a mu+, but a "
1119  << EvtPDL::name( Id_daug1 ) << std::endl;
1120  abort();
1121  }
1122  if ( Id_daug2 != EvtPDL::getId( "mu-" ) ) {
1123  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1124  << " Daughter2 is not a mu-, but a "
1125  << EvtPDL::name( Id_daug2 ) << std::endl;
1126  abort();
1127  }
1128  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1129  << " Decay mode successfully initialized : J/psi -> mu+ mu-"
1130  << std::endl;
1131  break;
1132 
1133  case VID::RHO:
1134  case VID::OMEGA:
1135  case VID::RHO_OMEGA_MIXING:
1136  if ( Id_daug1 != EvtPDL::getId( "pi+" ) ) {
1137  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1138  << " Daughter1 is not a pi+, but a "
1139  << EvtPDL::name( Id_daug1 ) << std::endl;
1140  abort();
1141  }
1142  if ( Id_daug2 != EvtPDL::getId( "pi-" ) ) {
1143  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1144  << " Daughter2 is not a pi-, but a "
1145  << EvtPDL::name( Id_daug2 ) << std::endl;
1146  abort();
1147  }
1148  if ( Vtype == VID::RHO )
1149  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1150  << " Decay mode successfully initialized : rho0 -> pi+ pi-"
1151  << std::endl;
1152  if ( Vtype == VID::OMEGA )
1153  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1154  << " Decay mode successfully initialized : omega -> pi+ pi-"
1155  << std::endl;
1156  if ( Vtype == VID::RHO_OMEGA_MIXING )
1157  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1158  << " Decay mode successfully initialized : rho-omega mixing -> pi+ pi-"
1159  << std::endl;
1160  break;
1161 
1162  default:
1163  EvtGenReport( EVTGEN_ERROR, fname.c_str() )
1164  << "No decay mode chosen ! " << std::endl;
1165  abort();
1166  break;
1167  }
1168 
1169  //fix dynamics parameters
1170  switch ( Vtype ) {
1171  case VID::JPSI:
1172  A = 0.66;
1173  break;
1174  case VID::RHO:
1175  case VID::OMEGA:
1176  case VID::RHO_OMEGA_MIXING:
1177  A = 0.79;
1178  break;
1179  default:
1180  A = 0;
1181  break;
1182  }
1183 
1184  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1185  << " V decay parameters : " << std::endl;
1186  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1187  << " - V density matrix rho00 A = " << A << std::endl;
1188 }
1189 
1190 //------------------------------------------------------------------------
1191 // Method 'decay'
1192 //------------------------------------------------------------------------
1193 void EvtV2VpVmForLambdaB2LambdaV::decay( EvtParticle* V )
1194 {
1195  V->makeDaughters( getNDaug(), getDaugs() );
1196 
1197  //get Vp and Vm particles
1198  EvtParticle* Vp = V->getDaug( 0 );
1199  EvtParticle* Vm = V->getDaug( 1 );
1200 
1201  //get resonance masses
1202  // - V -> mass given by EvtLambdaB2LambdaV class
1203  // - Vp & Vm -> nominal mass
1204  double MASS_V = V->mass();
1205  double MASS_VM = 0;
1206  switch ( Vtype ) {
1207  case VID::JPSI:
1208  MASS_VM = EvtPDL::getMeanMass( EvtPDL::getId( "mu-" ) );
1209  break;
1210  case VID::RHO:
1211  case VID::OMEGA:
1212  case VID::RHO_OMEGA_MIXING:
1213  MASS_VM = EvtPDL::getMeanMass( EvtPDL::getId( "pi-" ) );
1214  break;
1215  default:
1216  MASS_VM = 0;
1217  break;
1218  }
1219  double MASS_VP = MASS_VM;
1220 
1221  //generate random angles
1222  double phi = EvtRandom::Flat( 0, 2 * EvtConst::pi );
1223  double theta = acos( EvtRandom::Flat( -1, +1 ) );
1224  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1225  << " Angular angles : theta = " << theta << " ; phi = " << phi
1226  << std::endl;
1227 
1228  //computate resonance quadrivectors
1229  double E_Vp = ( MASS_V * MASS_V + MASS_VP * MASS_VP - MASS_VM * MASS_VM ) /
1230  ( 2 * MASS_V );
1231  double E_Vm = ( MASS_V * MASS_V + MASS_VM * MASS_VM - MASS_VP * MASS_VP ) /
1232  ( 2 * MASS_V );
1233  double P = sqrt( E_Vp * E_Vp - Vp->mass() * Vp->mass() );
1234 
1235  EvtVector4R P_V = V->getP4();
1236  EvtParticle* Mother_V = V->getParent();
1237  EvtVector4R lambdab = Mother_V->getP4();
1238 
1239  double E_V = ( P_V.get( 0 ) );
1240  double px_M = lambdab.get( 1 );
1241  double py_M = lambdab.get( 2 );
1242  double pz_M = lambdab.get( 3 );
1243  double E_M = lambdab.get( 0 );
1244 
1245  EvtVector4R q_lambdab2(
1246  E_M,
1247  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1248  ( ( px_M * ( px_M ) ) + ( py_M * ( py_M ) ) ) ),
1249  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1250  ( -( ( py_M ) * ( px_M ) ) + ( px_M * ( py_M ) ) ) ),
1251  ( pz_M ) );
1252 
1253  EvtVector4R q_lambdab3( E_M, q_lambdab2.get( 3 ), q_lambdab2.get( 1 ),
1254  q_lambdab2.get( 2 ) );
1255 
1256  EvtVector4R q_V1(
1257  E_V,
1258  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1259  ( ( px_M * ( P_V.get( 1 ) ) ) + ( py_M * ( P_V.get( 2 ) ) ) ) ),
1260  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1261  ( -( py_M * ( P_V.get( 1 ) ) ) + ( px_M * ( P_V.get( 2 ) ) ) ) ),
1262  P_V.get( 3 ) );
1263 
1264  EvtVector4R q_V2( E_V, q_V1.get( 3 ), q_V1.get( 1 ), q_V1.get( 2 ) );
1265 
1266  double px = -( q_V2.get( 1 ) );
1267  double py = -( q_V2.get( 2 ) );
1268  double pz = -( q_V2.get( 3 ) );
1269 
1270  EvtVector4R q_V4(
1271  q_V2.get( 0 ),
1272  ( ( 1 / ( sqrt( pow( q_V2.get( 1 ), 2 ) + pow( q_V2.get( 2 ), 2 ) +
1273  pow( q_V2.get( 3 ), 2 ) ) ) ) *
1274  ( 1 / ( sqrt( pow( q_V2.get( 1 ), 2 ) + pow( q_V2.get( 2 ), 2 ) ) ) ) *
1275  ( ( q_V2.get( 1 ) ) * ( q_V2.get( 1 ) ) * ( q_V2.get( 3 ) ) +
1276  ( ( q_V2.get( 2 ) ) * ( q_V2.get( 2 ) ) * ( q_V2.get( 3 ) ) ) -
1277  ( ( q_V2.get( 3 ) ) *
1278  ( pow( q_V2.get( 1 ), 2 ) + pow( q_V2.get( 2 ), 2 ) ) ) ) ),
1279  ( ( ( ( q_V2.get( 2 ) ) * ( q_V2.get( 1 ) ) ) -
1280  ( ( q_V2.get( 1 ) ) * ( q_V2.get( 2 ) ) ) ) /
1281  ( sqrt( pow( q_V2.get( 1 ), 2 ) + pow( q_V2.get( 2 ), 2 ) ) ) ),
1282  ( ( ( 1 / sqrt( pow( q_V2.get( 1 ), 2 ) + pow( q_V2.get( 2 ), 2 ) +
1283  pow( q_V2.get( 3 ), 2 ) ) ) *
1284  ( ( ( q_V2.get( 1 ) ) * ( q_V2.get( 1 ) ) ) +
1285  ( ( q_V2.get( 2 ) ) * ( q_V2.get( 2 ) ) ) +
1286  ( ( q_V2.get( 3 ) ) * ( q_V2.get( 3 ) ) ) ) ) ) );
1287 
1288  EvtVector4R q_Vp1( E_Vp, P * sin( theta ) * cos( phi ),
1289  P * sin( theta ) * sin( phi ), P * cos( theta ) );
1290  EvtVector4R q_Vm1( E_Vm, -P * sin( theta ) * cos( phi ),
1291  -P * sin( theta ) * sin( phi ), -P * cos( theta ) );
1292 
1293  EvtVector4R q_Vp3(
1294  q_Vp1.get( 0 ),
1295  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
1296  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
1297  ( ( q_Vp1.get( 1 ) ) * ( px ) * ( pz ) +
1298  ( ( q_Vp1.get( 2 ) ) * ( py ) *
1299  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) -
1300  ( ( ( q_Vp1.get( 3 ) ) ) * ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) *
1301  ( px ) ) ) ),
1302  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
1303  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
1304  ( ( ( q_Vp1.get( 1 ) ) ) * ( py ) * ( pz ) -
1305  ( ( q_Vp1.get( 2 ) ) * ( px ) *
1306  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) -
1307  ( ( ( q_Vp1.get( 3 ) ) ) * ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) *
1308  ( py ) ) ) ),
1309  ( ( -( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) *
1310  ( ( q_Vp1.get( 1 ) ) * ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) +
1311  ( ( q_Vp1.get( 3 ) ) * ( pz ) ) ) ) );
1312 
1313  EvtVector4R q_Vm3(
1314  q_Vm1.get( 0 ),
1315  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
1316  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
1317  ( ( q_Vm1.get( 1 ) ) * ( px ) * ( pz ) +
1318  ( ( q_Vm1.get( 2 ) ) * ( py ) *
1319  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) -
1320  ( ( ( q_Vm1.get( 3 ) ) ) * ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) *
1321  ( px ) ) ) ),
1322  ( ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) *
1323  ( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) *
1324  ( ( ( q_Vm1.get( 1 ) ) ) * ( py ) * ( pz ) -
1325  ( ( q_Vm1.get( 2 ) ) * ( px ) *
1326  ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) -
1327  ( ( ( q_Vm1.get( 3 ) ) ) * ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) *
1328  ( py ) ) ) ),
1329  ( ( -( 1 / ( sqrt( pow( px, 2 ) + pow( py, 2 ) + pow( pz, 2 ) ) ) ) ) *
1330  ( ( q_Vm1.get( 1 ) ) * ( ( sqrt( pow( px, 2 ) + pow( py, 2 ) ) ) ) +
1331  ( ( q_Vm1.get( 3 ) ) * ( pz ) ) ) ) );
1332 
1333  EvtVector4R q_Vp5( q_Vp3.get( 0 ), ( q_Vp3.get( 2 ) ), ( q_Vp3.get( 3 ) ),
1334  ( q_Vp3.get( 1 ) ) );
1335 
1336  EvtVector4R q_Vm5( q_Vm3.get( 0 ), ( q_Vm3.get( 2 ) ), ( q_Vm3.get( 3 ) ),
1337  ( q_Vm3.get( 1 ) ) );
1338 
1339  EvtVector4R q_Vp(
1340  q_Vp5.get( 0 ),
1341  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1342  ( ( px_M * ( q_Vp5.get( 1 ) ) ) - ( py_M * ( q_Vp5.get( 2 ) ) ) ) ),
1343  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1344  ( ( py_M * ( q_Vp5.get( 1 ) ) ) + ( px_M * ( q_Vp5.get( 2 ) ) ) ) ),
1345  ( q_Vp5.get( 3 ) ) );
1346 
1347  EvtVector4R q_Vm(
1348  q_Vm5.get( 0 ),
1349  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1350  ( ( px_M * ( q_Vm5.get( 1 ) ) ) - ( py_M * ( q_Vm5.get( 2 ) ) ) ) ),
1351  ( ( 1 / ( sqrt( pow( px_M, 2 ) + pow( py_M, 2 ) ) ) ) *
1352  ( ( py_M * ( q_Vm5.get( 1 ) ) ) + ( px_M * ( q_Vm5.get( 2 ) ) ) ) ),
1353  ( q_Vm5.get( 3 ) ) );
1354 
1355  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1356  << " Lambdab px: " << px_M << std::endl;
1357  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1358  << " Lambdab py: " << py_M << std::endl;
1359  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1360  << " Lambdab pz: " << pz_M << std::endl;
1361  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1362  << " Lambdab E: " << E_M << std::endl;
1363  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1364  << " Lambdab2 px: " << q_lambdab2.get( 1 ) << std::endl;
1365  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1366  << " Lambdab2 py: " << q_lambdab2.get( 2 ) << std::endl;
1367  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1368  << " Lambdab2 pz: " << q_lambdab2.get( 3 ) << std::endl;
1369  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1370  << " Lambdab2 E: " << q_lambdab2.get( 0 ) << std::endl;
1371  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1372  << " Lambdab3 px: " << q_lambdab3.get( 1 ) << std::endl;
1373  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1374  << " Lambdab3 py: " << q_lambdab3.get( 2 ) << std::endl;
1375  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1376  << " Lambdab3 pz: " << q_lambdab3.get( 3 ) << std::endl;
1377  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1378  << " Lambdab3 E: " << q_lambdab3.get( 0 ) << std::endl;
1379  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1380  << " V 0 px: " << P_V.get( 1 ) << std::endl;
1381  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1382  << " V 0 py: " << P_V.get( 2 ) << std::endl;
1383  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1384  << " V 0 pz: " << P_V.get( 3 ) << std::endl;
1385  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1386  << " V 0 E: " << P_V.get( 0 ) << std::endl;
1387  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1388  << " V 1 px: " << q_V1.get( 1 ) << std::endl;
1389  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1390  << " V 1 py: " << q_V1.get( 2 ) << std::endl;
1391  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1392  << " V 1 pz: " << q_V1.get( 3 ) << std::endl;
1393  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1394  << " V 1 E: " << q_V1.get( 0 ) << std::endl;
1395  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1396  << " V 2 px: " << q_V2.get( 1 ) << std::endl;
1397  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1398  << " V 2 py: " << q_V2.get( 2 ) << std::endl;
1399  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1400  << " V 2 pz: " << q_V2.get( 3 ) << std::endl;
1401  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1402  << " V 2 E: " << q_V2.get( 0 ) << std::endl;
1403  EvtGenReport( EVTGEN_INFO, fname.c_str() ) << " V px: " << px << std::endl;
1404  EvtGenReport( EVTGEN_INFO, fname.c_str() ) << " V py: " << py << std::endl;
1405  EvtGenReport( EVTGEN_INFO, fname.c_str() ) << " V pz: " << pz << std::endl;
1406  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1407  << " Vm 1 px: " << q_Vm1.get( 1 ) << std::endl;
1408  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1409  << " Vm 1 py: " << q_Vm1.get( 2 ) << std::endl;
1410  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1411  << " Vm 1 pz: " << q_Vm1.get( 3 ) << std::endl;
1412  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1413  << " Vm 1 E: " << q_Vm1.get( 0 ) << std::endl;
1414  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1415  << " Vm 3 px: " << q_Vm3.get( 1 ) << std::endl;
1416  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1417  << " Vm 3 px: " << q_Vm3.get( 1 ) << std::endl;
1418  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1419  << " Vm 3 py: " << q_Vm3.get( 2 ) << std::endl;
1420  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1421  << " Vm 3 pz: " << q_Vm3.get( 3 ) << std::endl;
1422  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1423  << " Vm 3 E: " << q_Vm3.get( 0 ) << std::endl;
1424  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1425  << " Vm 5 px: " << q_Vm5.get( 1 ) << std::endl;
1426  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1427  << " Vm 5 py: " << q_Vm5.get( 2 ) << std::endl;
1428  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1429  << " Vm 5 pz: " << q_Vm5.get( 3 ) << std::endl;
1430  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1431  << " Vm 5 E: " << q_Vm5.get( 0 ) << std::endl;
1432  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1433  << " Vp 1 px: " << q_Vp1.get( 1 ) << std::endl;
1434  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1435  << " Vp 1 py: " << q_Vp1.get( 2 ) << std::endl;
1436  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1437  << " Vp 1 pz: " << q_Vp1.get( 3 ) << std::endl;
1438  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1439  << " Vp 1 E: " << q_Vp1.get( 0 ) << std::endl;
1440  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1441  << " Vp 3 px: " << q_Vp3.get( 1 ) << std::endl;
1442  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1443  << " Vp 3 py: " << q_Vp3.get( 2 ) << std::endl;
1444  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1445  << " Vp 3 pz: " << q_Vp3.get( 3 ) << std::endl;
1446  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1447  << " Vp 3 E: " << q_Vp3.get( 0 ) << std::endl;
1448  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1449  << " Vp 5 px: " << q_Vp5.get( 1 ) << std::endl;
1450  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1451  << " Vp 5 py: " << q_Vp5.get( 2 ) << std::endl;
1452  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1453  << " Vp 5 pz: " << q_Vp5.get( 3 ) << std::endl;
1454  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1455  << " Vp 5 E: " << q_Vp5.get( 0 ) << std::endl;
1456  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1457  << " Vp px: " << q_Vp.get( 1 ) << std::endl;
1458  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1459  << " Vp py: " << q_Vp.get( 2 ) << std::endl;
1460  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1461  << "Vp pz: " << q_Vp.get( 3 ) << std::endl;
1462  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1463  << " Vm px: " << q_Vm.get( 1 ) << std::endl;
1464  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1465  << " Vm py: " << q_Vm.get( 2 ) << std::endl;
1466  EvtGenReport( EVTGEN_INFO, fname.c_str() )
1467  << " Vm pz: " << q_Vm.get( 3 ) << std::endl;
1468 
1469  //set quadrivectors to particles
1470  Vp->init( getDaugs()[0], q_Vp );
1471  Vm->init( getDaugs()[1], q_Vm );
1472 
1473  //computate pdf
1474  double pdf = 0;
1475  if ( Vtype == VID::JPSI ) {
1476  //leptonic case
1477  pdf = ( 1 - 3 * A ) * cos( theta ) * cos( theta ) + ( 1 + A );
1478  } else {
1479  //hadronic case
1480  pdf = ( 3 * A - 1 ) * cos( theta ) * cos( theta ) + ( 1 - A );
1481  }
1482  EvtGenReport( EVTGEN_DEBUG, fname.c_str() )
1483  << " V decay pdf value : " << pdf << std::endl;
1484 
1485  //set probability
1486  setProb( pdf );
1487  return;
1488 }
EvtDecayProb::setProb
void setProb(double prob)
Definition: EvtDecayProb.hh:32
EvtV2VpVmForLambdaB2LambdaV::decay
void decay(EvtParticle *V) override
Definition: EvtLambdaB2LambdaV.cpp:1193
EvtPDL::name
static std::string name(EvtId i)
Definition: EvtPDL.cpp:382
EvtParticle::getParent
EvtParticle * getParent() const
Definition: EvtParticle.cpp:96
EvtDecayBase::getArg
double getArg(unsigned int j)
Definition: EvtDecayBase.cpp:604
EvtGenReport
std::ostream & EvtGenReport(EvtGenSeverity severity, const char *facility=0)
Definition: EvtReport.cpp:33
EvtParticle::init
virtual void init(EvtId part_n, const EvtVector4R &p4)=0
EvtLambdaB2LambdaV::BreitWignerRelPDF
double BreitWignerRelPDF(double m, double _m0, double _g0)
Definition: EvtLambdaB2LambdaV.cpp:381
EvtPDL::getMeanMass
static double getMeanMass(EvtId i)
Definition: EvtPDL.cpp:314
EvtParticle::makeDaughters
void makeDaughters(unsigned int ndaug, EvtId *id)
Definition: EvtParticle.cpp:1239
EvtDecayBase::getDaugs
EvtId * getDaugs()
Definition: EvtDecayBase.hh:66
EvtLambdaB2LambdaV::RhoOmegaMixingPDF
double RhoOmegaMixingPDF(double m, double _mr, double _gr, double _mo, double _go)
Definition: EvtLambdaB2LambdaV.cpp:391
EvtDecayBase::setProbMax
void setProbMax(double prbmx)
Definition: EvtDecayBase.cpp:321
EvtLambda2PPiForLambdaB2LambdaV::decay
void decay(EvtParticle *lambda) override
Definition: EvtLambdaB2LambdaV.cpp:699
EvtId
Definition: EvtId.hh:27
EvtConst::pi
static const double pi
Definition: EvtConst.hh:26
EvtDecayBase::getParentId
EvtId getParentId() const
Definition: EvtDecayBase.hh:61
EvtVector4R::get
double get(int i) const
Definition: EvtVector4R.hh:162
EvtDecayBase::checkNDaug
void checkNDaug(int d1, int d2=-1)
Definition: EvtDecayBase.cpp:542
EvtRandom::Flat
static double Flat()
Definition: EvtRandom.cpp:72
EvtLambdaB2LambdaV::getVMass
double getVMass(double MASS_LAMBDAB, double MASS_LAMBDA)
Definition: EvtLambdaB2LambdaV.cpp:421
EvtDecayBase::checkNArg
void checkNArg(int a1, int a2=-1, int a3=-1, int a4=-1)
Definition: EvtDecayBase.cpp:518
EvtLambdaB2LambdaV::getName
std::string getName() override
Definition: EvtLambdaB2LambdaV.cpp:52
EvtLambda2PPiForLambdaB2LambdaV::clone
EvtDecayBase * clone() override
Definition: EvtLambdaB2LambdaV.cpp:511
EvtPDL::getId
static EvtId getId(const std::string &name)
Definition: EvtPDL.cpp:287
EvtParticle::getP4
const EvtVector4R & getP4() const
Definition: EvtParticle.cpp:150
EvtV2VpVmForLambdaB2LambdaV::clone
EvtDecayBase * clone() override
Definition: EvtLambdaB2LambdaV.cpp:1029
EvtLambdaB2LambdaV::Vtype
VID::VectorMesonType Vtype
Definition: EvtLambdaB2LambdaV.hh:79
lambda
double lambda(double q, double m1, double m2)
Definition: EvtFlatQ2.cpp:31
EvtParticle::mass
double mass() const
Definition: EvtParticle.cpp:165
EvtPDL::getWidth
static double getWidth(EvtId i)
Definition: EvtPDL.cpp:352
EvtDecayBase::getNDaug
int getNDaug() const
Definition: EvtDecayBase.hh:65
EvtLambdaB2LambdaV::decay
void decay(EvtParticle *lambdab) override
Definition: EvtLambdaB2LambdaV.cpp:221
EvtV2VpVmForLambdaB2LambdaV::getName
std::string getName() override
Definition: EvtLambdaB2LambdaV.cpp:1021
EvtLambdaB2LambdaV::initProbMax
void initProbMax() override
Definition: EvtLambdaB2LambdaV.cpp:68
EvtParticle::getDaug
EvtParticle * getDaug(int i)
Definition: EvtParticle.cpp:91
EvtPDL::getMass
static double getMass(EvtId i)
Definition: EvtPDL.cpp:319
real
double real(const EvtComplex &c)
Definition: EvtComplex.hh:230
EvtLambdaB2LambdaV::clone
EvtDecayBase * clone() override
Definition: EvtLambdaB2LambdaV.cpp:60
EvtDecayBase::getDaug
EvtId getDaug(int i) const
Definition: EvtDecayBase.hh:67
Generated by   doxygen 1.8.15