TEvtGen/EvtGenModels/EvtBtoXsgammaFlatEnergy.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 // Module: EvtBtoXsgammaKagan.cc
   9 //
  10 // Description:
  11 //       Routine to perform two-body B->Xs,gamma decays with a fixed hadronic
  12 //       mass. For spectrum measurements.
  13 //       The input parameters are 1: the hadronic mass
  14 // Modification history:
  15 //
  16 //      Jim Libby October 11 2002
  17 //------------------------------------------------------------------------
  18 //
  19 #include "EvtGenBase/EvtPatches.hh"
  20
  21 #include <stdlib.h>
  22 #include "EvtGenModels/EvtBtoXsgamma.hh"
  23 #include "EvtGenModels/EvtBtoXsgammaFlatEnergy.hh"
  24 #include "EvtGenBase/EvtRandom.hh"
  25 #include "EvtGenBase/EvtReport.hh"
  26 #include <fstream>
  27 using std::endl;
  28 using std::fstream;
  29
  30 EvtBtoXsgammaFlatEnergy::~EvtBtoXsgammaFlatEnergy(){
  31 }
  32
  33 void EvtBtoXsgammaFlatEnergy::init(int nArg, double* args){
  34
  35   if ((nArg) > 3 || (nArg > 1 && nArg <3)){
  36
  37   report(ERROR,"EvtGen") << "EvtBtoXsgamma generator model "
  38                          << "EvtBtoXsgammaFlatEnergy expected "
  39                          << "either 1(default config) or two arguments but found: "<<nArg<<endl;
  40   report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
  41     ::abort();
  42   }
  43   _mB0=5.2794;
  44   double mPi = 0.140;
  45   double mK = 0.494;
  46   if(nArg == 1){
  47     _eMin = 1.7;
  48     //Invariant mass of Xsd must be greater the m_pi+m_K leads to
  49     //Egamma < (m_B**2-(m_pi+m_k)**2)/(2m_B)
  50     _eMax = (pow(_mB0,2)-pow(mPi+mK,2))/(2.0*_mB0);
  51   }else{
  52     _eMin=args[1];
  53     _eMax=args[2];
  54   }
  55   if (_eMax>(pow(_mB0,2)-pow(mPi+mK,2))/(2.0*_mB0)){
  56     report(ERROR,"EvtGen") << "Emax greater than Kinematic limit" << endl;
  57     report(ERROR,"EvtGen") << "Reset to the kinematic limit" << endl;
  58     report(ERROR,"EvtGen") << "(m_B**2-(m_pi+m_k)**2)/(2m_B)" << endl;
  59     _eMax = (pow(_mB0,2)-pow(mPi+mK,2))/(2.0*_mB0);
  60   }
  61   _eRange=_eMax-_eMin;
  62 }
  63
  64 double EvtBtoXsgammaFlatEnergy::GetMass( int /*Xscode*/ ){
  65
  66   double eGamma = EvtRandom::Flat(_eRange)+_eMin;
  67   double mH = sqrt(pow(_mB0,2)-2.0*_mB0*eGamma);
  68   return mH;
  69 }
  70