TEvtGen/EvtGenModels/EvtDDalitz.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: EvtDDalitz.cc
  12 //
  13 // Description: Routine to handle three-body decays of D0/D0_bar or D+/D-
  14 //
  15 // Modification history:
  16 //
  17 //    NK     September 3, 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/EvtPDL.hh"
  26 #include "EvtGenBase/EvtReport.hh"
  27 #include "EvtGenBase/EvtResonance.hh"
  28 #include "EvtGenBase/EvtResonance2.hh"
  29 #include "EvtGenModels/EvtDDalitz.hh"
  30 #include <string>
  31 #include "EvtGenBase/EvtConst.hh"
  32 #include "EvtGenBase/EvtFlatte.hh"
  33 #include "EvtGenBase/EvtDecayTable.hh"
  34 using std::endl;
  35
  36 EvtDDalitz::~EvtDDalitz() {}
  37
  38 std::string EvtDDalitz::getName(){
  39
  40   return "D_DALITZ";
  41
  42 }
  43
  44
  45 EvtDecayBase* EvtDDalitz::clone(){
  46
  47   return new EvtDDalitz;
  48
  49 }
  50
  51 void EvtDDalitz::init(){
  52
  53   // check that there are 0 arguments
  54
  55   static EvtId DM=EvtPDL::getId("D-");
  56   static EvtId DP=EvtPDL::getId("D+");
  57   static EvtId D0=EvtPDL::getId("D0");
  58   static EvtId D0B=EvtPDL::getId("anti-D0");
  59   static EvtId DSP=EvtPDL::getId("D_s+");
  60   static EvtId DSM=EvtPDL::getId("D_s-");
  61   static EvtId KM=EvtPDL::getId("K-");
  62   static EvtId KP=EvtPDL::getId("K+");
  63   static EvtId K0=EvtPDL::getId("K0");
  64   static EvtId KB=EvtPDL::getId("anti-K0");
  65   static EvtId KL=EvtPDL::getId("K_L0");
  66   static EvtId KS=EvtPDL::getId("K_S0");
  67   static EvtId PIM=EvtPDL::getId("pi-");
  68   static EvtId PIP=EvtPDL::getId("pi+");
  69   static EvtId PI0=EvtPDL::getId("pi0");
  70
  71   static double MPI = EvtPDL::getMeanMass(PI0);
  72   static double MKP = EvtPDL::getMeanMass(KP);
  73
  74
  75   checkNArg(0);
  76   checkNDaug(3);
  77
  78   checkSpinParent(EvtSpinType::SCALAR);
  79
  80   checkSpinDaughter(0,EvtSpinType::SCALAR);
  81   checkSpinDaughter(1,EvtSpinType::SCALAR);
  82   checkSpinDaughter(2,EvtSpinType::SCALAR);
  83
  84   EvtId parnum=getParentId();
  85   EvtId d1=getDaug(0);
  86   EvtId d2=getDaug(1);
  87   EvtId d3=getDaug(2);
  88   _flag=0;
  89   if ( parnum == D0 ) {
  90     //look for either a K- pi+ pi0 or K0bar pi+ pi-
  91     if ( d1==KM && d2==PIP && d3==PI0 ) { _flag=4; _d1=0; _d2=1; _d3=2;}
  92     if ( d1==KM && d3==PIP && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;}
  93     if ( d2==KM && d1==PIP && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;}
  94     if ( d2==KM && d3==PIP && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;}
  95     if ( d3==KM && d1==PIP && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;}
  96     if ( d3==KM && d2==PIP && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;}
  97
  98     if ( d1==KB && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=2; _d3=1;}
  99     if ( d1==KB && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
 100     if ( d2==KB && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 101     if ( d2==KB && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 102     if ( d3==KB && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 103     if ( d3==KB && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 104
 105     if ( d1==KL && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=2; _d3=1;}
 106     if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
 107     if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 108     if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 109     if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 110     if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 111
 112     if ( d1==KS && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=2; _d3=1;}
 113     if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
 114     if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 115     if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 116     if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 117     if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 118
 119
 120     if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 121     if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 122     if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 123     if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 124     if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 125     if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 126
 127     if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 128     if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 129     if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 130     if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 131     if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 132     if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 133
 134     if ( d1==K0 && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 135     if ( d1==K0 && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 136     if ( d2==K0 && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 137     if ( d2==K0 && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 138     if ( d3==K0 && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 139     if ( d3==K0 && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 140   }
 141   if ( parnum == D0B ) {
 142     //look for either a K+ pi- pi0 or K0 pi+ pi-
 143     if ( d1==KP && d2==PIM && d3==PI0 )  { _flag=4; _d1=0; _d2=1; _d3=2;}
 144     if ( d1==KP && d3==PIM && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;}
 145     if ( d2==KP && d1==PIM && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;}
 146     if ( d2==KP && d3==PIM && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;}
 147     if ( d3==KP && d1==PIM && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;}
 148     if ( d3==KP && d2==PIM && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;}
 149
 150     if ( d1==K0 && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=1; _d3=2;}
 151     if ( d1==K0 && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
 152     if ( d2==K0 && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 153     if ( d2==K0 && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 154     if ( d3==K0 && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 155     if ( d3==K0 && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 156
 157     if ( d1==KL && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=1; _d3=2;}
 158     if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
 159     if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 160     if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 161     if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 162     if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 163
 164     if ( d1==KS && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=1; _d3=2;}
 165     if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
 166     if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 167     if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 168     if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 169     if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 170
 171     if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 172     if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 173     if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 174     if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 175     if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 176     if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 177
 178     if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 179     if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 180     if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 181     if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 182     if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 183     if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 184
 185     if ( d1==K0 && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 186     if ( d1==K0 && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 187     if ( d2==K0 && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 188     if ( d2==K0 && d3==KP && d1==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 189     if ( d3==K0 && d1==KP && d2==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 190     if ( d3==K0 && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 191
 192   }
 193
 194   if ( parnum == DP ) {
 195     //look for K- pi+ pi+
 196     if ( d1==KB && d2==PIP && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 197     if ( d1==KB && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 198     if ( d2==KB && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 199     if ( d2==KB && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 200     if ( d3==KB && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 201     if ( d3==KB && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 202
 203     if ( d1==KL && d2==PIP && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 204     if ( d1==KL && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 205     if ( d2==KL && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 206     if ( d2==KL && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 207     if ( d3==KL && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 208     if ( d3==KL && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 209
 210     if ( d1==KS && d2==PIP && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 211     if ( d1==KS && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 212     if ( d2==KS && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 213     if ( d2==KS && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 214     if ( d3==KS && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 215     if ( d3==KS && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 216
 217     if ( d1==KM && d2==PIP && d3==PIP )  { _flag=1; _d1=0; _d2=1; _d3=2;}
 218     if ( d2==KM && d1==PIP && d3==PIP ) { _flag=1; _d1=1; _d2=0; _d3=2;}
 219     if ( d3==KM && d1==PIP && d2==PIP ) { _flag=1; _d1=2; _d2=0; _d3=1;}
 220   }
 221
 222   if ( parnum == DM ) {
 223     //look for K- pi+ pi+
 224     if ( d1==K0 && d2==PIM && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 225     if ( d1==K0 && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 226     if ( d2==K0 && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 227     if ( d2==K0 && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 228     if ( d3==K0 && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 229     if ( d3==K0 && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 230
 231     if ( d1==KL && d2==PIM && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 232     if ( d1==KL && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 233     if ( d2==KL && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 234     if ( d2==KL && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 235     if ( d3==KL && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 236     if ( d3==KL && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 237
 238     if ( d1==KS && d2==PIM && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 239     if ( d1==KS && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 240     if ( d2==KS && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 241     if ( d2==KS && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 242     if ( d3==KS && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 243     if ( d3==KS && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 244
 245     if ( d1==KP && d2==PIM && d3==PIM )  { _flag=1; _d1=0; _d2=1; _d3=2;}
 246     if ( d2==KP && d1==PIM && d3==PIM ) { _flag=1; _d1=1; _d2=0; _d3=2;}
 247     if ( d3==KP && d1==PIM && d2==PIM ) { _flag=1; _d1=2; _d2=0; _d3=1;}
 248   }
 249
 250   if ( parnum == DSP ) {
 251      if ( d1==KM && d2==KP && d3==PIP ) { _flag=6; _d1=0; _d2=1; _d3=2; }
 252      if ( d1==KM && d3==KP && d2==PIP ) { _flag=6; _d1=0; _d2=2; _d3=1; }
 253      if ( d2==KM && d1==KP && d3==PIP ) { _flag=6; _d1=1; _d2=0; _d3=2; }
 254      if ( d2==KM && d3==KP && d1==PIP ) { _flag=6; _d1=1; _d2=2; _d3=0; }
 255      if ( d3==KM && d1==KP && d2==PIP ) { _flag=6; _d1=2; _d2=0; _d3=1; }
 256      if ( d3==KM && d2==KP && d1==PIP ) { _flag=6; _d1=2; _d2=1; _d3=0; }
 257   }
 258
 259   if ( parnum == DSM ) {
 260      if ( d1==KP && d2==KM && d3==PIM ) { _flag=6; _d1=0; _d2=1; _d3=2; }
 261      if ( d1==KP && d3==KM && d2==PIM ) { _flag=6; _d1=0; _d2=2; _d3=1; }
 262      if ( d2==KP && d1==KM && d3==PIM ) { _flag=6; _d1=1; _d2=0; _d3=2; }
 263      if ( d2==KP && d3==KM && d1==PIM ) { _flag=6; _d1=1; _d2=2; _d3=0; }
 264      if ( d3==KP && d1==KM && d2==PIM ) { _flag=6; _d1=2; _d2=0; _d3=1; }
 265      if ( d3==KP && d2==KM && d1==PIM ) { _flag=6; _d1=2; _d2=1; _d3=0; }
 266   }
 267
 268   if ( _flag==6) {
 269     _kkpi_params.push_back(EvtFlatteParam(MPI, MPI, 0.406));
 270     _kkpi_params.push_back(EvtFlatteParam(MKP, MKP, 0.800));
 271   }
 272
 273   if ( _flag==0) {
 274     report(ERROR,"EvtGen") << "EvtDDaltiz: Invalid mode."<<endl;
 275     assert(0);
 276   }
 277 }
 278
 279 void EvtDDalitz::initProbMax() {
 280
 281 //probmax different for different modes!
 282
 283   if ( _flag==1 ) {setProbMax(9.6);}
 284   if ( _flag==2 ) {setProbMax(147.9);}
 285   if ( _flag==3 ) {setProbMax(5000.0);}
 286   if ( _flag==4 ) {setProbMax(3000.0);}
 287   if ( _flag==5 ) {setProbMax(10000000.0);}
 288   if ( _flag==6 ) {setProbMax(50000.0);}
 289
 290 }
 291
 292 void EvtDDalitz::decay( EvtParticle *p){
 293
 294   static EvtId BP = EvtPDL::getId("B+");
 295   static EvtId BM = EvtPDL::getId("B-");
 296   static EvtId B0 = EvtPDL::getId("B0");
 297   static EvtId B0B = EvtPDL::getId("anti-B0");
 298   static EvtId DM=EvtPDL::getId("D-");
 299   static EvtId DP=EvtPDL::getId("D+");
 300   static EvtId D0=EvtPDL::getId("D0");
 301   static EvtId D0B=EvtPDL::getId("anti-D0");
 302   static EvtId KM=EvtPDL::getId("K-");
 303   static EvtId KP=EvtPDL::getId("K+");
 304   static EvtId K0=EvtPDL::getId("K0");
 305   static EvtId KB=EvtPDL::getId("anti-K0");
 306   static EvtId PIM=EvtPDL::getId("pi-");
 307   static EvtId PIP=EvtPDL::getId("pi+");
 308   static EvtId PI0=EvtPDL::getId("pi0");
 309
 310   double oneby2 = 0.707106782;
 311
 312   bool isBToDK=false;
 313   if ( p -> getParent () ) {
 314     EvtId parId = p -> getParent()->getId ();
 315     if ( ( BP == parId ) || ( BM == parId ) || ( B0 == parId ) ||
 316                        ( B0B == parId ) )
 317       if (EvtDecayTable::getDecayFunc(p->getParent())->getName() == "BTODDALITZCPK") isBToDK=true;
 318   }
 319
 320
 321 //same structure for all of these decays
 322
 323   p->initializePhaseSpace(getNDaug(),getDaugs());
 324   EvtVector4R moms1 = p->getDaug(_d1)->getP4();
 325   EvtVector4R moms2 = p->getDaug(_d2)->getP4();
 326   EvtVector4R moms3 = p->getDaug(_d3)->getP4();
 327
 328   EvtVector4R p4_p;
 329   p4_p.set(p->mass(),0.0,0.0,0.0);
 330
 331   EvtComplex amp(1.0,0.0);
 332
 333 //now determine which D and which decay
 334
 335 //data from Anjos et al, Phys.Rev.D 1993, v.48,num.1,p.56 (E691 resuls)
 336 //for D+ -> K- pi+ pi+, and from Adler et al, Phys.Lett. B196 (1987), 107
 337 //(Mark III results) for D+ -> K0bar pi+ pi0.
 338   //CLEO results for D0->k-pi+pi0
 339
 340   if ( _flag==1) {
 341
 342 //have a D+ -> K- pi+ pi+ decay, or charge conjugate
 343 //Anjos etal e691 - Phys Rev D48, 56 (1993)
 344     EvtResonance DplusRes11(p4_p,moms1,moms2,0.78,-60.0,0.0498,0.89610,1);
 345     EvtResonance DplusRes12(p4_p,moms3,moms1,0.78,-60.0,0.0498,0.89610,1);//K*(892)
 346
 347     EvtResonance DplusRes21(p4_p,moms1,moms2,0.53,132.0,0.287,1.429,0);
 348     EvtResonance DplusRes22(p4_p,moms3,moms1,0.53,132.0,0.287,1.429,0);//K*(1430)
 349
 350     EvtResonance DplusRes31(p4_p,moms1,moms2,0.47,-51.0,0.323,1.714,1);
 351     EvtResonance DplusRes32(p4_p,moms3,moms1,0.47,-51.0,0.323,1.714,1);//K*(1680)
 352
 353     amp = amp + oneby2*(-DplusRes11.resAmpl()+DplusRes12.resAmpl()) + oneby2*(DplusRes21.resAmpl() + DplusRes22.resAmpl()) + oneby2*(-DplusRes31.resAmpl()+ DplusRes32.resAmpl());
 354
 355  }
 356
 357   if ( _flag==2) {
 358
 359 //have a D+ -> K0bar pi+ pi0 decay
 360 //adler etal MarkIII - Phys Lett B196, 107 (1987)
 361 // Results in this paper:
 362 //   Kbar rho+    FitFraction = 68+/-8+/-12    Phase   0
 363 //   Kbar* pi+                  19+/-6+/-6            43+/-23
 364 //   nonres                     13+/-7+/-8           250+/-19
 365 // These numbers below seem not to be exactly the same
 366 // the phases are equiv to -106=254 and 41
 367 //
 368     EvtResonance DplusKpipi0Res1(p4_p,moms2,moms3,1.00,0.00,0.1512,0.7699,1); //rho+
 369     EvtResonance DplusKpipi0Res2(p4_p,moms3,moms1,0.8695,0.7191,0.0498,0.89159,1); //K*0
 370
 371     amp = 0.9522*EvtComplex(cos(-1.8565),sin(-1.8565)) + 1.00*DplusKpipi0Res1.relBrWig(0) + 0.8695*EvtComplex(cos(0.7191),sin(0.7191))*DplusKpipi0Res2.relBrWig(1);
 372
 373   }
 374
 375   if(_flag==3) {
 376     // D0 -> K0 pi+ pi- + CC
 377     // If it does not come from a B->DK, decay it as D0 or D0bar separatly
 378     // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+
 379     // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi-
 380
 381     if ( isBToDK ) {
 382       // Gamma angle in rad.
 383       double gamma = EvtDecayTable::getDecayFunc( p->getParent() )
 384         -> getArg( 0 )  ;
 385       // Strong phase in rad.
 386       double delta =  EvtDecayTable::getDecayFunc( p->getParent() )
 387         -> getArg( 1 )  ;
 388       // Ratio between B->D0K and B->D0barK
 389       double A     =  EvtDecayTable::getDecayFunc( p->getParent() )
 390         -> getArg( 2 )  ;
 391
 392       EvtComplex Factor( fabs( A ) * cos ( delta ) ,
 393                          fabs( A ) * sin ( delta ) ) ;
 394
 395       if ( ( p->getParent()->getId() == BP ) ||
 396            ( p->getParent()->getId() == B0 ) ) {
 397         // the ratio D/Dbar
 398         Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ;
 399         if ( p->getId() == D0 ) {
 400           // the flavor of the particle has no meaning. But we need
 401           // it to know which daughter is pi+ or pi-
 402           // M( B+ or B0 ) = f(Dbar) + factor * f(D)
 403           // f(Dbar) = amplDtoK0PiPi(pD, K0, pi+, pi-)
 404           // f(D)    = amplDtoK0PiPi(pD, K0, pi-, pi+)
 405           // Then ...
 406           amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) +
 407             Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
 408         }
 409         else {
 410           amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) +
 411             Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ;
 412         }
 413       }
 414       else if ( ( p->getParent() -> getId() == BM ) ||
 415                 ( p->getParent() -> getId() == B0B ) ) {
 416         Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ;
 417         // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ...
 418         if ( p->getId() == D0 ) {
 419           amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) +
 420             Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ;
 421         }
 422         else {
 423           amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) +
 424             Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
 425         }
 426       }
 427     }
 428     else {
 429       amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
 430     }
 431   }
 432
 433
 434   if(_flag==4) {
 435
 436     EvtResonance2 DKpipi0Res1(p4_p,moms2,moms3,1.0  ,0.0   ,0.1507,0.770 ,1); //rho
 437     EvtResonance2 DKpipi0Res2(p4_p,moms1,moms2,0.39, -0.2  ,0.0505,0.8961,1); //k*0
 438     EvtResonance2 DKpipi0Res3(p4_p,moms1,moms3,0.44, 163.0 ,0.050 ,0.8915,1); //k*-
 439
 440     EvtResonance2 DKpipi0Res4(p4_p,moms1,moms3,0.77 ,55.5  ,0.294 ,1.412 ,0); //k01430-
 441     EvtResonance2 DKpipi0Res5(p4_p,moms1,moms2,0.85 ,166.0 ,0.294 ,1.412 ,0); //k01430bar
 442     EvtResonance2 DKpipi0Res6(p4_p,moms2,moms3,2.5  ,171.0 ,0.240 ,1.700 ,1); //rho1700
 443     EvtResonance2 DKpipi0Res7(p4_p,moms1,moms3,2.5  ,103.0 ,0.322 ,1.717 ,1); //K*1680-
 444
 445
 446
 447     double pi180inv = 1.0/EvtConst::radToDegrees;
 448
 449     amp = EvtComplex(1.75*cos(31.2*pi180inv),1.75*sin(31.2*pi180inv))
 450       + DKpipi0Res1.resAmpl() + DKpipi0Res2.resAmpl() + DKpipi0Res3.resAmpl()
 451       + DKpipi0Res4.resAmpl() + DKpipi0Res5.resAmpl()
 452       + DKpipi0Res6.resAmpl()
 453       + DKpipi0Res7.resAmpl();
 454
 455   }
 456
 457   if(_flag==5) {
 458
 459     // D0 -> K0 K+ K- + CC
 460     // If it does not come from a B->DK, decay it as D0 or D0bar separatly
 461     // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+
 462     // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi-
 463
 464     if ( isBToDK ){
 465       // Gamma angle in rad.
 466       double gamma = EvtDecayTable::getDecayFunc( p->getParent() )
 467         -> getArg( 0 )  ;
 468       // Strong phase in rad.
 469       double delta =  EvtDecayTable::getDecayFunc( p->getParent() )
 470         -> getArg( 1 )  ;
 471       // Ratio between B->D0K and B->D0barK
 472       double A     =  EvtDecayTable::getDecayFunc( p->getParent() )
 473         -> getArg( 2 )  ;
 474
 475       EvtComplex Factor( fabs( A ) * cos ( delta ) ,
 476                          fabs( A ) * sin ( delta ) ) ;
 477
 478       if ( ( p->getParent()->getId() == BP ) ||
 479            ( p->getParent()->getId() == B0 ) ) {
 480         // the ratio D/Dbar
 481         Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ;
 482         if ( p->getId() == D0 ) {
 483           // the flavor of the particle has no meaning. But we need
 484           // it to know which daughter is pi+ or pi-
 485           // M( B+ or B0 ) = f(Dbar) + factor * f(D)
 486           // f(Dbar) = amplDtoK0PiPi(pD, K0, K+, K-)
 487           // f(D)    = amplDtoK0PiPi(pD, K0, K-, K+)
 488           // Then ...
 489           amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) +
 490             Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
 491         }
 492         else {
 493           amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) +
 494             Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ;
 495         }
 496       }
 497       else if ( ( p->getParent() -> getId() == BM ) ||
 498                 ( p->getParent() -> getId() == B0B ) ) {
 499         Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ;
 500         // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ...
 501         if ( p->getId() == D0 ) {
 502           amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) +
 503             Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ;
 504         }
 505         else {
 506           amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) +
 507             Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
 508         }
 509       }
 510     }
 511     else {
 512       amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
 513     }
 514   }
 515
 516
 517
 518
 519   // Ds -> K K pi
 520   if(_flag==6) {
 521      EvtResonance2 DsKKpiRes1(p4_p, moms3, moms1, 1.0, 0.0, 0.0455, 0.8944, 1, true); // K*(892)
 522      EvtResonance2 DsKKpiRes2(p4_p, moms3, moms1, 1.48, 138., 0.290, 1.414, 0); // K*_0(1430)
 523      EvtFlatte     DsKKpiRes3(p4_p, moms1, moms2, 5.07, 156., 0.965, _kkpi_params); // f_0(980)
 524      EvtResonance2 DsKKpiRes4(p4_p, moms1, moms2, 1.15, -10., 0.00426, 1.019455, 1, true); // phi(1020)
 525      EvtResonance2 DsKKpiRes5(p4_p, moms1, moms2, 1.28, 53., 0.265, 1.350, 0); // f_0(1370)
 526      EvtResonance2 DsKKpiRes6(p4_p, moms1, moms2, 1.19, 87., 0.137, 1.724, 0); // f_0(1710)
 527      amp = DsKKpiRes1.resAmpl() + DsKKpiRes2.resAmpl() + DsKKpiRes3.resAmpl()
 528         + DsKKpiRes4.resAmpl() + DsKKpiRes5.resAmpl() + DsKKpiRes6.resAmpl();
 529
 530   }
 531
 532
 533   vertex(amp);
 534
 535   return ;
 536 }
 537
 538 EvtComplex EvtDDalitz::amplDtoK0PiPi(EvtVector4R p4_p,  EvtVector4R moms1,
 539                                      EvtVector4R moms2, EvtVector4R moms3) {
 540
 541     //K*(892)-
 542     EvtResonance2 DK2piRes1(p4_p,moms1,moms2,1.418,-190.0,0.0508,0.89166,1);
 543     //K0*(1430)
 544     EvtResonance2 DK2piRes2(p4_p,moms1,moms2,1.818,-337.0,0.294 ,1.412  ,0);
 545     //K2*(1430)
 546     EvtResonance2 DK2piRes3(p4_p,moms1,moms2,0.909,  -5.0,0.0985,1.4256 ,2);
 547     //K*(1680)
 548     EvtResonance2 DK2piRes4(p4_p,moms1,moms2,5.091,-166.0,0.322 ,1.717  ,1);
 549     //DCS K*(892)
 550     EvtResonance2 DK2piRes5(p4_p,moms1,moms3,0.100, -19.0,0.0508,0.89166,1);
 551
 552     //Rho
 553     EvtResonance2 DK2piRes6(p4_p,moms3,moms2,0.909,-340.0,0.1502,0.7693,1);
 554     //Omega
 555     EvtResonance2 DK2piRes7(p4_p,moms3,moms2,.0336,-226.0,0.00844,0.78257,1);
 556     //f0(980)
 557     EvtResonance2 DK2piRes8(p4_p,moms3,moms2,0.309,-152.0,0.05,0.977,0);
 558     //f0(1370)
 559     EvtResonance2 DK2piRes9(p4_p,moms3,moms2,1.636,-255.0,0.272,1.31,0);
 560     //f2(1270)
 561     EvtResonance2 DK2piRes10(p4_p,moms3,moms2,0.636,-32.0,0.1851,1.2754,2);
 562
 563     return EvtComplex(1.0,0.0) +
 564       DK2piRes1.resAmpl() + DK2piRes2.resAmpl() +
 565       DK2piRes3.resAmpl() + DK2piRes4.resAmpl() +
 566       DK2piRes5.resAmpl() + DK2piRes6.resAmpl() +
 567       DK2piRes7.resAmpl() + DK2piRes8.resAmpl() +
 568       DK2piRes9.resAmpl() + DK2piRes10.resAmpl();
 569 }
 570
 571 //
 572 // BaBar decay amplitudes for D0->Ks K+ K-
 573 //
 574 // p4_p is D0
 575 // moms1 is K0s
 576 // moms2 is K+
 577 // moms3 is K-
 578 // Amplitudes and phases are taken from BaBar hep-ex/0207089
 579 // with convention : Non Resonant = Amp 1. / Phase 0.
 580
 581 EvtComplex EvtDDalitz::amplDtoK0KK(EvtVector4R p4_p,  EvtVector4R moms1,
 582                                      EvtVector4R moms2, EvtVector4R moms3) {
 583
 584     //phi
 585     EvtResonance DK0KKRes1( p4_p, moms2, moms3, 113.75, -40.0, 0.0043,
 586                             1.019456, 1 ) ;
 587     //a0(980)
 588     EvtResonance DK0KKRes2( p4_p, moms2, moms3, 152.25, 69.0, 0.1196 , 0.9847,
 589                             0 ) ;
 590     //f0(980)
 591     EvtResonance DK0KKRes3( p4_p, moms2, moms3, 30.5, -201.0, 0.05, 0.980 ,
 592                             0 ) ;
 593     //a0(980)+
 594     EvtResonance DK0KKRes4( p4_p, moms1, moms2, 85.75, -93.0, 0.1196 , 0.9847,
 595                             0 ) ;
 596     //a0(980)-
 597     EvtResonance DK0KKRes5( p4_p, moms3, moms1, 8. , -53.0 ,0.1196, 0.9847,
 598                             0 ) ;
 599
 600     return EvtComplex(1.0,0.0) +
 601       DK0KKRes1.resAmpl() + DK0KKRes2.resAmpl() +
 602       DK0KKRes3.resAmpl() + DK0KKRes4.resAmpl() +
 603       DK0KKRes5.resAmpl() ;
 604
 605 }