TEvtGen/EvtGenModels/EvtBTo4piCP.cxx

   1 //--------------------------------------------------------------------------
   2 //
   3 // Environment:
   4 //      This software is part of the EvtGen package developed jointly
   5 //      for the BaBar and CLEO collaborations.  If you use all or part
   6 //      of it, please give an appropriate acknowledgement.
   7 //
   8 // Copyright Information: See EvtGen/COPYRIGHT
   9 //      Copyright (C) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtBTo4piCP.cc
  12 //
  13 // Description: Routine to decay B->pi+ pi- pi+ pi-.
  14 //
  15 // Modification history:
  16 //
  17 //    RYD     March 2, 1997         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtGenKine.hh"
  25 #include "EvtGenBase/EvtCPUtil.hh"
  26 #include "EvtGenBase/EvtPDL.hh"
  27 #include "EvtGenBase/EvtReport.hh"
  28 #include "EvtGenModels/EvtBTo4piCP.hh"
  29 #include "EvtGenBase/EvtId.hh"
  30 #include "EvtGenBase/EvtConst.hh"
  31 #include <string>
  32
  33 EvtBTo4piCP::~EvtBTo4piCP() {}
  34
  35
  36 EvtComplex EvtAmpA2(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
  37                     const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
  38
  39   //added by Lange Jan4,2000
  40   static EvtId A2M=EvtPDL::getId("a_2-");
  41   static EvtId RHO0=EvtPDL::getId("rho0");
  42
  43   EvtVector4R p4a2,p4rho,p4b;
  44
  45   p4rho=p4pi1+p4pi2;
  46
  47   p4a2=p4rho+p4pi3;
  48
  49   p4b=p4a2+p4pi4;
  50
  51   EvtVector4R p4b_a2,p4rho_a2,p4pi1_a2,p4a2_a2;
  52
  53   p4b_a2=boostTo(p4b,p4a2);
  54   p4rho_a2=boostTo(p4rho,p4a2);
  55   p4pi1_a2=boostTo(p4pi1,p4a2);
  56   p4a2_a2=boostTo(p4a2,p4a2);
  57
  58   EvtVector4R p4pi1_rho;
  59
  60   p4pi1_rho=boostTo(p4pi1_a2,p4rho_a2);
  61
  62   EvtVector4R vb,vrho,vpi,t;
  63
  64   vb=p4b_a2/p4b_a2.d3mag();
  65   vrho=p4rho_a2/p4rho_a2.d3mag();
  66   vpi=p4pi1_rho/p4pi1_rho.d3mag();
  67
  68   t.set(1.0,0.0,0.0,0.0);
  69
  70   //  EvtComplex amp_a1,amp_a2;
  71   EvtComplex amp_a2;
  72
  73   //  double bwm_a1=EvtPDL::getMeanMass(A1M);
  74   //  double gamma_a1=EvtPDL::getWidth(A1M);
  75   double bwm_a2=EvtPDL::getMeanMass(A2M);
  76   double gamma_a2=EvtPDL::getWidth(A2M);
  77   double bwm_rho=EvtPDL::getMeanMass(RHO0);
  78   double gamma_rho=EvtPDL::getWidth(RHO0);
  79
  80   amp_a2=(sqrt(gamma_a2/EvtConst::twoPi)/
  81     ((p4a2).mass()-bwm_a2-EvtComplex(0.0,0.5*gamma_a2)))*
  82          (sqrt(gamma_rho/EvtConst::twoPi)/
  83     ((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
  84
  85   return amp_a2*
  86     (vb.get(1)*vrho.get(1)+vb.get(2)*vrho.get(2)+vb.get(3)*vrho.get(3))*
  87     (
  88      vpi.get(1)*(vb.get(2)*vrho.get(3)-vb.get(3)*vrho.get(2))+
  89      vpi.get(2)*(vb.get(3)*vrho.get(1)-vb.get(1)*vrho.get(3))+
  90      vpi.get(3)*(vb.get(1)*vrho.get(2)-vb.get(2)*vrho.get(1))
  91      );
  92
  93 }
  94
  95 EvtComplex EvtAmpA1(const EvtVector4R& p4pi1,const EvtVector4R& p4pi2,
  96                     const EvtVector4R& p4pi3,const EvtVector4R& p4pi4){
  97
  98   //added by Lange Jan4,2000
  99   static EvtId A1M=EvtPDL::getId("a_1-");
 100   static EvtId RHO0=EvtPDL::getId("rho0");
 101
 102   EvtVector4R p4a1,p4rho,p4b;
 103
 104   p4rho=p4pi1+p4pi2;
 105
 106   p4a1=p4rho+p4pi3;
 107
 108   p4b=p4a1+p4pi4;
 109
 110   EvtVector4R p4b_a1,p4rho_a1,p4pi1_a1,p4a1_a1;
 111
 112   p4b_a1=boostTo(p4b,p4a1);
 113   p4rho_a1=boostTo(p4rho,p4a1);
 114   p4pi1_a1=boostTo(p4pi1,p4a1);
 115   p4a1_a1=boostTo(p4a1,p4a1);
 116
 117   EvtVector4R p4pi1_rho;
 118
 119   p4pi1_rho=boostTo(p4pi1_a1,p4rho_a1);
 120
 121   EvtVector4R vb,vrho,vpi,t;
 122
 123   vb=p4b_a1/p4b_a1.d3mag();
 124   vrho=p4rho_a1/p4rho_a1.d3mag();
 125   vpi=p4pi1_rho/p4pi1_rho.d3mag();
 126
 127   t.set(1.0,0.0,0.0,0.0);
 128
 129   EvtComplex amp_a1;
 130
 131   double bwm_a1=EvtPDL::getMeanMass(A1M);
 132   double gamma_a1=EvtPDL::getWidth(A1M);
 133   //  double bwm_a2=EvtPDL::getMeanMass(A2M);
 134   //  double gamma_a2=EvtPDL::getWidth(A2M);
 135   double bwm_rho=EvtPDL::getMeanMass(RHO0);
 136   double gamma_rho=EvtPDL::getWidth(RHO0);
 137
 138   amp_a1=(sqrt(gamma_a1/EvtConst::twoPi)/
 139     ((p4a1).mass()-bwm_a1-EvtComplex(0.0,0.5*gamma_a1)))*
 140          (sqrt(gamma_rho/EvtConst::twoPi)/
 141     ((p4rho).mass()-bwm_rho-EvtComplex(0.0,0.5*gamma_rho)));
 142
 143   return amp_a1*
 144     (vb.get(1)*vpi.get(1)+vb.get(2)*vpi.get(2)+vb.get(3)*vpi.get(3));
 145
 146 }
 147
 148
 149 std::string EvtBTo4piCP::getName(){
 150
 151   return "BTO4PI_CP";
 152
 153 }
 154
 155
 156 EvtDecayBase* EvtBTo4piCP::clone(){
 157
 158   return new EvtBTo4piCP;
 159
 160 }
 161
 162 void EvtBTo4piCP::init(){
 163
 164   // check that there are 18 arguments
 165   checkNArg(18);
 166   checkNDaug(4);
 167
 168   checkSpinParent(EvtSpinType::SCALAR);
 169
 170   checkSpinDaughter(0,EvtSpinType::SCALAR);
 171   checkSpinDaughter(1,EvtSpinType::SCALAR);
 172   checkSpinDaughter(2,EvtSpinType::SCALAR);
 173   checkSpinDaughter(3,EvtSpinType::SCALAR);
 174 }
 175
 176 void EvtBTo4piCP::decay( EvtParticle *p){
 177
 178   //added by Lange Jan4,2000
 179   static EvtId B0=EvtPDL::getId("B0");
 180   static EvtId B0B=EvtPDL::getId("anti-B0");
 181
 182
 183   double t;
 184   EvtId other_b;
 185
 186   EvtCPUtil::OtherB(p,t,other_b);
 187
 188   p->initializePhaseSpace(getNDaug(),getDaugs());
 189   EvtVector4R mom1 = p->getDaug(0)->getP4();
 190   EvtVector4R mom2 = p->getDaug(1)->getP4();
 191   EvtVector4R mom3 = p->getDaug(2)->getP4();
 192   EvtVector4R mom4 = p->getDaug(3)->getP4();
 193
 194   //  double alpha=getArg(0);
 195   //double dm=getArg(1);
 196
 197    EvtComplex amp;
 198
 199    EvtComplex A,Abar;
 200
 201
 202    EvtComplex A_a1p,Abar_a1p,A_a2p,Abar_a2p;
 203    EvtComplex A_a1m,Abar_a1m,A_a2m,Abar_a2m;
 204
 205    A_a1p=EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
 206    Abar_a1p=EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
 207
 208    A_a2p=EvtComplex(getArg(6)*cos(getArg(7)),getArg(6)*sin(getArg(7)));
 209    Abar_a2p=EvtComplex(getArg(8)*cos(getArg(9)),getArg(8)*sin(getArg(9)));
 210
 211    A_a1m=EvtComplex(getArg(10)*cos(getArg(11)),getArg(10)*sin(getArg(11)));
 212    Abar_a1m=EvtComplex(getArg(12)*cos(getArg(13)),getArg(12)*sin(getArg(13)));
 213
 214    A_a2m=EvtComplex(getArg(14)*cos(getArg(15)),getArg(14)*sin(getArg(15)));
 215    Abar_a2m=EvtComplex(getArg(16)*cos(getArg(17)),getArg(16)*sin(getArg(17)));
 216
 217    EvtComplex a2p_amp=EvtAmpA2(mom1,mom2,mom3,mom4)+
 218                       EvtAmpA2(mom1,mom4,mom3,mom2)+
 219                       EvtAmpA2(mom3,mom2,mom1,mom4)+
 220                       EvtAmpA2(mom3,mom4,mom1,mom2);
 221
 222    EvtComplex a2m_amp=EvtAmpA2(mom2,mom3,mom4,mom1)+
 223                       EvtAmpA2(mom2,mom1,mom4,mom3)+
 224                       EvtAmpA2(mom4,mom3,mom2,mom1)+
 225                       EvtAmpA2(mom4,mom1,mom2,mom3);
 226
 227    EvtComplex a1p_amp=EvtAmpA1(mom1,mom2,mom3,mom4)+
 228                       EvtAmpA1(mom1,mom4,mom3,mom2)+
 229                       EvtAmpA1(mom3,mom2,mom1,mom4)+
 230                       EvtAmpA1(mom3,mom4,mom1,mom2);
 231
 232    EvtComplex a1m_amp=EvtAmpA1(mom2,mom3,mom4,mom1)+
 233                       EvtAmpA1(mom2,mom1,mom4,mom3)+
 234                       EvtAmpA1(mom4,mom3,mom2,mom1)+
 235                       EvtAmpA1(mom4,mom1,mom2,mom3);
 236
 237
 238    A=A_a2p*a2p_amp+A_a1p*a1p_amp+
 239      A_a2m*a2m_amp+A_a1m*a1m_amp;
 240    Abar=Abar_a2p*a2p_amp+Abar_a1p*a1p_amp+
 241         Abar_a2m*a2m_amp+Abar_a1m*a1m_amp;
 242
 243
 244    if (other_b==B0B){
 245      amp=A*cos(getArg(1)*t/(2*EvtConst::c))+
 246        EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
 247        getArg(2)*EvtComplex(0.0,1.0)*Abar*sin(getArg(1)*t/(2*EvtConst::c));
 248    }
 249    if (other_b==B0){
 250      amp=A*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
 251        EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
 252        getArg(2)*Abar*cos(getArg(1)*t/(2*EvtConst::c));
 253    }
 254
 255    vertex(amp);
 256
 257    return ;
 258 }
 259