TEvtGen/EvtGenModels/EvtBtoKD3P.cxx

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //      $Id: EvtBtoKD3P.cc,v 1.7 2009/02/19 03:23:29 ryd Exp $
   4 //
   5 // Environment:
   6 //      This software is part of the EvtGen package developed jointly
   7 //      for the BaBar and CLEO collaborations. If you use all or part
   8 //      of it, please give an appropriate acknowledgement.
   9 //
  10 // Copyright Information:
  11 //      Copyright (C) 2003, Colorado State University
  12 //
  13 // Module creator:
  14 //      Abi soffer, CSU, 2003
  15 //-----------------------------------------------------------------------
  16 #include "EvtGenBase/EvtPatches.hh"
  17
  18 // Decay model that does the decay B+->D0K, D0->3 psudoscalars
  19
  20 #include <assert.h>
  21
  22 #include "EvtGenModels/EvtBtoKD3P.hh"
  23 #include "EvtGenBase/EvtDecayTable.hh"
  24 #include "EvtGenBase/EvtParticle.hh"
  25 #include "EvtGenBase/EvtId.hh"
  26 #include "EvtGenBase/EvtRandom.hh"
  27 #include "EvtGenModels/EvtPto3P.hh"
  28 #include "EvtGenBase/EvtReport.hh"
  29 #include "EvtGenBase/EvtDalitzPoint.hh"
  30 #include "EvtGenBase/EvtCyclic3.hh"
  31 using std::endl;
  32
  33 //------------------------------------------------------------------
  34 EvtBtoKD3P::EvtBtoKD3P() :
  35   _model1(0),
  36   _model2(0),
  37   _decayedOnce(false)
  38 {
  39 }
  40
  41 //------------------------------------------------------------------
  42 EvtBtoKD3P::EvtBtoKD3P(const EvtBtoKD3P & other){
  43 }
  44
  45 //------------------------------------------------------------------
  46 EvtBtoKD3P::~EvtBtoKD3P(){
  47 }
  48
  49 //------------------------------------------------------------------
  50 EvtDecayBase * EvtBtoKD3P::clone(){
  51   return new EvtBtoKD3P();
  52 }
  53
  54 //------------------------------------------------------------------
  55 std::string EvtBtoKD3P::getName(){
  56   return "BTOKD3P";
  57 }
  58
  59 //------------------------------------------------------------------
  60 void EvtBtoKD3P::init(){
  61   checkNArg(2); // r, phase
  62   checkNDaug(3); // K, D0(allowed), D0(suppressed).
  63                  // The last two daughters are really one particle
  64
  65   // check that the mother and all daughters are scalars:
  66   checkSpinParent  (  EvtSpinType::SCALAR);
  67   checkSpinDaughter(0,EvtSpinType::SCALAR);
  68   checkSpinDaughter(1,EvtSpinType::SCALAR);
  69   checkSpinDaughter(2,EvtSpinType::SCALAR);
  70
  71   // Check that the B dtr types are K D D:
  72
  73   // get the parameters:
  74   _r = getArg(0);
  75   double phase = getArg(1);
  76   _exp = EvtComplex(cos(phase), sin(phase));
  77 }
  78
  79 //------------------------------------------------------------------
  80 void EvtBtoKD3P::initProbMax(){
  81   setProbMax(1); // this is later changed in decay()
  82 }
  83
  84 //------------------------------------------------------------------
  85 void EvtBtoKD3P::decay(EvtParticle *p){
  86   // tell the subclass that we decay the daughter:
  87   _daugsDecayedByParentModel = true;
  88
  89   // the K is the 1st daughter of the B EvtParticle.
  90   // The decay mode of the allowed D (the one produced in b->c decay) is 2nd
  91   // The decay mode of the suppressed D (the one produced in b->u decay) is 3rd
  92   const int KIND = 0;
  93   const int D1IND = 1;
  94   const int D2IND = 2;
  95
  96   // generate kinematics of daughters (K and D):
  97   EvtId tempDaug[2] = {getDaug(KIND), getDaug(D1IND)};
  98   p->initializePhaseSpace(2, tempDaug);
  99
 100   // Get the D daughter particle and the decay models of the allowed
 101   // and suppressed D modes:
 102   EvtParticle * theD = p->getDaug(D1IND);
 103   EvtPto3P * model1 = (EvtPto3P*)(EvtDecayTable::getDecayFunc(theD));
 104
 105   // for the suppressed mode, re-initialize theD as the suppressed D alias:
 106   theD->init(getDaug(D2IND), theD->getP4());
 107   EvtPto3P * model2 = (EvtPto3P*)(EvtDecayTable::getDecayFunc(theD));
 108
 109   // on the first call:
 110   if (false == _decayedOnce) {
 111     _decayedOnce = true;
 112
 113     // store the D decay model pointers:
 114     _model1 = model1;
 115     _model2 = model2;
 116
 117     // check the decay models of the first 2 daughters and that they
 118     // have the same final states:
 119     std::string name1=model1->getName();
 120     std::string name2=model2->getName();
 121
 122     if (name1 != "PTO3P") {
 123       report(ERROR,"EvtGen")
 124         << "D daughters of EvtBtoKD3P decay must decay via the \"PTO3P\" model"
 125         << endl
 126         << "    but found to decay via " << name1.c_str()
 127         << " or " << name2.c_str()
 128         << ". Will terminate execution!" << endl;
 129       assert(0);
 130     }
 131
 132     EvtId * daugs1 = model1->getDaugs();
 133     EvtId * daugs2 = model2->getDaugs();
 134
 135     bool idMatch = true;
 136     int d;
 137     for (d = 0; d < 2; ++d) {
 138       if (daugs1[d] != daugs2[d]) {
 139         idMatch = false;
 140       }
 141     }
 142     if (false == idMatch) {
 143       report(ERROR,"EvtGen")
 144         << "D daughters of EvtBtoKD3P decay must decay to the same final state"
 145         << endl
 146         << "   particles in the same order (not CP-conjugate order)," << endl
 147         << "   but they were found to decay to" << endl;
 148       for (d = 0; d < model1->getNDaug(); ++d) {
 149         report(ERROR,"") << "   " << EvtPDL::name(daugs1[d]).c_str() << " ";
 150       }
 151       report(ERROR,"") << endl;
 152       for (d = 0; d < model1->getNDaug(); ++d) {
 153         report(ERROR,"") << "   "  << EvtPDL::name(daugs2[d]).c_str() << " ";
 154       }
 155       report(ERROR,"") << endl << ". Will terminate execution!" << endl;
 156       assert(0);
 157     }
 158
 159     // estimate the probmax. Need to know the probmax's of the 2
 160     // models for this:
 161     setProbMax(model1->getProbMax(0)
 162                + _r * _r * model2->getProbMax(0)
 163                + 2 * _r * sqrt(model1->getProbMax(0) * model2->getProbMax(0)));
 164
 165   } // end of things to do on the first call
 166
 167   // make sure the models haven't changed since the first call:
 168   if (_model1 != model1 || _model2 != model2) {
 169     report(ERROR,"EvtGen")
 170       << "D daughters of EvtBtoKD3P decay should have only 1 decay modes, "
 171       << endl
 172       << "    but a new decay mode was found after the first call" << endl
 173       << "    Will terminate execution!" << endl;
 174     assert(0);
 175   }
 176
 177   // get the cover function for each of the models and add them up.
 178   // They are summed with coefficients 1 because we are willing to
 179   // take a small inefficiency (~50%) in order to ensure that the
 180   // cover function is large enough without getting into complications
 181   // associated with the smallness of _r:
 182   EvtPdfSum<EvtDalitzPoint> * pc1 = model1->getPC();
 183   EvtPdfSum<EvtDalitzPoint> * pc2 = model2->getPC();
 184   EvtPdfSum<EvtDalitzPoint> pc;
 185   pc.addTerm(1.0, *pc1);
 186   pc.addTerm(1.0, *pc2);
 187
 188   // from this combined cover function, generate the Dalitz point:
 189   EvtDalitzPoint x = pc.randomPoint();
 190
 191   // get the aptitude for each of the models on this point and add them up:
 192   EvtComplex amp1 = model1->amplNonCP(x);
 193   EvtComplex amp2 = model2->amplNonCP(x);
 194   EvtComplex amp = amp1 + amp2 * _r * _exp;
 195
 196   // get the value of the cover function for this point and set the
 197   // relative amplitude for this decay:
 198
 199   double comp = sqrt(pc.evaluate (x));
 200   vertex (amp/comp);
 201
 202   // Make the daughters of theD:
 203   theD->generateMassTree();
 204
 205   // Now generate the p4's of the daughters of theD:
 206   std::vector<EvtVector4R> v = model2->initDaughters(x);
 207
 208   if(v.size() != theD->getNDaug()) {
 209     report(ERROR,"EvtGen")
 210       << "Number of daughters " << theD->getNDaug()
 211       << " != " << "Momentum vector size " << v.size()
 212       << endl
 213       << "     Terminating execution." << endl;
 214     assert(0);
 215   }
 216
 217   // Apply the new p4's to the daughters:
 218   for(unsigned int i=0; i<theD->getNDaug(); ++i){
 219     theD->getDaug(i)->init(model2->getDaugs()[i], v[i]);
 220   }
 221 }
 222