TEvtGen/EvtGenModels/EvtDDalitz.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: 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
  89   _flag=0;
  90   if ( parnum == D0 ) {
  91     //look for either a K- pi+ pi0 or K0bar pi+ pi-
  92     if ( d1==KM && d2==PIP && d3==PI0 ) { _flag=4; _d1=0; _d2=1; _d3=2;}
  93     if ( d1==KM && d3==PIP && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;}
  94     if ( d2==KM && d1==PIP && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;}
  95     if ( d2==KM && d3==PIP && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;}
  96     if ( d3==KM && d1==PIP && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;}
  97     if ( d3==KM && d2==PIP && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;}
  98
  99     if ( d1==KB && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=2; _d3=1;}
 100     if ( d1==KB && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
 101     if ( d2==KB && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 102     if ( d2==KB && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 103     if ( d3==KB && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 104     if ( d3==KB && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 105
 106     if ( d1==KL && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=2; _d3=1;}
 107     if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
 108     if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 109     if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 110     if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 111     if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 112
 113     if ( d1==KS && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=2; _d3=1;}
 114     if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=1; _d3=2;}
 115     if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 116     if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 117     if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 118     if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 119
 120     if ( d1==KS && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 121     if ( d1==KS && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 122     if ( d2==KS && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 123     if ( d2==KS && d3==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 124     if ( d3==KS && d1==KP && d2==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 125     if ( d3==KS && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 126
 127     if ( d1==KL && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 128     if ( d1==KL && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 129     if ( d2==KL && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 130     if ( d2==KL && d3==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 131     if ( d3==KL && d1==KP && d2==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 132     if ( d3==KL && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 133
 134     if ( d1==KB && d2==KP && d3==KM ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 135     if ( d1==KB && d3==KP && d2==KM ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 136     if ( d2==KB && d1==KP && d3==KM ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 137     if ( d2==KB && d3==KP && d1==KM ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 138     if ( d3==KB && d1==KP && d2==KM ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 139     if ( d3==KB && d2==KP && d1==KM ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 140
 141     if ( d1==PIM && d2==PIP && d3==PI0 ) { _flag=12;_d1=0;_d2=1;_d3=2;}
 142     if ( d1==PIM && d3==PIP && d2==PI0 ) { _flag=12;_d1=0;_d2=2;_d3=1;}
 143     if ( d2==PIM && d1==PIP && d3==PI0 ) { _flag=12;_d1=1;_d2=0;_d3=2;}
 144     if ( d2==PIM && d3==PIP && d1==PI0 ) { _flag=12;_d1=2;_d2=0;_d3=1;}
 145     if ( d3==PIM && d1==PIP && d2==PI0 ) { _flag=12;_d1=1;_d2=2;_d3=0;}
 146     if ( d3==PIM && d2==PIP && d1==PI0 ) { _flag=12;_d1=2;_d2=1;_d3=0;}
 147   }
 148   if ( parnum == D0B ) {
 149     //look for either a K+ pi- pi0 or K0 pi+ pi-
 150     if ( d1==KP && d2==PIM && d3==PI0 )  { _flag=4; _d1=0; _d2=1; _d3=2;}
 151     if ( d1==KP && d3==PIM && d2==PI0 ) { _flag=4; _d1=0; _d2=2; _d3=1;}
 152     if ( d2==KP && d1==PIM && d3==PI0 ) { _flag=4; _d1=1; _d2=0; _d3=2;}
 153     if ( d2==KP && d3==PIM && d1==PI0 ) { _flag=4; _d1=1; _d2=2; _d3=0;}
 154     if ( d3==KP && d1==PIM && d2==PI0 ) { _flag=4; _d1=2; _d2=0; _d3=1;}
 155     if ( d3==KP && d2==PIM && d1==PI0 ) { _flag=4; _d1=2; _d2=1; _d3=0;}
 156
 157     if ( d1==K0 && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=1; _d3=2;}
 158     if ( d1==K0 && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
 159     if ( d2==K0 && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 160     if ( d2==K0 && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 161     if ( d3==K0 && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 162     if ( d3==K0 && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 163
 164     if ( d1==KL && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=1; _d3=2;}
 165     if ( d1==KL && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
 166     if ( d2==KL && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 167     if ( d2==KL && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 168     if ( d3==KL && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 169     if ( d3==KL && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 170
 171     if ( d1==KS && d2==PIP && d3==PIM )  { _flag=3; _d1=0; _d2=1; _d3=2;}
 172     if ( d1==KS && d3==PIP && d2==PIM ) { _flag=3; _d1=0; _d2=2; _d3=1;}
 173     if ( d2==KS && d1==PIP && d3==PIM ) { _flag=3; _d1=1; _d2=0; _d3=2;}
 174     if ( d2==KS && d3==PIP && d1==PIM ) { _flag=3; _d1=1; _d2=2; _d3=0;}
 175     if ( d3==KS && d1==PIP && d2==PIM ) { _flag=3; _d1=2; _d2=0; _d3=1;}
 176     if ( d3==KS && d2==PIP && d1==PIM ) { _flag=3; _d1=2; _d2=1; _d3=0;}
 177
 178     if ( d1==KS && d2==KM && d3==KP ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 179     if ( d1==KS && d3==KM && d2==KP ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 180     if ( d2==KS && d1==KM && d3==KP ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 181     if ( d2==KS && d3==KM && d1==KP ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 182     if ( d3==KS && d1==KM && d2==KP ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 183     if ( d3==KS && d2==KM && d1==KP ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 184
 185     if ( d1==KL && d2==KM && d3==KP ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 186     if ( d1==KL && d3==KM && d2==KP ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 187     if ( d2==KL && d1==KM && d3==KP ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 188     if ( d2==KL && d3==KM && d1==KP ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 189     if ( d3==KL && d1==KM && d2==KP ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 190     if ( d3==KL && d2==KM && d1==KP ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 191
 192     if ( d1==K0 && d2==KM && d3==KP ) {_flag=5;_d1=0;_d2=1;_d3=2;}
 193     if ( d1==K0 && d3==KM && d2==KP ) {_flag=5;_d1=0;_d2=2;_d3=1;}
 194     if ( d2==K0 && d1==KM && d3==KP ) {_flag=5;_d1=1;_d2=0;_d3=2;}
 195     if ( d2==K0 && d3==KM && d1==KP ) {_flag=5;_d1=2;_d2=0;_d3=1;}
 196     if ( d3==K0 && d1==KM && d2==KP ) {_flag=5;_d1=1;_d2=2;_d3=0;}
 197     if ( d3==K0 && d2==KM && d1==KP ) {_flag=5;_d1=2;_d2=1;_d3=0;}
 198
 199     if ( d1==PIP && d2==PIM && d3==PI0 ) { _flag=12;_d1=0;_d2=1;_d3=2;}
 200     if ( d1==PIP && d3==PIM && d2==PI0 ) { _flag=12;_d1=0;_d2=2;_d3=1;}
 201     if ( d2==PIP && d1==PIM && d3==PI0 ) { _flag=12;_d1=1;_d2=0;_d3=2;}
 202     if ( d2==PIP && d3==PIM && d1==PI0 ) { _flag=12;_d1=2;_d2=0;_d3=1;}
 203     if ( d3==PIP && d1==PIM && d2==PI0 ) { _flag=12;_d1=1;_d2=2;_d3=0;}
 204     if ( d3==PIP && d2==PIM && d1==PI0 ) { _flag=12;_d1=2;_d2=1;_d3=0;}
 205
 206   }
 207
 208   if ( parnum == DP ) {
 209     //look for K- pi+ pi+
 210     if ( d1==KB && d2==PIP && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 211     if ( d1==KB && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 212     if ( d2==KB && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 213     if ( d2==KB && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 214     if ( d3==KB && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 215     if ( d3==KB && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 216
 217     if ( d1==KL && d2==PIP && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 218     if ( d1==KL && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 219     if ( d2==KL && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 220     if ( d2==KL && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 221     if ( d3==KL && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 222     if ( d3==KL && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 223
 224     if ( d1==KS && d2==PIP && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 225     if ( d1==KS && d3==PIP && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 226     if ( d2==KS && d1==PIP && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 227     if ( d2==KS && d3==PIP && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 228     if ( d3==KS && d1==PIP && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 229     if ( d3==KS && d2==PIP && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 230
 231     if ( d1==KM && d2==PIP && d3==PIP )  { _flag=1; _d1=0; _d2=1; _d3=2;}
 232     if ( d2==KM && d1==PIP && d3==PIP ) { _flag=1; _d1=1; _d2=0; _d3=2;}
 233     if ( d3==KM && d1==PIP && d2==PIP ) { _flag=1; _d1=2; _d2=0; _d3=1;}
 234   }
 235
 236   if ( parnum == DM ) {
 237     //look for K- pi+ pi+
 238     if ( d1==K0 && d2==PIM && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 239     if ( d1==K0 && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 240     if ( d2==K0 && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 241     if ( d2==K0 && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 242     if ( d3==K0 && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 243     if ( d3==K0 && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 244
 245     if ( d1==KL && d2==PIM && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 246     if ( d1==KL && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 247     if ( d2==KL && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 248     if ( d2==KL && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 249     if ( d3==KL && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 250     if ( d3==KL && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 251
 252     if ( d1==KS && d2==PIM && d3==PI0 )  { _flag=2; _d1=0; _d2=1; _d3=2;}
 253     if ( d1==KS && d3==PIM && d2==PI0 ) { _flag=2; _d1=0; _d2=2; _d3=1;}
 254     if ( d2==KS && d1==PIM && d3==PI0 ) { _flag=2; _d1=1; _d2=0; _d3=2;}
 255     if ( d2==KS && d3==PIM && d1==PI0 ) { _flag=2; _d1=1; _d2=2; _d3=0;}
 256     if ( d3==KS && d1==PIM && d2==PI0 ) { _flag=2; _d1=2; _d2=0; _d3=1;}
 257     if ( d3==KS && d2==PIM && d1==PI0 ) { _flag=2; _d1=2; _d2=1; _d3=0;}
 258
 259     if ( d1==KP && d2==PIM && d3==PIM )  { _flag=1; _d1=0; _d2=1; _d3=2;}
 260     if ( d2==KP && d1==PIM && d3==PIM ) { _flag=1; _d1=1; _d2=0; _d3=2;}
 261     if ( d3==KP && d1==PIM && d2==PIM ) { _flag=1; _d1=2; _d2=0; _d3=1;}
 262   }
 263
 264   if ( parnum == DSP ) {
 265      if ( d1==KM && d2==KP && d3==PIP ) { _flag=6; _d1=0; _d2=1; _d3=2; }
 266      if ( d1==KM && d3==KP && d2==PIP ) { _flag=6; _d1=0; _d2=2; _d3=1; }
 267      if ( d2==KM && d1==KP && d3==PIP ) { _flag=6; _d1=1; _d2=0; _d3=2; }
 268      if ( d2==KM && d3==KP && d1==PIP ) { _flag=6; _d1=1; _d2=2; _d3=0; }
 269      if ( d3==KM && d1==KP && d2==PIP ) { _flag=6; _d1=2; _d2=0; _d3=1; }
 270      if ( d3==KM && d2==KP && d1==PIP ) { _flag=6; _d1=2; _d2=1; _d3=0; }
 271
 272      if ( d1==PIM && d2==PIP && d3==KP ) { _flag=9; _d1=0; _d2=1; _d3=2; }
 273      if ( d1==PIM && d3==PIP && d2==KP ) { _flag=9; _d1=0; _d2=2; _d3=1; }
 274      if ( d2==PIM && d1==PIP && d3==KP ) { _flag=9; _d1=1; _d2=0; _d3=2; }
 275      if ( d2==PIM && d3==PIP && d1==KP ) { _flag=9; _d1=1; _d2=2; _d3=0; }
 276      if ( d3==PIM && d1==PIP && d2==KP ) { _flag=9; _d1=2; _d2=0; _d3=1; }
 277      if ( d3==PIM && d2==PIP && d1==KP ) { _flag=9; _d1=2; _d2=1; _d3=0; }
 278
 279      if ( d1==PIM && d2==PIP && d3==PIP ) { _flag=11; _d1=0; _d2=1; _d3=2; }
 280      if ( d1==PIM && d3==PIP && d2==PIP ) { _flag=11; _d1=0; _d2=2; _d3=1; }
 281      if ( d2==PIM && d1==PIP && d3==PIP ) { _flag=11; _d1=1; _d2=0; _d3=2; }
 282      if ( d2==PIM && d3==PIP && d1==PIP ) { _flag=11; _d1=1; _d2=2; _d3=0; }
 283      if ( d3==PIM && d1==PIP && d2==PIP ) { _flag=11; _d1=2; _d2=0; _d3=1; }
 284      if ( d3==PIM && d2==PIP && d1==PIP ) { _flag=11; _d1=2; _d2=1; _d3=0; }
 285   }
 286
 287   if ( parnum == DSM ) {
 288      if ( d1==KP && d2==KM && d3==PIM ) { _flag=6; _d1=0; _d2=1; _d3=2; }
 289      if ( d1==KP && d3==KM && d2==PIM ) { _flag=6; _d1=0; _d2=2; _d3=1; }
 290      if ( d2==KP && d1==KM && d3==PIM ) { _flag=6; _d1=1; _d2=0; _d3=2; }
 291      if ( d2==KP && d3==KM && d1==PIM ) { _flag=6; _d1=1; _d2=2; _d3=0; }
 292      if ( d3==KP && d1==KM && d2==PIM ) { _flag=6; _d1=2; _d2=0; _d3=1; }
 293      if ( d3==KP && d2==KM && d1==PIM ) { _flag=6; _d1=2; _d2=1; _d3=0; }
 294
 295      if ( d1==PIP && d2==PIM && d3==KM ) { _flag=9; _d1=0; _d2=1; _d3=2; }
 296      if ( d1==PIP && d3==PIM && d2==KM ) { _flag=9; _d1=0; _d2=2; _d3=1; }
 297      if ( d2==PIP && d1==PIM && d3==KM ) { _flag=9; _d1=1; _d2=0; _d3=2; }
 298      if ( d2==PIP && d3==PIM && d1==KM ) { _flag=9; _d1=1; _d2=2; _d3=0; }
 299      if ( d3==PIP && d1==PIM && d2==KM ) { _flag=9; _d1=2; _d2=0; _d3=1; }
 300      if ( d3==PIP && d2==PIM && d1==KM ) { _flag=9; _d1=2; _d2=1; _d3=0; }
 301
 302      if ( d1==PIP && d2==PIM && d3==PIM ) { _flag=11; _d1=0; _d2=1; _d3=2; }
 303      if ( d1==PIP && d3==PIM && d2==PIM ) { _flag=11; _d1=0; _d2=2; _d3=1; }
 304      if ( d2==PIP && d1==PIM && d3==PIM ) { _flag=11; _d1=1; _d2=0; _d3=2; }
 305      if ( d2==PIP && d3==PIM && d1==PIM ) { _flag=11; _d1=1; _d2=2; _d3=0; }
 306      if ( d3==PIP && d1==PIM && d2==PIM ) { _flag=11; _d1=2; _d2=0; _d3=1; }
 307      if ( d3==PIP && d2==PIM && d1==PIM ) { _flag=11; _d1=2; _d2=1; _d3=0; }
 308
 309   }
 310   if ( parnum == DP ) {
 311      if ( d1==KM && d2==KP && d3==PIP ) { _flag=7; _d1=0; _d2=1; _d3=2; }
 312      if ( d1==KM && d3==KP && d2==PIP ) { _flag=7; _d1=0; _d2=2; _d3=1; }
 313      if ( d2==KM && d1==KP && d3==PIP ) { _flag=7; _d1=1; _d2=0; _d3=2; }
 314      if ( d2==KM && d3==KP && d1==PIP ) { _flag=7; _d1=2; _d2=0; _d3=1; }
 315      if ( d3==KM && d1==KP && d2==PIP ) { _flag=7; _d1=1; _d2=2; _d3=0; }
 316      if ( d3==KM && d2==KP && d1==PIP ) { _flag=7; _d1=2; _d2=1; _d3=0; }
 317
 318      if ( d1==PIM && d2==PIP && d3==KP ) { _flag=8; _d1=0; _d2=1; _d3=2; }
 319      if ( d1==PIM && d3==PIP && d2==KP ) { _flag=8; _d1=0; _d2=2; _d3=1; }
 320      if ( d2==PIM && d1==PIP && d3==KP ) { _flag=8; _d1=1; _d2=0; _d3=2; }
 321      if ( d2==PIM && d3==PIP && d1==KP ) { _flag=8; _d1=1; _d2=2; _d3=0; }
 322      if ( d3==PIM && d1==PIP && d2==KP ) { _flag=8; _d1=2; _d2=0; _d3=1; }
 323      if ( d3==PIM && d2==PIP && d1==KP ) { _flag=8; _d1=2; _d2=1; _d3=0; }
 324
 325      if ( d1==PIM && d2==PIP && d3==PIP ) { _flag=10; _d1=0; _d2=1; _d3=2; }
 326      if ( d1==PIM && d3==PIP && d2==PIP ) { _flag=10; _d1=0; _d2=2; _d3=1; }
 327      if ( d2==PIM && d1==PIP && d3==PIP ) { _flag=10; _d1=1; _d2=0; _d3=2; }
 328      if ( d2==PIM && d3==PIP && d1==PIP ) { _flag=10; _d1=1; _d2=2; _d3=0; }
 329      if ( d3==PIM && d1==PIP && d2==PIP ) { _flag=10; _d1=2; _d2=0; _d3=1; }
 330      if ( d3==PIM && d2==PIP && d1==PIP ) { _flag=10; _d1=2; _d2=1; _d3=0; }
 331
 332   }
 333   if ( parnum == DM ) {
 334      if ( d1==KP && d2==KM && d3==PIM ) { _flag=7; _d1=0; _d2=1; _d3=2; }
 335      if ( d1==KP && d3==KM && d2==PIM ) { _flag=7; _d1=0; _d2=2; _d3=1; }
 336      if ( d2==KP && d1==KM && d3==PIM ) { _flag=7; _d1=1; _d2=0; _d3=2; }
 337      if ( d2==KP && d3==KM && d1==PIM ) { _flag=7; _d1=2; _d2=0; _d3=1; }
 338      if ( d3==KP && d1==KM && d2==PIM ) { _flag=7; _d1=1; _d2=2; _d3=0; }
 339      if ( d3==KP && d2==KM && d1==PIM ) { _flag=7; _d1=2; _d2=1; _d3=0; }
 340
 341      if ( d1==PIP && d2==PIM && d3==KM ) { _flag=8; _d1=0; _d2=1; _d3=2; }
 342      if ( d1==PIP && d3==PIM && d2==KM ) { _flag=8; _d1=0; _d2=2; _d3=1; }
 343      if ( d2==PIP && d1==PIM && d3==KM ) { _flag=8; _d1=1; _d2=0; _d3=2; }
 344      if ( d2==PIP && d3==PIM && d1==KM ) { _flag=8; _d1=1; _d2=2; _d3=0; }
 345      if ( d3==PIP && d1==PIM && d2==KM ) { _flag=8; _d1=2; _d2=0; _d3=1; }
 346      if ( d3==PIP && d2==PIM && d1==KM ) { _flag=8; _d1=2; _d2=1; _d3=0; }
 347
 348      if ( d1==PIP && d2==PIM && d3==PIM ) { _flag=10; _d1=0; _d2=1; _d3=2; }
 349      if ( d1==PIP && d3==PIM && d2==PIM ) { _flag=10; _d1=0; _d2=2; _d3=1; }
 350      if ( d2==PIP && d1==PIM && d3==PIM ) { _flag=10; _d1=1; _d2=0; _d3=2; }
 351      if ( d2==PIP && d3==PIM && d1==PIM ) { _flag=10; _d1=1; _d2=2; _d3=0; }
 352      if ( d3==PIP && d1==PIM && d2==PIM ) { _flag=10; _d1=2; _d2=0; _d3=1; }
 353      if ( d3==PIP && d2==PIM && d1==PIM ) { _flag=10; _d1=2; _d2=1; _d3=0; }
 354   }
 355
 356   if ( _flag==6) {
 357     _kkpi_params.push_back(EvtFlatteParam(MPI, MPI, 0.406));
 358     _kkpi_params.push_back(EvtFlatteParam(MKP, MKP, 0.800));
 359   }
 360
 361   if ( _flag==0) {
 362     report(ERROR,"EvtGen") << "EvtDDaltiz: Invalid mode."<<endl;
 363     assert(0);
 364   }
 365 }
 366
 367 void EvtDDalitz::initProbMax() {
 368
 369   // probmax different for different modes!
 370
 371   if ( _flag==1 ) {setProbMax(2500.0);}
 372   if ( _flag==2 ) {setProbMax(150.0);}
 373   if ( _flag==3 ) {setProbMax(3000.0);}
 374   if ( _flag==4 ) {setProbMax(600.0);}
 375   if ( _flag==5 ) {setProbMax(2500000.0);}
 376   if ( _flag==6 ) {setProbMax(45000.0);}
 377   if ( _flag==7 ) {setProbMax(35000.0);}
 378   if ( _flag==8 ) {setProbMax(2500.0);}
 379   if ( _flag==9 ) {setProbMax(1700.0);}
 380   if ( _flag==10 ) {setProbMax(1300.0);}
 381   if ( _flag==11 ) {setProbMax(2200.0);}
 382   if ( _flag==12 ) {setProbMax(1000.0);}
 383
 384 }
 385
 386 void EvtDDalitz::decay( EvtParticle *p){
 387
 388   static EvtId BP = EvtPDL::getId("B+");
 389   static EvtId BM = EvtPDL::getId("B-");
 390   static EvtId B0 = EvtPDL::getId("B0");
 391   static EvtId B0B = EvtPDL::getId("anti-B0");
 392
 393   static EvtId D0=EvtPDL::getId("D0");
 394
 395   double oneby2 = 0.707106782;
 396
 397   bool isBToDK=false;
 398   if ( p -> getParent () ) {
 399     EvtId parId = p -> getParent()->getId ();
 400     if ( ( BP == parId ) || ( BM == parId ) || ( B0 == parId ) ||
 401                        ( B0B == parId ) )
 402       if (EvtDecayTable::getInstance()->getDecayFunc(p->getParent())->getName() == "BTODDALITZCPK") isBToDK=true;
 403   }
 404
 405
 406 //same structure for all of these decays
 407
 408   p->initializePhaseSpace(getNDaug(),getDaugs());
 409   EvtVector4R moms1 = p->getDaug(_d1)->getP4();
 410   EvtVector4R moms2 = p->getDaug(_d2)->getP4();
 411   EvtVector4R moms3 = p->getDaug(_d3)->getP4();
 412
 413   EvtVector4R p4_p;
 414   p4_p.set(p->mass(),0.0,0.0,0.0);
 415
 416   EvtComplex amp(1.0,0.0);
 417
 418 //now determine which D and which decay
 419
 420 //data from Anjos et al, Phys.Rev.D 1993, v.48,num.1,p.56 (E691 resuls)
 421 //for D+ -> K- pi+ pi+, and from Adler et al, Phys.Lett. B196 (1987), 107
 422 //(Mark III results) for D+ -> K0bar pi+ pi0.
 423   //CLEO results for D0->k-pi+pi0
 424
 425   if ( _flag==1) {
 426
 427    //  //have a D+ -> K- pi+ pi+ decay, or charge conjugate
 428 //     //Anjos etal e691 - Phys Rev D48, 56 (1993)
 429     // EvtResonance DplusRes11(p4_p,moms1,moms2,0.78,-60.0,0.0498,0.89610,1);
 430 //     EvtResonance DplusRes12(p4_p,moms3,moms1,0.78,-60.0,0.0498,0.89610,1);//K*(892)
 431
 432 //     EvtResonance DplusRes21(p4_p,moms1,moms2,0.53,132.0,0.287,1.429,0);
 433 //     EvtResonance DplusRes22(p4_p,moms3,moms1,0.53,132.0,0.287,1.429,0);//K*(1430)
 434
 435 //     EvtResonance DplusRes31(p4_p,moms1,moms2,0.47,-51.0,0.323,1.714,1);
 436 //     EvtResonance DplusRes32(p4_p,moms3,moms1,0.47,-51.0,0.323,1.714,1);//K*(1680)
 437
 438 //     amp = amp + oneby2*(-DplusRes11.resAmpl()+DplusRes12.resAmpl()) + oneby2*(DplusRes21.resAmpl() + DplusRes22.resAmpl()) + oneby2*(-DplusRes31.resAmpl()+ DplusRes32.resAmpl());
 439
 440
 441 //    EvtResonance DplusRes11(p4_p,moms1,moms2,amp,phase,width,mass,L);
 442     //CLEO-c p15,arxiv:0802.4214v2
 443     EvtResonance2 DplusRes11(p4_p,moms1,moms2,1.0, 0.0, 0.0503, 0.896, 1, true);
 444     EvtResonance2 DplusRes12(p4_p,moms3,moms1,1.0, 0.0, 0.0503, 0.896, 1, true);//K*(892)
 445     EvtResonance2 DplusRes21(p4_p,moms1,moms2,3.0, 49.7-180.0, 0.164, 1.463, 0);
 446     EvtResonance2 DplusRes22(p4_p,moms3,moms1,3.0, 49.7-180.0, 0.164, 1.463, 0);//K*(1430)
 447     EvtResonance2 DplusRes31(p4_p, moms1, moms2, 0.96, -29.9+180.0, 0.109, 1.4324, 2, true);
 448     EvtResonance2 DplusRes32(p4_p, moms3, moms1, 0.96, -29.9+180.0, 0.109, 1.4324, 2, true);// K*_2(1430)
 449     EvtResonance2 DplusRes41(p4_p,moms1,moms2, 6.5, 29.0, 0.323, 1.717, 1, true);
 450     EvtResonance2 DplusRes42(p4_p,moms3,moms1, 6.5, 29.0, 0.323, 1.717, 1, true);//K*(1680)
 451     EvtResonance2 DplusRes51(p4_p,moms1,moms2, 5.01, -163.7+180.0, 0.470, 0.809, 0);
 452     EvtResonance2 DplusRes52(p4_p,moms3,moms1, 5.01, -163.7+180.0, 0.470, 0.809, 0);//kappa(800)
 453     double pi180inv = 1.0/EvtConst::radToDegrees;
 454     amp = EvtComplex(7.4*cos((-18.4+180.0)*pi180inv),7.4*sin((-18.4+180.0)*pi180inv))+ oneby2*(-DplusRes11.resAmpl()+DplusRes12.resAmpl()) + oneby2*(DplusRes21.resAmpl() + DplusRes22.resAmpl()) + oneby2*(DplusRes31.resAmpl()+ DplusRes32.resAmpl()) + oneby2*(-DplusRes41.resAmpl()+ DplusRes42.resAmpl()) + oneby2*(DplusRes51.resAmpl()+ DplusRes52.resAmpl());
 455     //amp = amp+oneby2*(-DplusRes11.resAmpl()+DplusRes12.resAmpl());
 456
 457  }
 458
 459   if ( _flag==2) {
 460
 461 //have a D+ -> K0bar pi+ pi0 decay
 462 //adler etal MarkIII - Phys Lett B196, 107 (1987)
 463 // Results in this paper:
 464 //   Kbar rho+    FitFraction = 68+/-8+/-12    Phase   0
 465 //   Kbar* pi+                  19+/-6+/-6            43+/-23
 466 //   nonres                     13+/-7+/-8           250+/-19
 467 // These numbers below seem not to be exactly the same
 468 // the phases are equiv to -106=254 and 41
 469 //
 470     EvtResonance DplusKpipi0Res1(p4_p,moms2,moms3,1.00,0.00,0.1512,0.7699,1); //rho+
 471     EvtResonance DplusKpipi0Res2(p4_p,moms3,moms1,0.8695,0.7191,0.0498,0.89159,1); //K*0
 472
 473     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);
 474
 475   }
 476
 477   if(_flag==3) {
 478     // D0 -> K0 pi+ pi- + CC
 479     // If it does not come from a B->DK, decay it as D0 or D0bar separatly
 480     // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+
 481     // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi-
 482
 483     if ( isBToDK ) {
 484       // Gamma angle in rad.
 485       double gamma = EvtDecayTable::getInstance()->getDecayFunc( p->getParent() )
 486         -> getArg( 0 )  ;
 487       // Strong phase in rad.
 488       double delta =  EvtDecayTable::getInstance()->getDecayFunc( p->getParent() )
 489         -> getArg( 1 )  ;
 490       // Ratio between B->D0K and B->D0barK
 491       double A     =  EvtDecayTable::getInstance()->getDecayFunc( p->getParent() )
 492         -> getArg( 2 )  ;
 493
 494       EvtComplex Factor( fabs( A ) * cos ( delta ) ,
 495                          fabs( A ) * sin ( delta ) ) ;
 496
 497       if ( ( p->getParent()->getId() == BP ) ||
 498            ( p->getParent()->getId() == B0 ) ) {
 499         // the ratio D/Dbar
 500         Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ;
 501         if ( p->getId() == D0 ) {
 502           // the flavor of the particle has no meaning. But we need
 503           // it to know which daughter is pi+ or pi-
 504           // M( B+ or B0 ) = f(Dbar) + factor * f(D)
 505           // f(Dbar) = amplDtoK0PiPi(pD, K0, pi+, pi-)
 506           // f(D)    = amplDtoK0PiPi(pD, K0, pi-, pi+)
 507           // Then ...
 508           amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) +
 509             Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
 510         }
 511         else {
 512           amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) +
 513             Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ;
 514         }
 515       }
 516       else if ( ( p->getParent() -> getId() == BM ) ||
 517                 ( p->getParent() -> getId() == B0B ) ) {
 518         Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ;
 519         // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ...
 520         if ( p->getId() == D0 ) {
 521           amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) +
 522             Factor * amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) ;
 523         }
 524         else {
 525           amp = amplDtoK0PiPi( p4_p , moms1 , moms3 , moms2 ) +
 526             Factor * amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
 527         }
 528       }
 529     }
 530     else {
 531       amp = amplDtoK0PiPi( p4_p , moms1 , moms2 , moms3 ) ;
 532     }
 533   }
 534
 535
 536   if(_flag==4) {
 537
 538     EvtResonance2 DKpipi0Res1(p4_p,moms2,moms3,1.0  ,0.0   ,0.1507,0.770 ,1); //rho
 539     EvtResonance2 DKpipi0Res2(p4_p,moms1,moms2,0.39, -0.2  ,0.0505,0.8961,1); //k*0
 540     EvtResonance2 DKpipi0Res3(p4_p,moms1,moms3,0.44, 163.0 ,0.050 ,0.8915,1); //k*-
 541
 542     EvtResonance2 DKpipi0Res4(p4_p,moms1,moms3,0.77 ,55.5  ,0.294 ,1.412 ,0); //k01430-
 543     EvtResonance2 DKpipi0Res5(p4_p,moms1,moms2,0.85 ,166.0 ,0.294 ,1.412 ,0); //k01430bar
 544     EvtResonance2 DKpipi0Res6(p4_p,moms2,moms3,2.5  ,171.0 ,0.240 ,1.700 ,1); //rho1700
 545     EvtResonance2 DKpipi0Res7(p4_p,moms1,moms3,2.5  ,103.0 ,0.322 ,1.717 ,1); //K*1680-
 546
 547
 548
 549     double pi180inv = 1.0/EvtConst::radToDegrees;
 550
 551     amp = EvtComplex(1.75*cos(31.2*pi180inv),1.75*sin(31.2*pi180inv))
 552       + DKpipi0Res1.resAmpl() + DKpipi0Res2.resAmpl() + DKpipi0Res3.resAmpl()
 553       + DKpipi0Res4.resAmpl() + DKpipi0Res5.resAmpl()
 554       + DKpipi0Res6.resAmpl()
 555       + DKpipi0Res7.resAmpl();
 556
 557   }
 558
 559   if(_flag==5) {
 560
 561     // D0 -> K0 K+ K- + CC
 562     // If it does not come from a B->DK, decay it as D0 or D0bar separatly
 563     // if p4_p is D0, moms1 is K0, moms2 is pi-, moms3 is pi+
 564     // if p4_p is D0bar, moms1 is K0, moms2 is pi+, moms3 is pi-
 565
 566     if ( isBToDK ){
 567       // Gamma angle in rad.
 568       double gamma = EvtDecayTable::getInstance()->getDecayFunc( p->getParent() )
 569         -> getArg( 0 )  ;
 570       // Strong phase in rad.
 571       double delta =  EvtDecayTable::getInstance()->getDecayFunc( p->getParent() )
 572         -> getArg( 1 )  ;
 573       // Ratio between B->D0K and B->D0barK
 574       double A     =  EvtDecayTable::getInstance()->getDecayFunc( p->getParent() )
 575         -> getArg( 2 )  ;
 576
 577       EvtComplex Factor( fabs( A ) * cos ( delta ) ,
 578                          fabs( A ) * sin ( delta ) ) ;
 579
 580       if ( ( p->getParent()->getId() == BP ) ||
 581            ( p->getParent()->getId() == B0 ) ) {
 582         // the ratio D/Dbar
 583         Factor = Factor * EvtComplex( cos ( gamma ) , sin ( gamma ) ) ;
 584         if ( p->getId() == D0 ) {
 585           // the flavor of the particle has no meaning. But we need
 586           // it to know which daughter is pi+ or pi-
 587           // M( B+ or B0 ) = f(Dbar) + factor * f(D)
 588           // f(Dbar) = amplDtoK0PiPi(pD, K0, K+, K-)
 589           // f(D)    = amplDtoK0PiPi(pD, K0, K-, K+)
 590           // Then ...
 591           amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) +
 592             Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
 593         }
 594         else {
 595           amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) +
 596             Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ;
 597         }
 598       }
 599       else if ( ( p->getParent() -> getId() == BM ) ||
 600                 ( p->getParent() -> getId() == B0B ) ) {
 601         Factor = Factor * EvtComplex( cos ( gamma ) , - sin ( gamma ) ) ;
 602         // here M( B- or B0bar ) = f(D) + factor * f(Dbar) then ...
 603         if ( p->getId() == D0 ) {
 604           amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) +
 605             Factor * amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) ;
 606         }
 607         else {
 608           amp = amplDtoK0KK( p4_p , moms1 , moms3 , moms2 ) +
 609             Factor * amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
 610         }
 611       }
 612     }
 613     else {
 614       amp = amplDtoK0KK( p4_p , moms1 , moms2 , moms3 ) ;
 615     }
 616   }
 617
 618
 619
 620
 621   // Ds -> K K pi
 622   //Babar, arxiv:1011.4190
 623   if(_flag==6) {
 624      EvtResonance2 DsKKpiRes1(p4_p, moms3, moms1, 1.0, 0.0, 0.0455, 0.8944, 1, true); // K*(892)
 625      EvtResonance2 DsKKpiRes2(p4_p, moms3, moms1, 1.48, 138., 0.290, 1.414, 0); // K*_0(1430)
 626      EvtFlatte     DsKKpiRes3(p4_p, moms1, moms2, 5.07, 156., 0.965, _kkpi_params); // f_0(980)
 627      EvtResonance2 DsKKpiRes4(p4_p, moms1, moms2, 1.15, -10., 0.00426, 1.019455, 1, true); // phi(1020)
 628      EvtResonance2 DsKKpiRes5(p4_p, moms1, moms2, 1.28, 53., 0.265, 1.350, 0); // f_0(1370)
 629      EvtResonance2 DsKKpiRes6(p4_p, moms1, moms2, 1.19, 87., 0.137, 1.724, 0); // f_0(1710)
 630      amp = DsKKpiRes1.resAmpl() + DsKKpiRes2.resAmpl() + DsKKpiRes3.resAmpl()
 631         + DsKKpiRes4.resAmpl() + DsKKpiRes5.resAmpl() + DsKKpiRes6.resAmpl();
 632
 633   }
 634
 635   //D+ -> K K pi
 636   //CLEO PRD 78, 072003 (2008) Fit A
 637   if(_flag==7) {
 638     EvtResonance2 DpKKpiRes1(p4_p, moms3, moms1, 1.0, 0.0, 0.0503, 0.8960, 1, true); // K*(892)
 639     EvtResonance2 DpKKpiRes2(p4_p, moms3, moms1, 3.7, 73.0, 0.290, 1.414, 0); // K*_0(1430)
 640     EvtResonance2 DpKKpiRes3(p4_p, moms1, moms2, 1.189, -179.0+180.0, 0.00426, 1.019455, 1, true); // phi(1020)
 641     EvtResonance2 DpKKpiRes4(p4_p, moms1, moms2, 1.72, 123., 0.265, 1.474, 0); // a_0(1450)
 642     EvtResonance2 DpKKpiRes5(p4_p, moms1, moms2, 1.9, -52.0+180.0, 0.15, 1.68, 1, true); // phi(1680)
 643     EvtResonance2 DpKKpiRes6(p4_p, moms3, moms1, 6.4, 150., 0.109, 1.4324, 2, true); // K*_2(1430)
 644     double pi180inv = 1.0/EvtConst::radToDegrees;
 645     amp = EvtComplex(5.1*cos((53.0)*pi180inv),5.1*sin((53.0)*pi180inv)) +
 646       DpKKpiRes1.resAmpl() + DpKKpiRes2.resAmpl() + DpKKpiRes3.resAmpl()
 647       + DpKKpiRes4.resAmpl() + DpKKpiRes5.resAmpl() + DpKKpiRes6.resAmpl();
 648   }
 649
 650 //D+ -> K pi pi WS (DCS)
 651   //FOCUS PLB 601 10 (2004) ; amplitudes there are individually normalized (although not explicit in the paper)
 652   // thus the magnitudes appearing below come from dividing the ones appearing in the paper by the sqrt of the
 653   // integral over the DP of the corresponding squared amplitude. Writing as pi- pi+ K+ so pipi resonances are (12)
 654   // and Kpi resonances are (31); masses and widths corresponds to PDG 2010
 655   if(_flag==8) {
 656     EvtResonance2 DpKpipiDCSRes1(p4_p, moms1, moms2, 1.0, 0.0, 0.149, 0.775, 1, true); // rho(770)
 657     EvtResonance2 DpKpipiDCSRes2(p4_p, moms3, moms1, 1.0971, -167.1, 0.0487, 0.896, 1, true); // K*(890)
 658     EvtResonance2 DpKpipiDCSRes3(p4_p, moms1, moms2, 0.4738, -134.5, 0.059, 0.972, 0); // f0(980) as simple BW
 659     EvtResonance2 DpKpipiDCSRes4(p4_p, moms3, moms1, 2.2688, 54.4, 0.109, 1.432, 2, true); // K*2(1430)
 660     amp = DpKpipiDCSRes1.resAmpl() + DpKpipiDCSRes2.resAmpl() + DpKpipiDCSRes3.resAmpl()
 661       + DpKpipiDCSRes4.resAmpl();
 662   }
 663
 664   //Ds+ -> K pi pi WS (CS)
 665   //FOCUS PLB 601 10 (2004) ; amplitudes there are individually normalized (although not explicit in the paper)
 666   // thus the magnitudes appearing below come from dividing the ones appearing in the paper by the sqrt of the
 667   // integral over the DP of the corresponding squared amplitude. Writing as pi- pi+ K+ so pipi resonances are (12)
 668   // and Kpi resonances are (31); masses and widths corresponds to PDG 2010
 669   // PROBLEM: by simply doing the procedure for D+, the resulting DP and projections do not resemble what is
 670   // in the paper; the best model is by adding 180 to the vector Kpi resonances
 671   if(_flag==9) {
 672     EvtResonance2 DsKpipiCSRes1(p4_p, moms1, moms2, 1.0, 0.0, 0.149, 0.775, 1, true); // rho(770)
 673     EvtResonance2 DsKpipiCSRes2(p4_p, moms3, moms1, 0.7236, -18.3, 0.0487, 0.896, 1, true); // K*(890)
 674     EvtResonance2 DsKpipiCSRes3(p4_p, moms3, moms1, 2.711, 145.2, 0.232, 1.414, 1, true); // K*(1410)
 675     EvtResonance2 DsKpipiCSRes4(p4_p, moms3, moms1, 1.7549, 59.3, 0.270, 1.425, 0); // K*0(1430)
 676     EvtResonance2 DsKpipiCSRes5(p4_p, moms1, moms2, 7.0589, -151.7, 0.400, 1.465, 1, true); // rho(1450)
 677     double pi180inv = 1.0/EvtConst::radToDegrees;
 678     amp = EvtComplex(3.98*cos(43.1*pi180inv),3.98*sin(43.1*pi180inv)) + DsKpipiCSRes1.resAmpl()
 679          + DsKpipiCSRes2.resAmpl() + DsKpipiCSRes3.resAmpl() + DsKpipiCSRes4.resAmpl()
 680          + DsKpipiCSRes5.resAmpl();
 681   }
 682   // D+ -> pi- pi+ pi+  from E791  [PRL 86 770 (2001)]
 683   // masses and widths below correspond to what they used; there, the amplitudes were individually normalized
 684   // (although not explicit) so magnitudes here are obtained after correcting for that
 685   // Breit-Wigner has a factor of (-1) there which changes the relative phase of the NR wrt to the resonances
 686   // thus the NR magnitude is set as negative
 687   if(_flag==10) {
 688     EvtResonance2 DppipipiRes11(p4_p, moms1, moms2, 1.0, 0.0, 0.150, 0.769, 1, true); // rho(770)
 689     EvtResonance2 DppipipiRes12(p4_p, moms3, moms1, 1.0, 0.0, 0.150, 0.769, 1, true); // rho(770)
 690     EvtResonance2 DppipipiRes21(p4_p, moms1, moms2, 2.2811,  205.7, 0.324, 0.478, 0); // sigma(500)
 691     EvtResonance2 DppipipiRes22(p4_p, moms3, moms1, 2.2811,  205.7, 0.324, 0.478, 0); // sigma(500)
 692     EvtResonance2 DppipipiRes31(p4_p, moms1, moms2, 0.4265,  165.0, 0.044, 0.977, 0); // f0(980) simple BW
 693     EvtResonance2 DppipipiRes32(p4_p, moms3, moms1, 0.4265,  165.0, 0.044, 0.977, 0); // f0(980) simple BW
 694     EvtResonance2 DppipipiRes41(p4_p, moms1, moms2, 2.0321,   57.3, 0.185, 1.275, 2, true); // f2(1270)
 695     EvtResonance2 DppipipiRes42(p4_p, moms3, moms1, 2.0321,   57.3, 0.185, 1.275, 2, true); // f2(1270)
 696     EvtResonance2 DppipipiRes51(p4_p, moms1, moms2, 0.7888,  105.4, 0.173, 1.434, 0); // f0(1370)
 697     EvtResonance2 DppipipiRes52(p4_p, moms3, moms1, 0.7888,  105.4, 0.173, 1.434, 0); // f0(1370)
 698     EvtResonance2 DppipipiRes61(p4_p, moms1, moms2, 0.7363,  319.1, 0.310, 1.465, 1, true); // rho(1450)
 699     EvtResonance2 DppipipiRes62(p4_p, moms3, moms1, 0.7363,  319.1, 0.310, 1.465, 1, true); // rho(1450)
 700     double pi180inv = 1.0/EvtConst::radToDegrees;
 701     amp = EvtComplex(-3.98*cos(57.3*pi180inv),-3.98*sin(57.3*pi180inv))
 702         +  (DppipipiRes11.resAmpl() - DppipipiRes12.resAmpl())  //spin1
 703         +  (DppipipiRes21.resAmpl() + DppipipiRes22.resAmpl()) + (DppipipiRes31.resAmpl() + DppipipiRes32.resAmpl())
 704         +  (DppipipiRes41.resAmpl() + DppipipiRes42.resAmpl()) + (DppipipiRes51.resAmpl() + DppipipiRes52.resAmpl())
 705         +  (DppipipiRes61.resAmpl() - DppipipiRes62.resAmpl());  //spin1
 706   }
 707   // Ds+ -> pi- pi+ pi+  from E791  [PRL 86 765 (2001)]
 708   // masses and widths below correspond to what they used; there, the amplitudes were individually normalized
 709   // (although not explicit) so magnitudes here are obtained after correcting for that
 710   // Breit-Wigner has a factor of (-1) there which changes the relative phase of the NR wrt to the resonances
 711   // thus the NR magnitude is set as negative
 712   if(_flag==11) {
 713     EvtResonance2 DspipipiRes11(p4_p, moms1, moms2, 0.288, 109., 0.150, 0.769, 1, true); // rho(770)
 714     EvtResonance2 DspipipiRes12(p4_p, moms3, moms1, 0.288, 109., 0.150, 0.769, 1, true); // rho(770)
 715     EvtResonance2 DspipipiRes21(p4_p, moms1, moms2, 1.0, 0.0, 0.044, 0.977, 0); // f0(980) simple BW
 716     EvtResonance2 DspipipiRes22(p4_p, moms3, moms1, 1.0, 0.0, 0.044, 0.977, 0); // f0(980) simple BW
 717     EvtResonance2 DspipipiRes31(p4_p, moms1, moms2, 1.075, 133., 0.185, 1.275, 2, true); // f2(1270)
 718     EvtResonance2 DspipipiRes32(p4_p, moms3, moms1, 1.075, 133., 0.185, 1.275, 2, true); // f2(1270)
 719     EvtResonance2 DspipipiRes41(p4_p, moms1, moms2, 2.225, 198., 0.173, 1.434, 0); // f0(1370)
 720     EvtResonance2 DspipipiRes42(p4_p, moms3, moms1, 2.225, 198., 0.173, 1.434, 0); // f0(1370)
 721     EvtResonance2 DspipipiRes51(p4_p, moms1, moms2, 1.107, 162., 0.310, 1.465, 1, true); // rho(1450)
 722     EvtResonance2 DspipipiRes52(p4_p, moms3, moms1, 1.107, 162., 0.310, 1.465, 1, true); // rho(1450)
 723     double pi180inv = 1.0/EvtConst::radToDegrees;
 724     amp = EvtComplex(-0.723*cos(181.*pi180inv),-0.723*sin(181.*pi180inv))
 725         +  (DspipipiRes11.resAmpl() - DspipipiRes12.resAmpl())  //spin1
 726         +  (DspipipiRes21.resAmpl() + DspipipiRes22.resAmpl()) + (DspipipiRes31.resAmpl() + DspipipiRes32.resAmpl())
 727         +  (DspipipiRes41.resAmpl() + DspipipiRes42.resAmpl())
 728         +  (DspipipiRes51.resAmpl() - DspipipiRes52.resAmpl());  //spin1
 729   }
 730
 731   //D0 -> pi+pi-pi0
 732   //PRL 99, 251801 (2007)
 733   //arXiv:hep-ex/0703037
 734   if(_flag==12) {
 735     EvtResonance2 DpipipiRes1p(p4_p, moms2, moms3, 1.0, 0.0, 0.149, 0.775, 1, true);//rho+(770)
 736     EvtResonance2 DpipipiRes1(p4_p, moms1, moms2, 0.588, 16.2, 0.149, 0.775, 1, true);//rho0(770)
 737     EvtResonance2 DpipipiRes1m(p4_p, moms3, moms1, 0.714, -2.0, 0.149, 0.775, 1, true);//rho-(770)
 738     EvtResonance2 DpipipiRes2p(p4_p, moms2, moms3, 0.21, -146.0, 0.400, 1.465, 1, true);//rho+(1450)
 739     EvtResonance2 DpipipiRes2(p4_p, moms1, moms2, 0.33, 10.0, 0.400, 1.465, 1, true);//rho0(1450)
 740     EvtResonance2 DpipipiRes2m(p4_p, moms3, moms1, 0.82, 16.0, 0.400, 1.465, 1, true);//rho-(1450)
 741     EvtResonance2 DpipipiRes3p(p4_p, moms2, moms3, 2.25, -17.0, 0.250, 1.720, 1, true);//rho+(1700)
 742     EvtResonance2 DpipipiRes3(p4_p, moms1, moms2, 2.51, -17.0, 0.250, 1.720, 1, true);//rho0(1700)
 743     EvtResonance2 DpipipiRes3m(p4_p, moms3, moms1, 2.00, -50.0, 0.250, 1.720, 1, true);//rho-(1700)
 744     EvtResonance2 DpipipiRes4(p4_p, moms1, moms2, 0.015, -59.0, 0.07, 0.980, 0);//f0(980)
 745     EvtResonance2 DpipipiRes5(p4_p, moms1, moms2, 0.063, 156.0, 0.350, 1.370, 0);//f0(1370)
 746     EvtResonance2 DpipipiRes6(p4_p, moms1, moms2, 0.058, 12.0, 0.109, 1.505, 0);//f0(1500)
 747     EvtResonance2 DpipipiRes7(p4_p, moms1, moms2, 0.112, 51.0, 0.135, 1.720, 0);//f0(1720)
 748     EvtResonance2 DpipipiRes8(p4_p, moms1, moms2, 1.04, -171.0, 0.185, 1.275, 2, true);//f2(1270)
 749     EvtResonance2 DpipipiRes9(p4_p, moms1, moms2, 0.069, 8.0, 0.600, 0.400, 0);//sigma(400)
 750
 751     double pi180inv = 1.0/EvtConst::radToDegrees;
 752     amp = EvtComplex(0.57*cos(-11.0*pi180inv),0.57*sin(-11.0*pi180inv))
 753       + DpipipiRes1p.resAmpl() + DpipipiRes1.resAmpl() + DpipipiRes1m.resAmpl()
 754       + DpipipiRes2p.resAmpl() + DpipipiRes2.resAmpl() + DpipipiRes2m.resAmpl()
 755       + DpipipiRes3p.resAmpl() + DpipipiRes3.resAmpl() + DpipipiRes3m.resAmpl()
 756       + DpipipiRes4.resAmpl() + DpipipiRes5.resAmpl() + DpipipiRes6.resAmpl()
 757       + DpipipiRes7.resAmpl() + DpipipiRes8.resAmpl() + DpipipiRes9.resAmpl();
 758
 759   }
 760
 761   vertex(amp);
 762
 763   return ;
 764 }
 765
 766 EvtComplex EvtDDalitz::amplDtoK0PiPi(EvtVector4R p4_p,  EvtVector4R moms1,
 767                                      EvtVector4R moms2, EvtVector4R moms3) {
 768
 769     //K*(892)-
 770     EvtResonance2 DK2piRes1(p4_p,moms1,moms2,1.418,-190.0,0.0508,0.89166,1);
 771     //K0*(1430)
 772     EvtResonance2 DK2piRes2(p4_p,moms1,moms2,1.818,-337.0,0.294 ,1.412  ,0);
 773     //K2*(1430)
 774     EvtResonance2 DK2piRes3(p4_p,moms1,moms2,0.909,  -5.0,0.0985,1.4256 ,2);
 775     //K*(1680)
 776     EvtResonance2 DK2piRes4(p4_p,moms1,moms2,5.091,-166.0,0.322 ,1.717  ,1);
 777     //DCS K*(892)
 778     EvtResonance2 DK2piRes5(p4_p,moms1,moms3,0.100, -19.0,0.0508,0.89166,1);
 779
 780     //Rho
 781     EvtResonance2 DK2piRes6(p4_p,moms3,moms2,0.909,-340.0,0.1502,0.7693,1);
 782     //Omega
 783     EvtResonance2 DK2piRes7(p4_p,moms3,moms2,.0336,-226.0,0.00844,0.78257,1);
 784     //f0(980)
 785     EvtResonance2 DK2piRes8(p4_p,moms3,moms2,0.309,-152.0,0.05,0.977,0);
 786     //f0(1370)
 787     EvtResonance2 DK2piRes9(p4_p,moms3,moms2,1.636,-255.0,0.272,1.31,0);
 788     //f2(1270)
 789     EvtResonance2 DK2piRes10(p4_p,moms3,moms2,0.636,-32.0,0.1851,1.2754,2);
 790
 791     return EvtComplex(1.0,0.0) +
 792       DK2piRes1.resAmpl() + DK2piRes2.resAmpl() +
 793       DK2piRes3.resAmpl() + DK2piRes4.resAmpl() +
 794       DK2piRes5.resAmpl() + DK2piRes6.resAmpl() +
 795       DK2piRes7.resAmpl() + DK2piRes8.resAmpl() +
 796       DK2piRes9.resAmpl() + DK2piRes10.resAmpl();
 797 }
 798
 799 //
 800 // BaBar decay amplitudes for D0->Ks K+ K-
 801 //
 802 // p4_p is D0
 803 // moms1 is K0s
 804 // moms2 is K+
 805 // moms3 is K-
 806 // Amplitudes and phases are taken from BaBar hep-ex/0207089
 807 // with convention : Non Resonant = Amp 1. / Phase 0.
 808
 809 EvtComplex EvtDDalitz::amplDtoK0KK(EvtVector4R p4_p,  EvtVector4R moms1,
 810                                      EvtVector4R moms2, EvtVector4R moms3) {
 811
 812     //phi
 813     EvtResonance DK0KKRes1( p4_p, moms2, moms3, 113.75, -40.0, 0.0043,
 814                             1.019456, 1 ) ;
 815     //a0(980)
 816     EvtResonance DK0KKRes2( p4_p, moms2, moms3, 152.25, 69.0, 0.1196 , 0.9847,
 817                             0 ) ;
 818     //f0(980)
 819     EvtResonance DK0KKRes3( p4_p, moms2, moms3, 30.5, -201.0, 0.05, 0.980 ,
 820                             0 ) ;
 821     //a0(980)+
 822     EvtResonance DK0KKRes4( p4_p, moms1, moms2, 85.75, -93.0, 0.1196 , 0.9847,
 823                             0 ) ;
 824     //a0(980)-
 825     EvtResonance DK0KKRes5( p4_p, moms3, moms1, 8. , -53.0 ,0.1196, 0.9847,
 826                             0 ) ;
 827
 828     return EvtComplex(1.0,0.0) +
 829       DK0KKRes1.resAmpl() + DK0KKRes2.resAmpl() +
 830       DK0KKRes3.resAmpl() + DK0KKRes4.resAmpl() +
 831       DK0KKRes5.resAmpl() ;
 832
 833 }