TEvtGen/EvtGenModels/EvtSVSCPiso.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: EvtSVSCPiso.cc
  12 //
  13 // Description: Routine to decay scalar -> vectors scalar
  14 //              with CP violation and isospin amplitudes.
  15 //              More specifically, it is indended to handle
  16 //              decays like B->rho pi and B->a1 pi.
  17 //
  18 // Modification history:
  19 //
  20 //    RYD/NK    Febuary 16, 1998          Module created
  21 //
  22 //------------------------------------------------------------------------
  23 //
  24 #include "EvtGenBase/EvtPatches.hh"
  25 #include <stdlib.h>
  26 #include "EvtGenBase/EvtParticle.hh"
  27 #include "EvtGenBase/EvtRandom.hh"
  28 #include "EvtGenBase/EvtGenKine.hh"
  29 #include "EvtGenBase/EvtCPUtil.hh"
  30 #include "EvtGenBase/EvtPDL.hh"
  31 #include "EvtGenBase/EvtReport.hh"
  32 #include "EvtGenBase/EvtVector4C.hh"
  33 #include "EvtGenModels/EvtSVSCPiso.hh"
  34 #include "EvtGenBase/EvtId.hh"
  35 #include <string>
  36 #include "EvtGenBase/EvtConst.hh"
  37
  38 EvtSVSCPiso::~EvtSVSCPiso() {}
  39
  40 std::string EvtSVSCPiso::getName(){
  41
  42   return "SVS_CP_ISO";
  43
  44 }
  45
  46
  47 EvtDecayBase* EvtSVSCPiso::clone(){
  48
  49   return new EvtSVSCPiso;
  50
  51 }
  52
  53 void EvtSVSCPiso::init(){
  54
  55   // check that there are 26 arguments
  56   checkNArg(26);
  57   checkNDaug(2);
  58
  59   checkSpinParent(EvtSpinType::SCALAR);
  60
  61   checkSpinDaughter(0,EvtSpinType::VECTOR);
  62   checkSpinDaughter(1,EvtSpinType::SCALAR);
  63
  64 }
  65
  66
  67 void EvtSVSCPiso::initProbMax(){
  68
  69 //this might need some revision..
  70
  71 if ((EvtPDL::chg3(getDaug(0)) > 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
  72   setProbMax(2.0*(getArg(3)*getArg(3) + 4.0*getArg(23)*getArg(23)));
  73 }
  74
  75 if ((EvtPDL::chg3(getDaug(0)) < 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
  76   setProbMax(2.0*(getArg(5)*getArg(5) + 4.0*getArg(25)*getArg(25)));
  77 }
  78
  79 if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) > 0)) {
  80   setProbMax(2.0*(getArg(7)*getArg(7) + 4.0*getArg(23)*getArg(23)));
  81 }
  82
  83 if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) < 0)) {
  84   setProbMax(2.0*(getArg(9)*getArg(9) + 4.0*getArg(25)*getArg(25)));
  85 }
  86
  87 if ((EvtPDL::chg3(getDaug(0)) > 0) && (EvtPDL::chg3(getDaug(1)) < 0)) {
  88   setProbMax(2.0*(getArg(11)*getArg(11) + getArg(23)*getArg(23) +
  89                   getArg(19)*getArg(19) + getArg(13)*getArg(13) +
  90                   getArg(25)*getArg(25) + getArg(21)*getArg(21)));
  91 }
  92
  93 if ((EvtPDL::chg3(getDaug(0)) < 0) && (EvtPDL::chg3(getDaug(1)) > 0)) {
  94   setProbMax(2.0*(getArg(15)*getArg(15) + getArg(23)*getArg(23) +
  95                   getArg(19)*getArg(19) + getArg(17)*getArg(17) +
  96                   getArg(25)*getArg(25) + getArg(21)*getArg(21)));
  97 }
  98
  99 if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
 100    setProbMax(2.0*(getArg(7)*getArg(7) + getArg(3)*getArg(3) + getArg(11)*getArg(11) +
 101                   getArg(15)*getArg(15) + 4.0*getArg(19)*getArg(19) + getArg(9)*getArg(9)+
 102                    getArg(5)*getArg(5) + getArg(13)*getArg(13) + getArg(17)*getArg(17) +
 103                    4.0*getArg(21)*getArg(21)));
 104 }
 105
 106 }
 107
 108
 109 void EvtSVSCPiso::decay( EvtParticle *p){
 110
 111   //added by Lange Jan4,2000
 112   static EvtId B0=EvtPDL::getId("B0");
 113   static EvtId B0B=EvtPDL::getId("anti-B0");
 114
 115   double t;
 116   EvtId other_b;
 117   int charged(0);
 118
 119   int first_time=0;
 120   int flip=0;
 121   EvtId ds[2];
 122
 123
 124 //randomly generate the tag (B0 or B0B)
 125
 126    double tag = EvtRandom::Flat(0.0,1.0);
 127    if (tag < 0.5) {
 128
 129     EvtCPUtil::OtherB(p,t,other_b,1.0);
 130     other_b = B0;
 131    }
 132    else {
 133
 134     EvtCPUtil::OtherB(p,t,other_b,0.0);
 135     other_b = B0B;
 136    }
 137
 138   if (p->getNDaug()==0) first_time=1;
 139
 140   if (first_time){
 141     if (EvtRandom::Flat(0.0,1.0)<getArg(3)) flip=1;
 142   }
 143   else{
 144     if (getDaug(0)!=p->getDaug(0)->getId()) flip=1;
 145   }
 146
 147   if (!flip) {
 148     ds[0]=getDaug(0);
 149     ds[1]=getDaug(1);
 150   }
 151   else{
 152     ds[0]=EvtPDL::chargeConj(getDaug(0));
 153     ds[1]=EvtPDL::chargeConj(getDaug(1));
 154   }
 155
 156   p->initializePhaseSpace(getNDaug(),ds);
 157
 158   EvtParticle *v,*s;
 159   v=p->getDaug(0);
 160   s=p->getDaug(1);
 161
 162    EvtComplex amp;
 163
 164    EvtComplex A_f,Abar_f;
 165    EvtComplex A_fbar,Abar_fbar;
 166    EvtComplex Apm, Apm_bar, Amp, Amp_bar;
 167
 168    EvtComplex Tp0, Tp0_bar, T0p, T0p_bar,Tpm, Tpm_bar, Tmp, Tmp_bar;
 169    EvtComplex P1, P1_bar, P0, P0_bar;
 170
 171    Tp0 = EvtComplex(getArg(3)*cos(getArg(4)),getArg(3)*sin(getArg(4)));
 172    Tp0_bar = EvtComplex(getArg(5)*cos(getArg(6)),getArg(5)*sin(getArg(6)));
 173    T0p = EvtComplex(getArg(7)*cos(getArg(8)),getArg(7)*sin(getArg(8)));
 174    T0p_bar = EvtComplex(getArg(9)*cos(getArg(10)),getArg(9)*sin(getArg(10)));
 175    Tpm = EvtComplex(getArg(11)*cos(getArg(12)),getArg(11)*sin(getArg(12)));
 176    Tpm_bar = EvtComplex(getArg(13)*cos(getArg(14)),getArg(13)*sin(getArg(14)));
 177    Tmp = EvtComplex(getArg(15)*cos(getArg(16)),getArg(15)*sin(getArg(16)));
 178    Tmp_bar = EvtComplex(getArg(17)*cos(getArg(18)),getArg(17)*sin(getArg(18)));
 179    P0 = EvtComplex(getArg(19)*cos(getArg(20)),getArg(19)*sin(getArg(20)));
 180    P0_bar = EvtComplex(getArg(21)*cos(getArg(22)),getArg(21)*sin(getArg(22)));
 181    P1 = EvtComplex(getArg(23)*cos(getArg(24)),getArg(23)*sin(getArg(24)));
 182    P1_bar = EvtComplex(getArg(25)*cos(getArg(26)),getArg(25)*sin(getArg(26)));
 183
 184
 185 //***********************charged modes****************************
 186
 187  if ((EvtPDL::chg3(getDaug(0)) > 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
 188
 189 //V+ S0, so T+0 + 2 P1
 190
 191    charged = 1;
 192    A_f = Tp0 + 2.0*P1;
 193  }
 194
 195  if ((EvtPDL::chg3(getDaug(0)) < 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
 196
 197 //V- S0, so T+0_bar + 2P1_bar
 198
 199    charged = 1;
 200    A_f = Tp0_bar + 2.0*P1_bar;
 201  }
 202
 203  if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) > 0)) {
 204
 205 //V0 S+, so T0+ - 2 P1
 206
 207    charged = 1;
 208    A_f = T0p - 2.0*P1;
 209  }
 210
 211  if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) < 0)) {
 212
 213 //V0 S-, so T0+_bar - 2 P1_bar
 214
 215    charged = 1;
 216    A_f = T0p_bar - 2.0*P1_bar;
 217  }
 218
 219
 220 //***********************neutral modes***************************
 221
 222
 223 //V+ S-, so Af = T+- + P1 + P0
 224 Apm = Tpm + P1 + P0;
 225 Apm_bar = Tpm_bar + P1_bar + P0_bar;
 226
 227 //V- S+, so Af = T-+ - P1 + P0
 228 Amp = Tmp - P1 + P0;
 229 Amp_bar = Tmp_bar - P1_bar + P0;
 230
 231
 232  if ((EvtPDL::chg3(getDaug(0)) > 0 ) && (EvtPDL::chg3(getDaug(1)) < 0)) {
 233
 234 //V+ S-
 235    charged = 0;
 236    A_f = Apm;
 237    Abar_f = Apm_bar;
 238    A_fbar = Amp;
 239    Abar_fbar = Amp_bar;
 240
 241  }
 242
 243  if ((EvtPDL::chg3(getDaug(0)) < 0 ) && (EvtPDL::chg3(getDaug(1)) > 0)) {
 244
 245 //V- S+
 246    charged = 0;
 247    A_f = Amp;
 248    Abar_f = Amp_bar;
 249    A_fbar = Apm;
 250    Abar_fbar = Apm_bar;
 251
 252  }
 253
 254  if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
 255
 256 //V0 S0
 257    charged = 0;
 258    A_f = T0p + Tp0 - Tpm - Tmp - 2.0*P0 ;
 259    Abar_f = T0p_bar + Tp0_bar - Tpm_bar - Tmp_bar - 2.0*P0_bar;
 260    A_fbar = A_f;
 261    Abar_fbar = Abar_f;
 262
 263  }
 264
 265 if (charged==0) {
 266
 267    if (!flip) {
 268      if (other_b==B0B){
 269
 270        amp=A_f*cos(getArg(1)*t/(2*EvtConst::c))+
 271          EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
 272          EvtComplex(0.0,1.0)*Abar_f*sin(getArg(1)*t/(2*EvtConst::c));
 273      }
 274      if (other_b==B0){
 275
 276        amp=A_f*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
 277          EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
 278          Abar_f*cos(getArg(1)*t/(2*EvtConst::c));
 279      }
 280    }
 281    else{
 282      if (other_b==B0B){
 283
 284        amp=A_fbar*cos(getArg(1)*t/(2*EvtConst::c))+
 285          EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
 286          EvtComplex(0.0,1.0)*Abar_fbar*sin(getArg(1)*t/(2*EvtConst::c));
 287      }
 288      if (other_b==B0){
 289
 290        amp=A_fbar*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
 291          EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
 292          Abar_fbar*cos(getArg(1)*t/(2*EvtConst::c));
 293      }
 294    }
 295
 296 }
 297 else amp = A_f;
 298
 299   EvtVector4R p4_parent;
 300
 301   p4_parent=v->getP4()+s->getP4();
 302
 303   double norm=1.0/v->getP4().d3mag();
 304
 305    vertex(0,amp*norm*p4_parent*(v->epsParent(0)));
 306    vertex(1,amp*norm*p4_parent*(v->epsParent(1)));
 307    vertex(2,amp*norm*p4_parent*(v->epsParent(2)));
 308
 309   return ;
 310 }
 311