TEvtGen/EvtGen/EvtGenModels/EvtBto2piCPiso.cpp

   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: EvtBto2piCPiso.cc
  12 //
  13 // Description: Routine to decay B -> pi pi with isospin amplitudes
  14 //
  15 // Modification history:
  16 //
  17 //    RYD,NK     Febuary 7, 1998         Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include <stdlib.h>
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtRandom.hh"
  25 #include "EvtGenBase/EvtGenKine.hh"
  26 #include "EvtGenBase/EvtCPUtil.hh"
  27 #include "EvtGenBase/EvtPDL.hh"
  28 #include "EvtGenModels/EvtBto2piCPiso.hh"
  29 #include "EvtGenBase/EvtReport.hh"
  30 #include "EvtGenBase/EvtId.hh"
  31 #include <string>
  32 #include "EvtGenBase/EvtConst.hh"
  33
  34 EvtBto2piCPiso::~EvtBto2piCPiso() {}
  35
  36 std::string EvtBto2piCPiso::getName(){
  37
  38   return "BTO2PI_CP_ISO";
  39
  40 }
  41
  42
  43 EvtDecayBase* EvtBto2piCPiso::clone(){
  44
  45   return new EvtBto2piCPiso;
  46
  47 }
  48
  49 void EvtBto2piCPiso::init(){
  50
  51   // check that there are 10 arguments
  52
  53   checkNArg(10);
  54   checkNDaug(2);
  55
  56   checkSpinParent(EvtSpinType::SCALAR);
  57
  58   checkSpinDaughter(0,EvtSpinType::SCALAR);
  59   checkSpinDaughter(1,EvtSpinType::SCALAR);
  60
  61 }
  62
  63
  64 void EvtBto2piCPiso::initProbMax() {
  65
  66   //added by Lange Jan4,2000
  67   static EvtId PI0=EvtPDL::getId("pi0");
  68   static EvtId PIP=EvtPDL::getId("pi+");
  69   static EvtId PIM=EvtPDL::getId("pi-");
  70
  71 //this may need to be revised
  72
  73 if (((getDaugs()[0]==PIP) && (getDaugs()[1]==PIM)) || ((getDaugs()[0]==PIM) && (getDaugs()[1]==PIP))) {
  74    setProbMax(4.0*(getArg(2)*getArg(2)+getArg(4)*getArg(4)));
  75 }
  76
  77 if ((getDaugs()[0]==PI0) && (getDaugs()[1]==PI0)) {
  78    setProbMax(2.0*(4.0*getArg(2)*getArg(2)+getArg(4)*getArg(4)));
  79 }
  80
  81 if (((getDaugs()[0]==PIP) && (getDaugs()[1]==PI0)) || ((getDaugs()[0]==PI0) && (getDaugs()[1]==PIP))) {
  82    setProbMax(6.0*getArg(2)*getArg(2));
  83 }
  84
  85 if (((getDaugs()[0]==PI0) && (getDaugs()[1]==PIM)) || ((getDaugs()[0]==PIM) && (getDaugs()[1]==PI0))) {
  86    setProbMax(6.0*getArg(4)*getArg(4));
  87 }
  88
  89 }
  90
  91 void EvtBto2piCPiso::decay( EvtParticle *p ){
  92
  93   p->initializePhaseSpace(getNDaug(),getDaugs());
  94
  95   //added by Lange Jan4,2000
  96   static EvtId B0=EvtPDL::getId("B0");
  97   static EvtId B0B=EvtPDL::getId("anti-B0");
  98   static EvtId PI0=EvtPDL::getId("pi0");
  99   static EvtId PIP=EvtPDL::getId("pi+");
 100   static EvtId PIM=EvtPDL::getId("pi-");
 101
 102   double t;
 103   EvtId other_b;
 104   int charged=0;
 105
 106 //randomly generate the tag (B0 or B0B)
 107
 108   double tag = EvtRandom::Flat(0.0,1.0);
 109   if (tag < 0.5) {
 110
 111     EvtCPUtil::getInstance()->OtherB(p,t,other_b,1.0);
 112     other_b = B0;
 113   }
 114   else {
 115
 116     EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.0);
 117     other_b = B0B;
 118   }
 119
 120   EvtComplex amp;
 121
 122   EvtComplex A,Abar;
 123   EvtComplex A2, A2_bar, A0, A0_bar;
 124
 125   A2 = EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
 126   A2_bar = EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
 127
 128   A0 = EvtComplex(getArg(6)*cos(getArg(7)),getArg(6)*sin(getArg(7)));
 129   A0_bar = EvtComplex(getArg(8)*cos(getArg(9)),getArg(8)*sin(getArg(9)));
 130
 131   //depending on what combination of pi pi we have, there will be different
 132   //A and Abar
 133
 134   if (((getDaugs()[0]==PIP) && (getDaugs()[1]==PI0)) || ((getDaugs()[0]==PI0) && (getDaugs()[1]==PIP))) {
 135
 136     //pi+ pi0, so just A_2
 137
 138     charged = 1;
 139     A = 3.0*A2;
 140
 141   }
 142
 143   if (((getDaugs()[0]==PI0) && (getDaugs()[1]==PIM)) || ((getDaugs()[0]==PIM) && (getDaugs()[1]==PI0))) {
 144
 145     //pi- pi0, so just A2_bar
 146
 147     charged = 1;
 148     A = 3.0*A2_bar;
 149
 150   }
 151
 152   if (((getDaugs()[0]==PIP) && (getDaugs()[1]==PIM)) || ((getDaugs()[0]==PIM) && (getDaugs()[1]==PIP))) {
 153
 154 //pi+ pi-, so A_2 - A_0
 155
 156     charged = 0;
 157     A=sqrt(2.0)*(A2-A0);
 158     Abar=sqrt(2.0)*(A2_bar-A0_bar);
 159   }
 160
 161   if ((getDaugs()[0]==PI0) && (getDaugs()[1]==PI0)) {
 162
 163     //pi0 pi0, so 2*A_2 + A_0
 164
 165     charged = 0;
 166     A=2.0*A2 + A0;
 167     Abar=2.0*A2_bar + A0_bar;
 168   }
 169
 170   if(charged == 0) {
 171
 172     if (other_b==B0B){
 173       amp=A*cos(getArg(1)*t/(2*EvtConst::c))+
 174         EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
 175         EvtComplex(0.0,1.0)*Abar*sin(getArg(1)*t/(2*EvtConst::c));
 176     }
 177     if (other_b==B0){
 178       amp=A*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
 179         EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+
 180         Abar*cos(getArg(1)*t/(2*EvtConst::c));
 181     }
 182   }
 183   else amp = A;
 184
 185   vertex(amp);
 186
 187   return ;
 188 }
 189
 190 std::string EvtBto2piCPiso::getParamName(int i) {
 191   switch(i) {
 192   case 0:
 193     return "weakPhase";
 194   case 1:
 195     return "deltaM";
 196   case 2:
 197     return "A2";
 198   case 3:
 199     return "A2Phase";
 200   case 4:
 201     return "A2bar";
 202   case 5:
 203     return "A2barPhase";
 204   case 6:
 205     return "A0";
 206   case 7:
 207     return "A0Phase";
 208   case 8:
 209     return "A0bar";
 210   case 9:
 211     return "A0barPhase";
 212   default:
 213     return "";
 214   }
 215 }
 216
 217 std::string EvtBto2piCPiso::getParamDefault(int i) {
 218   switch(i) {
 219   case 2:
 220     return "1.0";
 221   case 3:
 222     return "0.0";
 223   case 4:
 224     return "1.0";
 225   case 5:
 226     return "0.0";
 227   case 6:
 228     return "1.0";
 229   case 7:
 230     return "0.0";
 231   case 8:
 232     return "1.0";
 233   case 9:
 234     return "0.0";
 235   default:
 236     return "";
 237   }
 238 }