TEvtGen/EvtGenModels/EvtBtoKD3P.cpp

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