TEvtGen/EvtGenModels/EvtLambdaB2LambdaV.cpp

   1 #include "EvtGenModels/EvtLambdaB2LambdaV.hh"
   2 #include "EvtGenBase/EvtRandom.hh"
   3 #include "EvtGenBase/EvtPatches.hh"
   4 #include <stdlib.h>
   5 #include <fstream>
   6 #include <stdio.h>
   7 #include <string>
   8 #include "EvtGenBase/EvtGenKine.hh"
   9 #include "EvtGenBase/EvtParticle.hh"
  10 #include "EvtGenBase/EvtPDL.hh"
  11 #include "EvtGenBase/EvtReport.hh"
  12
  13 using std::fstream ;
  14 //************************************************************************
  15 //*                                                                      *
  16 //*                      Class EvtLambdaB2LambdaV                        *
  17 //*                                                                      *
  18 //************************************************************************
  19 //DECLARE_ALGORITHM_FACTORY( EvtLambdaB2LambdaV );
  20
  21 EvtLambdaB2LambdaV::EvtLambdaB2LambdaV()
  22 {
  23   //set facility name
  24   fname="EvtGen.EvtLambdaB2LambdaV";
  25 }
  26
  27
  28 //------------------------------------------------------------------------
  29 // Destructor
  30 //------------------------------------------------------------------------
  31 EvtLambdaB2LambdaV::~EvtLambdaB2LambdaV()
  32 {}
  33
  34
  35 //------------------------------------------------------------------------
  36 // Method 'getName'
  37 //------------------------------------------------------------------------
  38 std::string EvtLambdaB2LambdaV::getName()
  39 {
  40   return  "LAMBDAB2LAMBDAV";
  41 }
  42
  43
  44 //------------------------------------------------------------------------
  45 // Method 'clone'
  46 //------------------------------------------------------------------------
  47 EvtDecayBase* EvtLambdaB2LambdaV::clone()
  48 {
  49   return new EvtLambdaB2LambdaV;
  50
  51 }
  52
  53
  54 //------------------------------------------------------------------------
  55 // Method 'initProbMax'
  56 //------------------------------------------------------------------------
  57 void EvtLambdaB2LambdaV::initProbMax()
  58 {
  59   //maximum (case where C=0)
  60   double Max = 1+fabs(A*B);
  61   report(DEBUG,fname.c_str())<<" PDF max value : "<<Max<<std::endl;
  62   setProbMax(Max);
  63 }
  64
  65
  66 //------------------------------------------------------------------------
  67 // Method 'init'
  68 //------------------------------------------------------------------------
  69 void EvtLambdaB2LambdaV::init()
  70 {
  71   bool antiparticle=false;
  72
  73   //introduction
  74   report(DEBUG,fname.c_str())<< "*************************************************"<<std::endl;
  75   report(DEBUG,fname.c_str())<< "*    Event Model Class : EvtLambdaB2LambdaV     *"<<std::endl;
  76   report(DEBUG,fname.c_str())<< "*************************************************"<<std::endl;
  77
  78   //check the number of arguments
  79   checkNArg(2);
  80
  81
  82   //check the number of daughters
  83   checkNDaug(2);
  84
  85   const EvtId Id_mother = getParentId();
  86   const EvtId Id_daug1  = getDaug(0);
  87   const EvtId Id_daug2  = getDaug(1);
  88
  89
  90   //lambdab identification
  91   if (Id_mother==EvtPDL::getId("Lambda_b0"))
  92   {
  93     antiparticle=false;
  94   }
  95   else if (Id_mother==EvtPDL::getId("anti-Lambda_b0"))
  96   {
  97     antiparticle=true;
  98   }
  99   else
 100   {
 101     report(ERROR,fname.c_str())<<" Mother is not a Lambda_b0 or an anti-Lambda_b0, but a "
 102                           <<EvtPDL::name(Id_mother)<<std::endl;
 103     abort();
 104   }
 105
 106   //lambda
 107   if ( !(Id_daug1==EvtPDL::getId("Lambda0") && !antiparticle ) && !(Id_daug1==EvtPDL::getId("anti-Lambda0") && antiparticle) )
 108   {
 109     if (!antiparticle)
 110     {
 111       report(ERROR,fname.c_str()) << " Daughter1 is not a Lambda0, but a "
 112                                                 << EvtPDL::name(Id_daug1)<<std::endl;
 113     }
 114     else
 115     { report(ERROR,fname.c_str()) << " Daughter1 is not an anti-Lambda0, but a "
 116                                                 << EvtPDL::name(Id_daug1)<<std::endl;
 117     }
 118     abort();
 119   }
 120
 121
 122   //identification meson V
 123   if (getArg(1)==1) Vtype=VID::JPSI;
 124   else if (getArg(1)==2) Vtype=VID::RHO;
 125   else if (getArg(1)==3) Vtype=VID::OMEGA;
 126   else if (getArg(1)==4) Vtype=VID::RHO_OMEGA_MIXING;
 127   else
 128   {
 129     report(ERROR,fname.c_str()) << " Vtype " <<getArg(1)<<" is unknown"<<std::endl;
 130     abort();
 131   }
 132
 133
 134   //check vector meson V
 135   if (Id_daug2==EvtPDL::getId("J/psi") && Vtype==VID::JPSI)
 136   {
 137     if (!antiparticle) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : Lambda_b0 -> Lambda J/psi"<<std::endl;
 138     else report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda J/psi"<<std::endl;
 139   }
 140   else if (Id_daug2==EvtPDL::getId("rho0") && Vtype==VID::RHO )
 141   {
 142     if (!antiparticle) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : Lambda_b0 -> Lambda rho0"<<std::endl;
 143     else report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda rho0"<<std::endl;
 144   }
 145   else if (Id_daug2==EvtPDL::getId("omega") && Vtype==VID::OMEGA)
 146   {
 147     if (!antiparticle) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : Lambda_b0 -> Lambda omega"<<std::endl;
 148     else report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda omega"<<std::endl;
 149   }
 150   else if ((Id_daug2==EvtPDL::getId("omega") ||  Id_daug2==EvtPDL::getId("rho0") )&& Vtype==VID::RHO_OMEGA_MIXING)
 151   {
 152      if (!antiparticle) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : "
 153                                                    <<"Lambda_b0 -> Lambda rho-omega-mixing"<<std::endl;
 154     else report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : "
 155                                     <<"anti-Lambda_b0 -> anti-Lambda rho-omega-mixing"<<std::endl;
 156   }
 157
 158   else
 159   {
 160     report(ERROR,fname.c_str())<<" Daughter2 is not a J/psi, phi or rho0 but a "
 161                           <<EvtPDL::name(Id_daug2)<<std::endl;
 162     abort();
 163   }
 164
 165   //fix dynamics parameters
 166   B = (double) getArg(0);
 167   C = EvtComplex((sqrt(2.)/2.),(sqrt(2.)/2.));
 168   switch(Vtype)
 169   {
 170     case VID::JPSI :             A = 0.490; break;
 171     case VID::RHO :
 172     case VID::OMEGA :
 173     case VID::RHO_OMEGA_MIXING : A = 0.194; break;
 174     default :                    A = 0;     break;
 175   }
 176   report(DEBUG,fname.c_str())<<" LambdaB decay parameters : "<<std::endl;
 177   report(DEBUG,fname.c_str())<<"   - lambda asymmetry A = "<<A<<std::endl;
 178   report(DEBUG,fname.c_str())<<"   - lambdab polarisation B = "<<B<<std::endl;
 179   report(DEBUG,fname.c_str())<<"   - lambdab density matrix rho+- C = "<<C<<std::endl;
 180
 181
 182
 183
 184 }
 185
 186
 187 //------------------------------------------------------------------------
 188 // Method 'decay'
 189 //------------------------------------------------------------------------
 190 void EvtLambdaB2LambdaV::decay( EvtParticle *lambdab)
 191 {
 192   lambdab->makeDaughters(getNDaug(),getDaugs());
 193
 194   //get lambda and V particles
 195   EvtParticle * lambda = lambdab->getDaug(0);
 196   EvtParticle * V      = lambdab->getDaug(1);
 197
 198   //get resonance masses
 199   // - LAMBDAB          -> mass given by the preceding class
 200   // - LAMBDA           -> nominal mass
 201   // - V                -> getVmass method
 202   double MASS_LAMBDAB   = lambdab->mass();
 203   double MASS_LAMBDA    = EvtPDL::getMeanMass(EvtPDL::getId("Lambda0"));
 204   double MASS_V         = getVMass(MASS_LAMBDAB,MASS_LAMBDA);
 205
 206   //generate random angles
 207   double phi   = EvtRandom::Flat(0,2*EvtConst::pi);
 208   double theta = acos( EvtRandom::Flat(-1,+1));
 209   report(DEBUG,fname.c_str())<<" Angular angles  : theta = "<<theta
 210                                            <<" ; phi = "<<phi<<std::endl;
 211   //computate resonance quadrivectors
 212   double E_lambda = (MASS_LAMBDAB*MASS_LAMBDAB + MASS_LAMBDA*MASS_LAMBDA - MASS_V*MASS_V)
 213                     /(2*MASS_LAMBDAB);
 214   double E_V      = (MASS_LAMBDAB*MASS_LAMBDAB + MASS_V*MASS_V - MASS_LAMBDA*MASS_LAMBDA)
 215                     /(2*MASS_LAMBDAB);
 216   double P = sqrt(E_lambda*E_lambda-lambda->mass()*lambda->mass());
 217
 218
 219
 220
 221   EvtVector4R P_lambdab=lambdab->getP4();
 222
 223     double px = P_lambdab.get(1);
 224     double py = P_lambdab.get(2);
 225     double pz = P_lambdab.get(3);
 226     double E  = P_lambdab.get(0);
 227     report(INFO,fname.c_str())<<"E of lambdab:  "<< P_lambdab.get(0)<<std::endl;
 228     report(INFO,fname.c_str())<<"E of lambdab:  "<< E<<std::endl;
 229
 230
 231     EvtVector4R q_lambdab2 (E,
 232                             ((1/(sqrt(pow(px,2)+pow(py,2))))*((px*(px))+(py*(py)))),
 233                             ((1/(sqrt(pow(px,2)+pow(py,2))))*(-((py)*(px))+(px*(py)))),
 234                             (pz));
 235
 236     EvtVector4R q_lambdab3 (E,
 237                             q_lambdab2.get(3),
 238                             q_lambdab2.get(1),
 239                             q_lambdab2.get(2));
 240
 241
 242     EvtVector4R q_lambda0 (E_lambda,
 243                            P*sin(theta)*cos(phi),
 244                            P*sin(theta)*sin(phi),
 245                            P*cos(theta) );
 246
 247     EvtVector4R q_V0      (E_V,
 248                            -P*sin(theta)*cos(phi),
 249                            -P*sin(theta)*sin(phi),
 250                            -P*cos(theta) );
 251
 252
 253     EvtVector4R q_lambda1 (E_lambda,
 254                            q_lambda0.get(2),
 255                            q_lambda0.get(3),
 256                            q_lambda0.get(1) );
 257
 258     EvtVector4R q_V1      (E_V,
 259                            q_V0.get(2),
 260                            q_V0.get(3),
 261                            q_V0.get(1) );
 262
 263      EvtVector4R q_lambda (E_lambda,
 264                           ((1/(sqrt(pow(px,2)+pow(py,2))))*((px*(q_lambda1.get(1))) - (py*(q_lambda1.get(2))))),
 265                           ((1/(sqrt(pow(px,2)+pow(py,2))))*((py*(q_lambda1.get(1))) + (px*(q_lambda1.get(2))))),
 266                           (q_lambda1.get(3)));
 267
 268
 269     EvtVector4R q_V     (E_V,
 270                           ((1/(sqrt(pow(px,2)+pow(py,2))))*((px*(q_V1.get(1))) - (py*(q_V1.get(2))))),
 271                           ((1/(sqrt(pow(px,2)+pow(py,2))))*((py*(q_V1.get(1))) + (px*(q_V1.get(2))))),
 272                           (q_V1.get(3)));
 273
 274
 275
 276
 277
 278    lambda->getP4();
 279    V->getP4();
 280    report(INFO,fname.c_str())<<" LambdaB  px: "<<px<<std::endl;
 281    report(INFO,fname.c_str())<<" LambdaB  py: "<<py<<std::endl;
 282    report(INFO,fname.c_str())<<" LambdaB  pz: "<<pz<<std::endl;
 283    report(INFO,fname.c_str())<<" LambdaB  E: "<<E<<std::endl;
 284
 285    report(INFO,fname.c_str())<<" Lambdab3  E:  "<<q_lambdab3.get(0)<<std::endl;
 286    report(INFO,fname.c_str())<<" Lambda 0 px:  "<<q_lambda0.get(1)<<std::endl;
 287    report(INFO,fname.c_str())<<" Lambda 0 py:  "<<q_lambda0.get(2)<<std::endl;
 288    report(INFO,fname.c_str())<<" Lambda 0 pz:  "<<q_lambda0.get(3)<<std::endl;
 289    report(INFO,fname.c_str())<<" Lambda 0 E:  "<<q_lambda0.get(0)<<std::endl;
 290    report(INFO,fname.c_str())<<" Lambda 1 px:  "<<q_lambda1.get(1)<<std::endl;
 291    report(INFO,fname.c_str())<<" Lambda 1 py:  "<<q_lambda1.get(2)<<std::endl;
 292    report(INFO,fname.c_str())<<" Lambda 1 pz:  "<<q_lambda1.get(3)<<std::endl;
 293    report(INFO,fname.c_str())<<" Lambda 1 E:  "<<q_lambda1.get(0)<<std::endl;
 294    report(INFO,fname.c_str())<<" Lambda  px:  "<<q_lambda.get(1)<<std::endl;
 295    report(INFO,fname.c_str())<<" Lambda  py:  "<<q_lambda.get(2)<<std::endl;
 296    report(INFO,fname.c_str())<<" Lambda  pz:  "<<q_lambda.get(3)<<std::endl;
 297    report(INFO,fname.c_str())<<" Lambda  E:  "<<q_lambda0.get(3)<<std::endl;
 298    report(INFO,fname.c_str())<<" V 0 px:  "<<q_V0.get(1)<<std::endl;
 299    report(INFO,fname.c_str())<<" V 0 py:  "<<q_V0.get(2)<<std::endl;
 300    report(INFO,fname.c_str())<<" V 0 pz:  "<<q_V0.get(3)<<std::endl;
 301    report(INFO,fname.c_str())<<" V 0 E:  "<<q_V0.get(0)<<std::endl;
 302    report(INFO,fname.c_str())<<" V 1 px:  "<<q_V1.get(1)<<std::endl;
 303    report(INFO,fname.c_str())<<" V 1 py:  "<<q_V1.get(2)<<std::endl;
 304    report(INFO,fname.c_str())<<" V 1 pz:  "<<q_V1.get(3)<<std::endl;
 305    report(INFO,fname.c_str())<<" V 1 E:  "<<q_V1.get(0)<<std::endl;
 306    report(INFO,fname.c_str())<<" V  px:  "<<q_V.get(1)<<std::endl;
 307    report(INFO,fname.c_str())<<" V  py:  "<<q_V.get(2)<<std::endl;
 308    report(INFO,fname.c_str())<<" V  pz:  "<<q_V.get(3)<<std::endl;
 309    report(INFO,fname.c_str())<<" V  E:  "<<q_V0.get(3)<<std::endl;
 310   //set quadrivectors to particles
 311   lambda ->init(getDaugs()[0],q_lambda);
 312   V      ->init(getDaugs()[1],q_V     );
 313
 314   //computate pdf
 315   double pdf = 1 + A*B*cos(theta) + 2*A*real(C*EvtComplex(cos(phi),sin(phi)))*sin(theta);
 316
 317   report(DEBUG,fname.c_str())<<" LambdaB decay pdf value : "<<pdf<<std::endl;
 318   //set probability
 319   setProb(pdf);
 320
 321   return;
 322 }
 323
 324
 325 //------------------------------------------------------------------------
 326 // Method 'BreitWignerRelPDF'
 327 //------------------------------------------------------------------------
 328 double EvtLambdaB2LambdaV::BreitWignerRelPDF(double m,double _m0, double _g0)
 329 {
 330   double a = (_m0 * _g0) * (_m0 * _g0);
 331   double b = (m*m - _m0*_m0)*(m*m - _m0*_m0);
 332   return a/(b+a);
 333 }
 334
 335
 336 //------------------------------------------------------------------------
 337 // Method 'RhoOmegaMixingPDF'
 338 //------------------------------------------------------------------------
 339 double EvtLambdaB2LambdaV::RhoOmegaMixingPDF(double m, double _mr, double _gr, double _mo, double _go)
 340 {
 341   double a     = m*m - _mr*_mr;
 342   double b     = m*m - _mo*_mo;
 343   double c     = _gr*_mr;
 344   double d     = _go*_mo;
 345   double re_pi = -3500e-6; //GeV^2
 346   double im_pi = -300e-6;  //GeV^2
 347   double va_pi = re_pi+im_pi;
 348
 349   //computate pdf value
 350   double f =  1/(a*a+c*c) * ( 1+
 351               (va_pi*va_pi+2*b*re_pi+2*d*im_pi)/(b*b+d*d));
 352
 353   //computate pdf max value
 354   a = 0;
 355   b = _mr*_mr - _mo*_mo;
 356
 357   double maxi = 1/(a*a+c*c) * ( 1+
 358               (va_pi*va_pi+2*b*re_pi+2*d*im_pi)/(b*b+d*d));
 359
 360   return f/maxi;
 361 }
 362
 363
 364 //------------------------------------------------------------------------
 365 // Method 'getVMass'
 366 //------------------------------------------------------------------------
 367 double EvtLambdaB2LambdaV::getVMass(double MASS_LAMBDAB, double MASS_LAMBDA)
 368 {
 369   //JPSI case
 370   if (Vtype==VID::JPSI)
 371   {
 372     return EvtPDL::getMass(EvtPDL::getId("J/psi"));
 373   }
 374
 375   //RHO,OMEGA,MIXING case
 376   else
 377   {
 378     //parameters
 379     double MASS_RHO     = EvtPDL::getMeanMass(EvtPDL::getId("rho0"));
 380     double MASS_OMEGA   = EvtPDL::getMeanMass(EvtPDL::getId("omega"));
 381     double WIDTH_RHO    = EvtPDL::getWidth(EvtPDL::getId("rho0"));
 382     double WIDTH_OMEGA  = EvtPDL::getWidth(EvtPDL::getId("omega"));
 383     double MASS_PION    = EvtPDL::getMeanMass(EvtPDL::getId("pi-"));
 384     double _max         = MASS_LAMBDAB - MASS_LAMBDA;
 385     double _min         = 2*MASS_PION;
 386
 387     double mass=0; bool test=false; int ntimes=10000;
 388     do
 389     {
 390       double y   = EvtRandom::Flat(0,1);
 391
 392       //generate mass
 393       mass = (_max-_min) * EvtRandom::Flat(0,1) + _min;
 394
 395       //pdf calculate
 396       double f=0;
 397       switch(Vtype)
 398       {
 399         case VID::RHO :              f = BreitWignerRelPDF(mass,MASS_RHO,WIDTH_RHO);     break;
 400         case VID::OMEGA :            f = BreitWignerRelPDF(mass,MASS_OMEGA,WIDTH_OMEGA); break;
 401         case VID::RHO_OMEGA_MIXING : f = RhoOmegaMixingPDF(mass,MASS_RHO,WIDTH_RHO,
 402                                                                 MASS_OMEGA,WIDTH_OMEGA); break;
 403         default :                    f = 1;                                              break;
 404       }
 405
 406       ntimes--;
 407       if (y<f) test=true;
 408     }while(ntimes && !test);
 409
 410   //looping 10000 times
 411   if (ntimes==0)
 412   {
 413       report(INFO,fname.c_str()) << "Tried accept/reject:10000"
 414                            <<" times, and rejected all the times!"<<std::endl;
 415       report(INFO,fname.c_str()) << "Is therefore accepting the last event!"<<std::endl;
 416   }
 417
 418   //return V particle mass
 419   return mass;
 420   }
 421 }
 422
 423
 424
 425
 426
 427 //************************************************************************
 428 //*                                                                      *
 429 //*                Class EvtLambda2PPiForLambdaB2LambdaV                 *
 430 //*                                                                      *
 431 //************************************************************************
 432
 433
 434 //------------------------------------------------------------------------
 435 // Constructor
 436 //------------------------------------------------------------------------
 437 EvtLambda2PPiForLambdaB2LambdaV::EvtLambda2PPiForLambdaB2LambdaV()
 438 {
 439   //set facility name
 440   fname="EvtGen.EvtLambda2PPiForLambdaB2LambdaV";
 441 }
 442
 443
 444 //------------------------------------------------------------------------
 445 // Destructor
 446 //------------------------------------------------------------------------
 447 EvtLambda2PPiForLambdaB2LambdaV::~EvtLambda2PPiForLambdaB2LambdaV()
 448 {
 449 }
 450
 451
 452 //------------------------------------------------------------------------
 453 // Method 'getName'
 454 //------------------------------------------------------------------------
 455 std::string EvtLambda2PPiForLambdaB2LambdaV::getName()
 456 {
 457   return "LAMBDA2PPIFORLAMBDAB2LAMBDAV";
 458 }
 459
 460
 461 //------------------------------------------------------------------------
 462 // Method 'clone'
 463 //------------------------------------------------------------------------
 464 EvtDecayBase* EvtLambda2PPiForLambdaB2LambdaV::clone()
 465 {
 466   return new EvtLambda2PPiForLambdaB2LambdaV;
 467 }
 468
 469
 470 //------------------------------------------------------------------------
 471 // Method 'initProbMax'
 472 //------------------------------------------------------------------------
 473 void EvtLambda2PPiForLambdaB2LambdaV::initProbMax()
 474 {
 475   //maximum (case where D is real)
 476   double Max=0;
 477   if (A==0) Max=1;
 478   else if (C==0 || real(D)==0) Max=1+fabs(A*B);
 479   else if (B==0) Max=1+ EvtConst::pi/2.0*fabs(C*A*real(D));
 480   else
 481     {
 482       double theta_max = atan(- EvtConst::pi/2.0*C*real(D)/B);
 483       double max1 = 1 + fabs(A*B*cos(theta_max)
 484                       - EvtConst::pi/2.0*C*A*real(D)*sin(theta_max));
 485       double max2 = 1 + fabs(A*B);
 486       if (max1>max2) Max=max1; else Max=max2;
 487     }
 488   report(DEBUG,fname.c_str())<<" PDF max value : "<<Max<<std::endl;
 489   setProbMax(Max);
 490 }
 491
 492
 493 //------------------------------------------------------------------------
 494 // Method 'init'
 495 //------------------------------------------------------------------------
 496 void EvtLambda2PPiForLambdaB2LambdaV::init()
 497 {
 498   bool antiparticle=false;
 499
 500   //introduction
 501   report(DEBUG,fname.c_str())<<" ***********************************************************"<<std::endl;
 502   report(DEBUG,fname.c_str())<<" *   Event Model Class : EvtLambda2PPiForLambdaB2LambdaV   *"<<std::endl;
 503   report(DEBUG,fname.c_str())<<" ***********************************************************"<<std::endl;
 504
 505   //check the number of arguments
 506   checkNArg(2);
 507
 508   //check the number of daughters
 509   checkNDaug(2);
 510
 511   const EvtId Id_mother = getParentId();
 512   const EvtId Id_daug1  = getDaug(0);
 513   const EvtId Id_daug2  = getDaug(1);
 514
 515   //lambda  identification
 516   if (Id_mother==EvtPDL::getId("Lambda0"))
 517   {
 518     antiparticle=false;
 519   }
 520   else if (Id_mother==EvtPDL::getId("anti-Lambda0"))
 521   {
 522     antiparticle=true;
 523   }
 524   else
 525   {
 526     report(ERROR,fname.c_str())<<" Mother is not a Lambda0 or an anti-Lambda0, but a "
 527                           <<EvtPDL::name(Id_mother)<<std::endl;
 528     abort();
 529   }
 530
 531   //proton identification
 532   if ( !(Id_daug1==EvtPDL::getId("p+") && !antiparticle ) && !(Id_daug1==EvtPDL::getId("anti-p-") && antiparticle) )
 533   {
 534     if (!antiparticle)
 535     {
 536       report(ERROR,fname.c_str()) << " Daughter1 is not a p+, but a "
 537                                                 << EvtPDL::name(Id_daug1)<<std::endl;
 538     }
 539     else
 540     { report(ERROR,fname.c_str()) << " Daughter1 is not an anti-p-, but a "
 541                                                 << EvtPDL::name(Id_daug1)<<std::endl;
 542     }
 543     abort();
 544   }
 545
 546   //pion identification
 547   if ( !(Id_daug2==EvtPDL::getId("pi-") && !antiparticle ) && !(Id_daug2==EvtPDL::getId("pi+") && antiparticle) )
 548   {
 549     if (!antiparticle)
 550     {
 551       report(ERROR,fname.c_str()) << " Daughter2 is not a p-, but a "
 552                                                 << EvtPDL::name(Id_daug1)<<std::endl;
 553     }
 554     else
 555     { report(ERROR,fname.c_str()) << " Daughter2 is not an p+, but a "
 556                                                 << EvtPDL::name(Id_daug1)<<std::endl;
 557     }
 558     abort();
 559   }
 560   if (!antiparticle) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : Lambda0 -> p+ pi-"<<std::endl;
 561   else report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : Anti-Lambda0 -> anti-p- pi+"<<std::endl;
 562
 563   //identification meson V
 564   if (getArg(1)==1)
 565   {
 566     Vtype=VID::JPSI;
 567     if (!antiparticle) report(DEBUG,fname.c_str())<<" From : Lambda_b0 -> Lambda J/psi"<<std::endl;
 568     else report(DEBUG,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda J/psi"<<std::endl;
 569   }
 570   else if (getArg(1)==2)
 571   {
 572     Vtype=VID::RHO;
 573     if (!antiparticle) report(DEBUG,fname.c_str())<<" From : Lambda_b0 -> Lambda rho0"<<std::endl;
 574     else report(DEBUG,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda rho0"<<std::endl;
 575   }
 576   else if (getArg(1)==3)
 577   {
 578     Vtype=VID::OMEGA;
 579     if (!antiparticle) report(DEBUG,fname.c_str())<<" From : Lambda_b0 -> Lambda omega"<<std::endl;
 580     else report(DEBUG,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda omega"<<std::endl;
 581   }
 582   else if (getArg(1)==4)
 583   {
 584     Vtype=VID::RHO_OMEGA_MIXING;
 585   }
 586   else
 587   {
 588     report(ERROR,fname.c_str()) << " Vtype " <<getArg(1)<<" is unknown"<<std::endl;
 589     if (!antiparticle) report(DEBUG,fname.c_str())<<" From : Lambda_b0 -> Lambda rho-omega-mixing"<<std::endl;
 590     else report(DEBUG,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda rho-omega-mixing"<<std::endl;    abort();
 591   }
 592
 593   //constants
 594   A = 0.642;
 595   C = (double) getArg(0);
 596   switch(Vtype)
 597   {
 598     case VID::JPSI :             B = -0.167; D = EvtComplex(0.25,0.0); break;
 599     case VID::RHO :
 600     case VID::OMEGA :
 601     case VID::RHO_OMEGA_MIXING : B = -0.21;  D = EvtComplex(0.31,0.0); break;
 602     default :                    B = 0;      D = EvtComplex(0,0);      break;
 603   }
 604
 605
 606   report(DEBUG,fname.c_str())<<" Lambda decay parameters : "<<std::endl;
 607   report(DEBUG,fname.c_str())<<"   - proton asymmetry A = "<<A<<std::endl;
 608   report(DEBUG,fname.c_str())<<"   - lambda polarisation B = "<<B<<std::endl;
 609   report(DEBUG,fname.c_str())<<"   - lambdaB polarisation C = "<<C<<std::endl;
 610   report(DEBUG,fname.c_str())<<"   - lambda density matrix rho+- D = "<<D<<std::endl;
 611 }
 612
 613
 614 //------------------------------------------------------------------------
 615 // Method 'decay'
 616 //------------------------------------------------------------------------
 617 void EvtLambda2PPiForLambdaB2LambdaV::decay( EvtParticle *lambda )
 618 {
 619   lambda->makeDaughters(getNDaug(),getDaugs());
 620
 621   //get proton and pion particles
 622   EvtParticle * proton = lambda->getDaug(0);
 623   EvtParticle * pion = lambda->getDaug(1);
 624
 625   //get resonance masses
 626   // - LAMBDA        -> mass given by EvtLambdaB2LambdaV class
 627   // - PROTON & PION -> nominal mass
 628   double MASS_LAMBDA    = lambda->mass();
 629   double MASS_PROTON    = EvtPDL::getMeanMass(EvtPDL::getId("p+"));
 630   double MASS_PION      = EvtPDL::getMeanMass(EvtPDL::getId("pi-"));
 631
 632   //generate random angles
 633   double phi   = EvtRandom::Flat(0,2*EvtConst::pi);
 634   double theta = acos( EvtRandom::Flat(-1,+1));
 635   report(DEBUG,fname.c_str())<<" Angular angles  : theta = "<<theta<<" ; phi = "<<phi<<std::endl;
 636
 637   //computate resonance quadrivectors
 638   double E_proton = (MASS_LAMBDA*MASS_LAMBDA + MASS_PROTON*MASS_PROTON - MASS_PION*MASS_PION)
 639                     /(2*MASS_LAMBDA);
 640   double E_pion   = (MASS_LAMBDA*MASS_LAMBDA + MASS_PION*MASS_PION - MASS_PROTON*MASS_PROTON)
 641                     /(2*MASS_LAMBDA);
 642   double P = sqrt(E_proton*E_proton-proton->mass()*proton->mass());
 643
 644   EvtVector4R P_lambda=lambda->getP4();
 645   EvtParticle *Mother_lambda=lambda->getParent();
 646   EvtVector4R lambdab=Mother_lambda->getP4();
 647
 648
 649
 650   double E_lambda =P_lambda.get(0);
 651   double px_M     =lambdab.get(1);
 652   double py_M     =lambdab.get(2);
 653   double pz_M     =lambdab.get(3);
 654   double E_M     =lambdab.get(0);
 655
 656   EvtVector4R q_lambdab2 (E_M,
 657                             ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(px_M))+(py_M*(py_M)))),
 658                             ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-((py_M)*(px_M))+(px_M*(py_M)))),
 659                             (pz_M));
 660
 661   EvtVector4R q_lambdab3 (E_M,
 662                             q_lambdab2.get(3),
 663                             q_lambdab2.get(1),
 664                             q_lambdab2.get(2));
 665
 666
 667
 668   EvtVector4R q_lambda1 (E_lambda,
 669                          ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(P_lambda.get(1))) + (py_M*(P_lambda.get(2))))),
 670                          ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-(py_M*(P_lambda.get(1))) + (px_M*(P_lambda.get(2))))),
 671                          P_lambda.get(3));
 672
 673   EvtVector4R q_lambda2 (E_lambda,
 674                          q_lambda1.get(3),
 675                          q_lambda1.get(1),
 676                          q_lambda1.get(2));
 677
 678
 679
 680
 681
 682    double px=q_lambda2.get(1);
 683   double py=q_lambda2.get(2);
 684   double pz=q_lambda2.get(3);
 685
 686
 687
 688
 689     EvtVector4R q_lambda4 (q_lambda2.get(0),
 690                           ((1/(sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2) + pow(q_lambda2.get(3),2))))* (1/(sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2))))*((q_lambda2.get(1))*(q_lambda2.get(1))*(q_lambda2.get(3))+((q_lambda2.get(2))*(q_lambda2.get(2))*(q_lambda2.get(3))) - ((q_lambda2.get(3))*(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2))))),
 691                           ((((q_lambda2.get(2))*(q_lambda2.get(1)))-((q_lambda2.get(1))*(q_lambda2.get(2))))/(sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2)))),
 692                           (((1/sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2) + pow(q_lambda2.get(3),2)))*( ((q_lambda2.get(1))*(q_lambda2.get(1))) +((q_lambda2.get(2))*(q_lambda2.get(2))) +   ((q_lambda2.get(3))*(q_lambda2.get(3)))))) );
 693
 694
 695
 696
 697    EvtVector4R q_proton1 (E_proton,
 698                         P*sin(theta)*cos(phi),
 699                         P*sin(theta)*sin(phi),
 700                         P*cos(theta));
 701    EvtVector4R q_pion1   (E_pion,
 702                         -P*sin(theta)*cos(phi),
 703                         -P*sin(theta)*sin(phi),
 704                         -P*cos(theta));
 705
 706
 707    EvtVector4R q_proton3     (E_proton,
 708                        ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_proton1.get(1))*(px)*(pz) - ((q_proton1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_proton1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
 709                        (((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))* (1/(sqrt(pow(px,2) + pow(py,2))))*(((q_proton1.get(1)))*(py)*(pz) + ((q_proton1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_proton1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
 710                         (((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))*((-(q_proton1.get(1)))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_proton1.get(3))*(pz)))));
 711
 712  EvtVector4R q_pion3     (E_pion,
 713                          ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_pion1.get(1))*(px)*(pz) - ((q_pion1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_pion1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
 714                          (((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_pion1.get(1))*(py)*(pz) + ((q_pion1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_pion1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
 715                          ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))*((-(q_pion1.get(1)))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_pion1.get(3))*(pz)))));
 716
 717  EvtVector4R q_proton5 (q_proton3.get(0),
 718                         (q_proton3.get(2)),
 719                         (q_proton3.get(3)),
 720                         (q_proton3.get(1)));
 721
 722  EvtVector4R q_pion5      (q_pion3.get(0),
 723                            (q_pion3.get(2)),
 724                            (q_pion3.get(3)),
 725                            (q_pion3.get(1)));
 726
 727  EvtVector4R q_proton          (q_proton5.get(0),
 728                                 ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_proton5.get(1)))-(py_M*(q_proton5.get(2))))),
 729                                 ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_proton5.get(1)))+(px_M*(q_proton5.get(2))))),
 730                                 (q_proton5.get(3)));
 731
 732
 733  EvtVector4R q_pion      (q_pion5.get(0),
 734                           ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_pion5.get(1)))-(py_M*(q_pion5.get(2))))),
 735                           ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_pion5.get(1)))+(px_M*(q_pion5.get(2))))),
 736                           (q_pion5.get(3)));
 737
 738 report(INFO,fname.c_str())<<" Lambdab  px: "<<px_M<<std::endl;
 739 report(INFO,fname.c_str())<<" Lambdab  py: "<<py_M<<std::endl;
 740 report(INFO,fname.c_str())<<" Lambdab  pz: "<<pz_M<<std::endl;
 741 report(INFO,fname.c_str())<<" Lambdab  E: "<<E_M<<std::endl;
 742 report(INFO,fname.c_str())<<" Lambdab2  px:  "<<q_lambdab2.get(1)<<std::endl;
 743 report(INFO,fname.c_str())<<" Lambdab2  py:  "<<q_lambdab2.get(2)<<std::endl;
 744 report(INFO,fname.c_str())<<" Lambdab2  pz:  "<<q_lambdab2.get(3)<<std::endl;
 745 report(INFO,fname.c_str())<<" Lambdab2  E:   "<<q_lambdab2.get(0)<<std::endl;
 746 report(INFO,fname.c_str())<<" Lambdab3  px:  "<<q_lambdab3.get(1)<<std::endl;
 747 report(INFO,fname.c_str())<<" Lambdab3  py:  "<<q_lambdab3.get(2)<<std::endl;
 748 report(INFO,fname.c_str())<<" Lambdab3  pz:  "<<q_lambdab3.get(3)<<std::endl;
 749 report(INFO,fname.c_str())<<" Lambdab3  E:   "<<q_lambdab3.get(0)<<std::endl;
 750 report(INFO,fname.c_str())<<" Lambda 0  px:  "<<P_lambda.get(1)<<std::endl;
 751 report(INFO,fname.c_str())<<" Lambda 0  py:  "<<P_lambda.get(2)<<std::endl;
 752 report(INFO,fname.c_str())<<" Lambda 0  pz:  "<<P_lambda.get(3)<<std::endl;
 753 report(INFO,fname.c_str())<<" Lambda 0  E:   "<<P_lambda.get(0)<<std::endl;
 754 report(INFO,fname.c_str())<<" Lambda 1  px:  "<<q_lambda1.get(1)<<std::endl;
 755 report(INFO,fname.c_str())<<" Lambda 1  py:  "<<q_lambda1.get(2)<<std::endl;
 756 report(INFO,fname.c_str())<<" Lambda 1  pz:  "<<q_lambda1.get(3)<<std::endl;
 757 report(INFO,fname.c_str())<<" Lambda 1  E:   "<<q_lambda1.get(0)<<std::endl;
 758 report(INFO,fname.c_str())<<" Lambda 2  px:  "<<q_lambda2.get(1)<<std::endl;
 759 report(INFO,fname.c_str())<<" Lambda 2  py:  "<<q_lambda2.get(2)<<std::endl;
 760 report(INFO,fname.c_str())<<" Lambda 2  pz:  "<<q_lambda2.get(3)<<std::endl;
 761 report(INFO,fname.c_str())<<" Lambda 2  E:   "<<q_lambda2.get(0)<<std::endl;
 762
 763 report(INFO,fname.c_str())<<" Lambda  px: "<<px<<std::endl;
 764 report(INFO,fname.c_str())<<" Lambda  py: "<<py<<std::endl;
 765 report(INFO,fname.c_str())<<" Lambda  pz: "<<pz<<std::endl;
 766
 767  report(INFO,fname.c_str())<<" pion 1 px:  "<<q_pion1.get(1)<<std::endl;
 768  report(INFO,fname.c_str())<<" pion 1 py:  "<<q_pion1.get(2)<<std::endl;
 769  report(INFO,fname.c_str())<<" pion 1 pz:  "<<q_pion1.get(3)<<std::endl;
 770  report(INFO,fname.c_str())<<" pion 1 E:   "<<q_pion1.get(0)<<std::endl;
 771
 772 report(INFO,fname.c_str())<<" pion 3 px:  "<<q_pion3.get(1)<<std::endl;
 773 report(INFO,fname.c_str())<<" pion 3 px:  "<<q_pion3.get(1)<<std::endl;
 774 report(INFO,fname.c_str())<<" pion 3 py:  "<<q_pion3.get(2)<<std::endl;
 775 report(INFO,fname.c_str())<<" pion 3 pz:  "<<q_pion3.get(3)<<std::endl;
 776 report(INFO,fname.c_str())<<" pion 3 E:   "<<q_pion3.get(0)<<std::endl;
 777
 778 report(INFO,fname.c_str())<<" pion 5 px:  "<<q_pion5.get(1)<<std::endl;
 779 report(INFO,fname.c_str())<<" pion 5 py:  "<<q_pion5.get(2)<<std::endl;
 780  report(INFO,fname.c_str())<<" pion 5 pz:  "<<q_pion5.get(3)<<std::endl;
 781  report(INFO,fname.c_str())<<" pion 5 E:   "<<q_pion5.get(0)<<std::endl;
 782
 783
 784
 785  report(INFO,fname.c_str())<<" proton 1  px:  "<<q_proton1.get(1)<<std::endl;
 786  report(INFO,fname.c_str())<<" proton 1  py:  "<<q_proton1.get(2)<<std::endl;
 787  report(INFO,fname.c_str())<<" proton 1  pz:  "<<q_proton1.get(3)<<std::endl;
 788  report(INFO,fname.c_str())<<" proton 1  E:   "<<q_proton1.get(0)<<std::endl;
 789
 790 report(INFO,fname.c_str())<<" proton 3  px:  "<<q_proton3.get(1)<<std::endl;
 791  report(INFO,fname.c_str())<<" proton 3  py:  "<<q_proton3.get(2)<<std::endl;
 792  report(INFO,fname.c_str())<<" proton 3  pz:  "<<q_proton3.get(3)<<std::endl;
 793  report(INFO,fname.c_str())<<" proton 3  E:   "<<q_proton3.get(0)<<std::endl;
 794
 795 report(INFO,fname.c_str())<<" proton 5  px:  "<<q_proton5.get(1)<<std::endl;
 796  report(INFO,fname.c_str())<<" proton 5  py:  "<<q_proton5.get(2)<<std::endl;
 797  report(INFO,fname.c_str())<<" proton 5  pz:  "<<q_proton5.get(3)<<std::endl;
 798  report(INFO,fname.c_str())<<" proton 5  E:   "<<q_proton5.get(0)<<std::endl;
 799
 800
 801 report(INFO,fname.c_str())<<" proton  px:  "<<q_proton.get(1)<<std::endl;
 802    report(INFO,fname.c_str())<<" proton  py:  "<<q_proton.get(2)<<std::endl;
 803    report(INFO,fname.c_str())<<"proton  pz:  "<<q_proton.get(3)<<std::endl;
 804  report(INFO,fname.c_str())<<" pion px:  "<<q_pion.get(1)<<std::endl;
 805    report(INFO,fname.c_str())<<" pion py:  "<<q_pion.get(2)<<std::endl;
 806    report(INFO,fname.c_str())<<" pion pz:  "<<q_pion.get(3)<<std::endl;
 807
 808
 809
 810
 811
 812   ;
 813
 814
 815
 816
 817
 818
 819  ///////////*******NEW********//////////////////////
 820
 821   //set quadrivectors to particles
 822   proton->init(getDaugs()[0],q_proton);
 823   pion  ->init(getDaugs()[1],q_pion  );
 824
 825   //computate pdf
 826   //double pdf = 1 + A*B*cos(theta) - EvtConst::pi/2.0*C*A*real(D*EvtComplex(cos(phi),sin(phi)))*sin(theta);
 827   double pdf = 1 + A*B*cos(theta) + 2*A*real(D*EvtComplex(cos(phi),sin(phi)))*sin(theta);
 828   report(DEBUG,fname.c_str())<<" Lambda decay pdf value : "<<pdf<<std::endl;
 829   //set probability
 830   setProb(pdf);
 831
 832   return;
 833 }
 834
 835
 836
 837
 838 //************************************************************************
 839 //*                                                                      *
 840 //*                Class EvtLambda2PPiForLambdaB2LambdaV                 *
 841 //*                                                                      *
 842 //************************************************************************
 843
 844
 845 //------------------------------------------------------------------------
 846 // Constructor
 847 //------------------------------------------------------------------------
 848 EvtV2VpVmForLambdaB2LambdaV::EvtV2VpVmForLambdaB2LambdaV()
 849 {
 850   //set facility name
 851   fname="EvtGen.EvtV2VpVmForLambdaB2LambdaV";
 852 }
 853
 854
 855 //------------------------------------------------------------------------
 856 // Destructor
 857 //------------------------------------------------------------------------
 858 EvtV2VpVmForLambdaB2LambdaV::~EvtV2VpVmForLambdaB2LambdaV()
 859 {}
 860
 861
 862 //------------------------------------------------------------------------
 863 // Method 'getName'
 864 //------------------------------------------------------------------------
 865 std::string EvtV2VpVmForLambdaB2LambdaV::getName()
 866 {
 867   return "V2VPVMFORLAMBDAB2LAMBDAV";
 868 }
 869
 870 //------------------------------------------------------------------------
 871 // Method 'clone'
 872 //------------------------------------------------------------------------
 873 EvtDecayBase* EvtV2VpVmForLambdaB2LambdaV::clone()
 874 {
 875   return new EvtV2VpVmForLambdaB2LambdaV;
 876 }
 877
 878
 879 //------------------------------------------------------------------------
 880 // Method 'initProbMax'
 881 //------------------------------------------------------------------------
 882 void EvtV2VpVmForLambdaB2LambdaV::initProbMax()
 883 {
 884   //maximum
 885   double Max = 0;
 886   if (Vtype==VID::JPSI)
 887   {
 888     if ((1-3*A)>0) Max=2*(1-A);
 889     else Max=1+A;
 890   }
 891   else
 892   {
 893     if ((3*A-1)>=0) Max=2*A;
 894     else Max=1-A;
 895   }
 896
 897   report(DEBUG,fname.c_str())<<" PDF max value : "<<Max<<std::endl;
 898   setProbMax(Max);
 899 }
 900
 901
 902 //------------------------------------------------------------------------
 903 // Method 'init'
 904 //------------------------------------------------------------------------
 905 void EvtV2VpVmForLambdaB2LambdaV::init()
 906 {
 907   //introduction
 908   report(DEBUG,fname.c_str())<<" ***********************************************************"<<std::endl;
 909   report(DEBUG,fname.c_str())<<" *     Event Model Class : EvtV2VpVmForLambdaB2LambdaV     *"<<std::endl;
 910   report(DEBUG,fname.c_str())<<" ***********************************************************"<<std::endl;
 911
 912   //check the number of arguments
 913   checkNArg(2);
 914
 915   //check the number of daughters
 916   checkNDaug(2);
 917
 918   const EvtId Id_mother = getParentId();
 919   const EvtId Id_daug1  = getDaug(0);
 920   const EvtId Id_daug2  = getDaug(1);
 921
 922   //identification meson V
 923   if (getArg(1)==1) Vtype=VID::JPSI;
 924   else if (getArg(1)==2) Vtype=VID::RHO;
 925   else if (getArg(1)==3) Vtype=VID::OMEGA;
 926   else if (getArg(1)==4) Vtype=VID::RHO_OMEGA_MIXING;
 927   else
 928   {
 929     report(ERROR,fname.c_str()) << " Vtype " <<getArg(1)<<" is unknown"<<std::endl;
 930     abort();
 931   }
 932
 933   //vector meson V
 934   if (Id_mother==EvtPDL::getId("J/psi") && Vtype==VID::JPSI)
 935   {
 936   }
 937   else if (Id_mother==EvtPDL::getId("omega") &&  Vtype==VID::OMEGA)
 938   {
 939   }
 940   else if (Id_mother==EvtPDL::getId("rho0") &&  Vtype==VID::RHO)
 941   {
 942   }
 943   else if ((Id_mother==EvtPDL::getId("rho0") || Id_mother==EvtPDL::getId("omega")) && Vtype==VID::RHO_OMEGA_MIXING)
 944   {
 945   }
 946   else
 947   {
 948     report(ERROR,fname.c_str())<<" Mother is not a J/psi, phi or rho0 but a "
 949                           <<EvtPDL::name(Id_mother)<<std::endl;
 950     abort();
 951   }
 952
 953   //daughters for each V possibility
 954   switch(Vtype)
 955     {
 956     case VID::JPSI :
 957       if (Id_daug1!=EvtPDL::getId("mu+"))
 958       {
 959         report(ERROR,fname.c_str()) << " Daughter1 is not a mu+, but a "
 960                                                            << EvtPDL::name(Id_daug1)<<std::endl;
 961         abort();
 962       }
 963       if (Id_daug2!=EvtPDL::getId("mu-"))
 964       {
 965         report(ERROR,fname.c_str()) << " Daughter2 is not a mu-, but a "
 966                                                            << EvtPDL::name(Id_daug2)<<std::endl;
 967         abort();
 968       }
 969       report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : J/psi -> mu+ mu-"<<std::endl;
 970       break;
 971
 972     case VID::RHO :
 973     case VID::OMEGA :
 974     case VID::RHO_OMEGA_MIXING :
 975       if (Id_daug1!=EvtPDL::getId("pi+"))
 976       {
 977         report(ERROR,fname.c_str()) << " Daughter1 is not a pi+, but a "
 978                                                            << EvtPDL::name(Id_daug1)<<std::endl;
 979         abort();
 980       }
 981       if (Id_daug2!=EvtPDL::getId("pi-"))
 982       {
 983         report(ERROR,fname.c_str()) << " Daughter2 is not a pi-, but a "
 984                                                            << EvtPDL::name(Id_daug2)<<std::endl;
 985         abort();
 986       }
 987       if (Vtype==VID::RHO) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : rho0 -> pi+ pi-"<<std::endl;
 988       if (Vtype==VID::OMEGA) report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : omega -> pi+ pi-"<<std::endl;
 989       if (Vtype==VID::RHO_OMEGA_MIXING)
 990               report(DEBUG,fname.c_str())<<" Decay mode successfully initialized : rho-omega mixing -> pi+ pi-"<<std::endl; break;
 991
 992     default :
 993       report(ERROR,fname.c_str()) << "No decay mode chosen ! "<<std::endl;
 994       abort();
 995       break;
 996     }
 997
 998   //fix dynamics parameters
 999   switch(Vtype)
1000   {
1001   case VID::JPSI :             A = 0.66; break;
1002   case VID::RHO :
1003   case VID::OMEGA :
1004   case VID::RHO_OMEGA_MIXING : A = 0.79; break;
1005   default :                    A = 0;    break;
1006   }
1007
1008   report(DEBUG,fname.c_str())<<" V decay parameters : "<<std::endl;
1009   report(DEBUG,fname.c_str())<<"   - V density matrix rho00 A = "<<A<<std::endl;
1010
1011
1012 }
1013
1014 //------------------------------------------------------------------------
1015 // Method 'decay'
1016 //------------------------------------------------------------------------
1017 void EvtV2VpVmForLambdaB2LambdaV::decay( EvtParticle *V )
1018 {
1019   V->makeDaughters(getNDaug(),getDaugs());
1020
1021   //get Vp and Vm particles
1022   EvtParticle * Vp = V->getDaug(0);
1023   EvtParticle * Vm = V->getDaug(1);
1024
1025   //get resonance masses
1026   // - V         -> mass given by EvtLambdaB2LambdaV class
1027   // - Vp & Vm   -> nominal mass
1028   double MASS_V   = V->mass();
1029   double MASS_VM  = 0;
1030   switch(Vtype)
1031   {
1032   case VID::JPSI :             MASS_VM=EvtPDL::getMeanMass(EvtPDL::getId("mu-")); break;
1033   case VID::RHO :
1034   case VID::OMEGA :
1035   case VID::RHO_OMEGA_MIXING : MASS_VM=EvtPDL::getMeanMass(EvtPDL::getId("pi-")); break;
1036   default :                    MASS_VM=0;                                         break;
1037   }
1038   double MASS_VP  = MASS_VM;
1039
1040   //generate random angles
1041   double phi   = EvtRandom::Flat(0,2*EvtConst::pi);
1042   double theta = acos( EvtRandom::Flat(-1,+1));
1043   report(DEBUG,fname.c_str())<<" Angular angles  : theta = "<<theta<<" ; phi = "<<phi<<std::endl;
1044
1045   //computate resonance quadrivectors
1046   double E_Vp = (MASS_V*MASS_V + MASS_VP*MASS_VP - MASS_VM*MASS_VM)
1047                  /(2*MASS_V);
1048   double E_Vm = (MASS_V*MASS_V + MASS_VM*MASS_VM - MASS_VP*MASS_VP)
1049                  /(2*MASS_V);
1050   double P = sqrt(E_Vp*E_Vp-Vp->mass()*Vp->mass());
1051
1052   EvtVector4R P_V=V->getP4();
1053   EvtParticle *Mother_V=V->getParent();
1054   EvtVector4R lambdab=Mother_V->getP4();
1055
1056
1057   double E_V=(P_V.get(0));
1058    double px_M=lambdab.get(1);
1059    double py_M=lambdab.get(2);
1060    double pz_M=lambdab.get(3);
1061    double E_M=lambdab.get(0);
1062
1063    EvtVector4R q_lambdab2 (E_M,
1064                             ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(px_M))+(py_M*(py_M)))),
1065                             ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-((py_M)*(px_M))+(px_M*(py_M)))),
1066                             (pz_M));
1067
1068   EvtVector4R q_lambdab3 (E_M,
1069                             q_lambdab2.get(3),
1070                             q_lambdab2.get(1),
1071                             q_lambdab2.get(2));
1072
1073
1074  EvtVector4R q_V1 (E_V,
1075                          ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(P_V.get(1))) + (py_M*(P_V.get(2))))),
1076                          ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-(py_M*(P_V.get(1))) + (px_M*(P_V.get(2))))),
1077                          P_V.get(3));
1078
1079   EvtVector4R q_V2 (E_V,
1080                          q_V1.get(3),
1081                          q_V1.get(1),
1082                          q_V1.get(2));
1083
1084
1085
1086      double px= -(q_V2.get(1));
1087   double py=-(q_V2.get(2));
1088   double pz=-(q_V2.get(3));
1089
1090
1091
1092
1093   EvtVector4R q_V4 (q_V2.get(0),
1094                           ((1/(sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2) + pow(q_V2.get(3),2))))* (1/(sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2))))*((q_V2.get(1))*(q_V2.get(1))*(q_V2.get(3))+((q_V2.get(2))*(q_V2.get(2))*(q_V2.get(3))) - ((q_V2.get(3))*(pow(q_V2.get(1),2) + pow(q_V2.get(2),2))))),
1095                           ((((q_V2.get(2))*(q_V2.get(1)))-((q_V2.get(1))*(q_V2.get(2))))/(sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2)))),
1096                           (((1/sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2) + pow(q_V2.get(3),2)))*( ((q_V2.get(1))*(q_V2.get(1))) +((q_V2.get(2))*(q_V2.get(2))) +   ((q_V2.get(3))*(q_V2.get(3)))))) );
1097
1098
1099
1100    EvtVector4R q_Vp1     (E_Vp,
1101                         P*sin(theta)*cos(phi),
1102                         P*sin(theta)*sin(phi),
1103                         P*cos(theta));
1104    EvtVector4R q_Vm1     (E_Vm,
1105                         -P*sin(theta)*cos(phi),
1106                         -P*sin(theta)*sin(phi),
1107                         -P*cos(theta));
1108
1109     EvtVector4R q_Vp3     (q_Vp1.get(0),
1110                        ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_Vp1.get(1))*(px)*(pz)+((q_Vp1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vp1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
1111                        ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*(((q_Vp1.get(1)))*(py)*(pz) - ((q_Vp1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vp1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
1112                         ((-(1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))*((q_Vp1.get(1))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_Vp1.get(3))*(pz)))));
1113
1114    EvtVector4R q_Vm3     (q_Vm1.get(0),
1115                           ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_Vm1.get(1))*(px)*(pz)+((q_Vm1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vm1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
1116                           ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*(((q_Vm1.get(1)))*(py)*(pz) - ((q_Vm1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vm1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
1117                           ((-(1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))*((q_Vm1.get(1))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_Vm1.get(3))*(pz)))));
1118
1119
1120
1121
1122
1123
1124  EvtVector4R q_Vp5 (q_Vp3.get(0),
1125                         (q_Vp3.get(2)),
1126                         (q_Vp3.get(3)),
1127                         (q_Vp3.get(1)));
1128
1129    EvtVector4R q_Vm5      (q_Vm3.get(0),
1130                            (q_Vm3.get(2)),
1131                            (q_Vm3.get(3)),
1132                            (q_Vm3.get(1)));
1133
1134
1135  EvtVector4R q_Vp          (q_Vp5.get(0),
1136                                 ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_Vp5.get(1)))-(py_M*(q_Vp5.get(2))))),
1137                                 ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_Vp5.get(1)))+(px_M*(q_Vp5.get(2))))),
1138                                 (q_Vp5.get(3)));
1139
1140
1141  EvtVector4R q_Vm      (q_Vm5.get(0),
1142                           ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_Vm5.get(1)))-(py_M*(q_Vm5.get(2))))),
1143                           ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_Vm5.get(1)))+(px_M*(q_Vm5.get(2))))),
1144                           (q_Vm5.get(3)));
1145
1146   report(INFO,fname.c_str())<<" Lambdab  px: "<<px_M<<std::endl;
1147 report(INFO,fname.c_str())<<" Lambdab  py: "<<py_M<<std::endl;
1148 report(INFO,fname.c_str())<<" Lambdab  pz: "<<pz_M<<std::endl;
1149 report(INFO,fname.c_str())<<" Lambdab  E: "<<E_M<<std::endl;
1150 report(INFO,fname.c_str())<<" Lambdab2  px:  "<<q_lambdab2.get(1)<<std::endl;
1151 report(INFO,fname.c_str())<<" Lambdab2  py:  "<<q_lambdab2.get(2)<<std::endl;
1152 report(INFO,fname.c_str())<<" Lambdab2  pz:  "<<q_lambdab2.get(3)<<std::endl;
1153 report(INFO,fname.c_str())<<" Lambdab2  E:   "<<q_lambdab2.get(0)<<std::endl;
1154 report(INFO,fname.c_str())<<" Lambdab3  px:  "<<q_lambdab3.get(1)<<std::endl;
1155 report(INFO,fname.c_str())<<" Lambdab3  py:  "<<q_lambdab3.get(2)<<std::endl;
1156 report(INFO,fname.c_str())<<" Lambdab3  pz:  "<<q_lambdab3.get(3)<<std::endl;
1157 report(INFO,fname.c_str())<<" Lambdab3  E:   "<<q_lambdab3.get(0)<<std::endl;
1158 report(INFO,fname.c_str())<<" V 0  px:  "<<P_V.get(1)<<std::endl;
1159 report(INFO,fname.c_str())<<" V 0  py:  "<<P_V.get(2)<<std::endl;
1160 report(INFO,fname.c_str())<<" V 0  pz:  "<<P_V.get(3)<<std::endl;
1161 report(INFO,fname.c_str())<<" V 0  E:   "<<P_V.get(0)<<std::endl;
1162 report(INFO,fname.c_str())<<" V 1  px:  "<<q_V1.get(1)<<std::endl;
1163 report(INFO,fname.c_str())<<" V 1  py:  "<<q_V1.get(2)<<std::endl;
1164 report(INFO,fname.c_str())<<" V 1  pz:  "<<q_V1.get(3)<<std::endl;
1165 report(INFO,fname.c_str())<<" V 1  E:   "<<q_V1.get(0)<<std::endl;
1166 report(INFO,fname.c_str())<<" V 2  px:  "<<q_V2.get(1)<<std::endl;
1167 report(INFO,fname.c_str())<<" V 2  py:  "<<q_V2.get(2)<<std::endl;
1168 report(INFO,fname.c_str())<<" V 2  pz:  "<<q_V2.get(3)<<std::endl;
1169 report(INFO,fname.c_str())<<" V 2  E:   "<<q_V2.get(0)<<std::endl;
1170 report(INFO,fname.c_str())<<" V  px: "<<px<<std::endl;
1171 report(INFO,fname.c_str())<<" V  py: "<<py<<std::endl;
1172 report(INFO,fname.c_str())<<" V  pz: "<<pz<<std::endl;
1173  report(INFO,fname.c_str())<<" Vm 1 px:  "<<q_Vm1.get(1)<<std::endl;
1174  report(INFO,fname.c_str())<<" Vm 1 py:  "<<q_Vm1.get(2)<<std::endl;
1175  report(INFO,fname.c_str())<<" Vm 1 pz:  "<<q_Vm1.get(3)<<std::endl;
1176  report(INFO,fname.c_str())<<" Vm 1 E:   "<<q_Vm1.get(0)<<std::endl;
1177 report(INFO,fname.c_str())<<" Vm 3 px:  "<<q_Vm3.get(1)<<std::endl;
1178 report(INFO,fname.c_str())<<" Vm 3 px:  "<<q_Vm3.get(1)<<std::endl;
1179 report(INFO,fname.c_str())<<" Vm 3 py:  "<<q_Vm3.get(2)<<std::endl;
1180 report(INFO,fname.c_str())<<" Vm 3 pz:  "<<q_Vm3.get(3)<<std::endl;
1181 report(INFO,fname.c_str())<<" Vm 3 E:   "<<q_Vm3.get(0)<<std::endl;
1182 report(INFO,fname.c_str())<<" Vm 5 px:  "<<q_Vm5.get(1)<<std::endl;
1183 report(INFO,fname.c_str())<<" Vm 5 py:  "<<q_Vm5.get(2)<<std::endl;
1184  report(INFO,fname.c_str())<<" Vm 5 pz:  "<<q_Vm5.get(3)<<std::endl;
1185  report(INFO,fname.c_str())<<" Vm 5 E:   "<<q_Vm5.get(0)<<std::endl;
1186  report(INFO,fname.c_str())<<" Vp 1  px:  "<<q_Vp1.get(1)<<std::endl;
1187  report(INFO,fname.c_str())<<" Vp 1  py:  "<<q_Vp1.get(2)<<std::endl;
1188  report(INFO,fname.c_str())<<" Vp 1  pz:  "<<q_Vp1.get(3)<<std::endl;
1189  report(INFO,fname.c_str())<<" Vp 1  E:   "<<q_Vp1.get(0)<<std::endl;
1190 report(INFO,fname.c_str())<<" Vp 3  px:  "<<q_Vp3.get(1)<<std::endl;
1191  report(INFO,fname.c_str())<<" Vp 3  py:  "<<q_Vp3.get(2)<<std::endl;
1192  report(INFO,fname.c_str())<<" Vp 3  pz:  "<<q_Vp3.get(3)<<std::endl;
1193  report(INFO,fname.c_str())<<" Vp 3  E:   "<<q_Vp3.get(0)<<std::endl;
1194  report(INFO,fname.c_str())<<" Vp 5  px:  "<<q_Vp5.get(1)<<std::endl;
1195  report(INFO,fname.c_str())<<" Vp 5  py:  "<<q_Vp5.get(2)<<std::endl;
1196  report(INFO,fname.c_str())<<" Vp 5  pz:  "<<q_Vp5.get(3)<<std::endl;
1197  report(INFO,fname.c_str())<<" Vp 5  E:   "<<q_Vp5.get(0)<<std::endl;
1198  report(INFO,fname.c_str())<<" Vp  px:  "<<q_Vp.get(1)<<std::endl;
1199  report(INFO,fname.c_str())<<" Vp  py:  "<<q_Vp.get(2)<<std::endl;
1200  report(INFO,fname.c_str())<<"Vp  pz:  "<<q_Vp.get(3)<<std::endl;
1201  report(INFO,fname.c_str())<<" Vm px:  "<<q_Vm.get(1)<<std::endl;
1202  report(INFO,fname.c_str())<<" Vm py:  "<<q_Vm.get(2)<<std::endl;
1203  report(INFO,fname.c_str())<<" Vm pz:  "<<q_Vm.get(3)<<std::endl;
1204
1205   //set quadrivectors to particles
1206   Vp->init(getDaugs()[0],q_Vp);
1207   Vm->init(getDaugs()[1],q_Vm);
1208
1209   //computate pdf
1210   double pdf = 0;
1211   if (Vtype==VID::JPSI)
1212   {
1213     //leptonic case
1214      pdf = (1-3*A)*cos(theta)*cos(theta) + (1+A);
1215   }
1216   else
1217   {
1218     //hadronic case
1219     pdf = (3*A-1)*cos(theta)*cos(theta) + (1-A);
1220
1221   }
1222   report(DEBUG,fname.c_str())<<" V decay pdf value : "<<pdf<<std::endl;
1223
1224   //set probability
1225   setProb(pdf);
1226   return;
1227 }