TEvtGen/EvtGenBase/EvtDalitzFlatPdf.cxx

   1 #include "EvtGenBase/EvtPatches.hh"
   2 /*******************************************************************************
   3  * Project: BaBar detector at the SLAC PEP-II B-factory
   4  * Package: EvtGenBase
   5  *    File: $Id: EvtDalitzFlatPdf.cc,v 1.4 2004/12/21 19:58:42 ryd Exp $
   6  *  Author: Alexei Dvoretskii, dvoretsk@slac.stanford.edu, 2001-2002
   7  *
   8  * Copyright (C) 2002 Caltech
   9  *******************************************************************************/
  10
  11 #include "EvtGenBase/EvtPatches.hh"
  12 #include "EvtGenBase/EvtDalitzFlatPdf.hh"
  13
  14 EvtDalitzFlatPdf::EvtDalitzFlatPdf(const EvtDalitzPlot& dp)
  15   : EvtPdf<EvtDalitzPoint>(), _dp(dp)
  16 {}
  17
  18 EvtDalitzFlatPdf::EvtDalitzFlatPdf(const EvtDalitzFlatPdf& other)
  19   : EvtPdf<EvtDalitzPoint>(other), _dp(other._dp)
  20 {}
  21
  22 EvtDalitzFlatPdf::~EvtDalitzFlatPdf()
  23 {}
  24
  25 EvtPdf<EvtDalitzPoint>* EvtDalitzFlatPdf::clone() const
  26 {
  27   return new EvtDalitzFlatPdf(*this);
  28 }
  29
  30 double EvtDalitzFlatPdf::pdf(const EvtDalitzPoint&) const
  31 {
  32   return 1.;
  33 }
  34
  35 EvtValError EvtDalitzFlatPdf::compute_integral(int N) const
  36 {
  37   return EvtValError(_dp.getArea(N),0.);
  38 }
  39
  40 EvtDalitzPoint EvtDalitzFlatPdf::randomPoint()
  41 {
  42   // To obtain a uniform distribution generate
  43   // in terms of q's. Generate in a box that circumscribes the
  44   // Dalitz plot. Accept points inside. If there are two
  45   // many unsuccessful attempts it's a hint that the Dalitz plot
  46   // area is tiny compared to the box. It's a pathological
  47   // case. Abort.
  48
  49   EvtCyclic3::Pair pair1 = EvtCyclic3::BC;
  50   EvtCyclic3::Pair pair2 = EvtCyclic3::CA;
  51
  52   int n = 0;
  53   int maxTries = 1000;
  54   while(n++ < maxTries) {
  55
  56     double q1 = EvtRandom::Flat(_dp.qAbsMin(pair1),_dp.qAbsMax(pair2));
  57     double q2 = EvtRandom::Flat(_dp.qAbsMin(pair2),_dp.qAbsMax(pair2));
  58
  59     EvtDalitzCoord point(pair1,q1,pair2,q2);
  60     EvtDalitzPoint x(_dp,point);
  61
  62     if(x.isValid()) return x;
  63   }
  64
  65   printf("No point generated for dalitz plot after %d tries\n",maxTries);
  66   assert(0);
  67 }