TEvtGen/EvtGenModels/EvtIntervalDecayAmp.hh

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //      $Id: EvtIntervalDecayAmp.hh,v 1.12 2009/02/15 18:21:51 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) 1998 Caltech, UCSB
  12 //
  13 // Module creator:
  14 //      Alexei Dvoretskii, Caltech, 2001-2002.
  15 //-----------------------------------------------------------------------
  16
  17 // Decay model that uses the "amplitude on an interval"
  18 // templatization
  19
  20 #ifndef EVT_INTERVAL_DECAY_AMP
  21 #define EVT_INTERVAL_DECAY_AMP
  22
  23 #define VERBOSE true
  24 #include <iostream>
  25 #include <vector>
  26 #include <string>
  27 #include "EvtGenBase/EvtDecayAmp.hh"
  28 #include "EvtGenBase/EvtParticle.hh"
  29 #include "EvtGenBase/EvtMacros.hh"
  30 #include "EvtGenBase/EvtPdf.hh"
  31 #include "EvtGenBase/EvtAmpFactory.hh"
  32 #include "EvtGenBase/EvtMultiChannelParser.hh"
  33 #include "EvtGenBase/EvtAmpPdf.hh"
  34 #include "EvtGenBase/EvtCPUtil.hh"
  35 #include "EvtGenBase/EvtPDL.hh"
  36 #include "EvtGenBase/EvtCyclic3.hh"
  37 #include "EvtGenBase/EvtReport.hh"
  38
  39 template <class T>
  40 class EvtIntervalDecayAmp : public  EvtDecayAmp {
  41
  42 public:
  43
  44   EvtIntervalDecayAmp()
  45     : _probMax(0.), _nScan(0), _fact(0)
  46   {}
  47
  48   EvtIntervalDecayAmp(const EvtIntervalDecayAmp<T>& other)
  49     : _probMax(other._probMax), _nScan(other._nScan),
  50       COPY_PTR(_fact)
  51   {}
  52
  53   virtual ~EvtIntervalDecayAmp()
  54   {
  55     delete _fact;
  56   }
  57
  58
  59   // Initialize model
  60
  61   virtual void init()
  62   {
  63     // Collect model parameters and parse them
  64
  65     vector<std::string> args;
  66     int i;
  67     for(i=0;i<getNArg();i++) args.push_back(getArgStr(i));
  68     EvtMultiChannelParser parser;
  69     parser.parse(args);
  70
  71     // Create factory and interval
  72
  73     if(VERBOSE) report(INFO,"EvtGen") << "Create factory and interval" << std::endl;
  74     _fact = createFactory(parser);
  75
  76     // Maximum PDF value over the Dalitz plot can be specified, or a scan
  77     // can be performed.
  78
  79     _probMax = parser.pdfMax();
  80     _nScan = parser.nScan();
  81     if(VERBOSE) report(INFO,"EvtGen") << "Pdf maximum " << _probMax << std::endl;
  82     if(VERBOSE) report(INFO,"EvtGen") << "Scan number " << _nScan << std::endl;
  83   }
  84
  85
  86   virtual void initProbMax()
  87   {
  88     if(0 == _nScan) {
  89
  90       if(_probMax > 0) setProbMax(_probMax);
  91       else assert(0);
  92     }
  93     else {
  94
  95       double factor = 1.2; // increase maximum probability by 20%
  96       EvtAmpPdf<T> pdf(*_fact->getAmp());
  97       EvtPdfSum<T>* pc = _fact->getPC();
  98       EvtPdfDiv<T> pdfdiv(pdf,*pc);
  99       printf("Sampling %d points to find maximum\n",_nScan);
 100       EvtPdfMax<T> x = pdfdiv.findMax(*pc,_nScan);
 101       _probMax = factor * x.value();
 102       printf("Found maximum %f\n",x.value());
 103       printf("Increase to   %f\n",_probMax);
 104       setProbMax(_probMax);
 105     }
 106   }
 107
 108   virtual void decay(EvtParticle *p)
 109   {
 110     // Set things up in most general way
 111
 112     //static EvtId B0=EvtPDL::getId("B0");
 113     //static EvtId B0B=EvtPDL::getId("anti-B0");
 114     EvtId B0=EvtPDL::getId("B0");
 115     EvtId B0B=EvtPDL::getId("anti-B0");
 116     double t;
 117     EvtId other_b;
 118     EvtComplex ampl(0.,0.);
 119
 120     // Sample using pole-compensator pdf
 121
 122     EvtPdfSum<T>* pc = getPC();
 123     _x = pc->randomPoint();
 124
 125     if(_fact->isCPModel()) {
 126
 127       // Time-dependent Dalitz plot changes
 128       // Dec 2005 (ddujmic@slac.stanford.edu)
 129
 130       EvtComplex A    = _fact->getAmp()->evaluate(_x);
 131       EvtComplex Abar = _fact->getAmpConj()->evaluate(_x);
 132
 133       EvtCPUtil::OtherB(p,t,other_b);
 134
 135       double dm = _fact->dm();
 136       double mixAmpli = _fact->mixAmpli();
 137       double mixPhase = _fact->mixPhase();
 138       EvtComplex qoverp( cos(mixPhase)*mixAmpli,  sin(mixPhase)*mixAmpli);
 139       EvtComplex poverq( cos(mixPhase)/mixAmpli, -sin(mixPhase)/mixAmpli);
 140
 141
 142       if (other_b==B0B) ampl = A*cos(dm*t/(2*EvtConst::c))  +
 143                           EvtComplex(0.,1.)*Abar*sin(dm*t/(2*EvtConst::c))*qoverp;
 144       if (other_b==B0)  ampl = Abar*cos(dm*t/(2*EvtConst::c))  +
 145                           EvtComplex(0.,1.)*A*sin(dm*t/(2*EvtConst::c))*poverq;
 146
 147
 148     }
 149     else {
 150
 151       ampl = amplNonCP(_x);
 152     }
 153
 154     // Pole-compensate
 155
 156     double comp = sqrt(pc->evaluate(_x));
 157     assert(comp > 0);
 158     vertex(ampl/comp);
 159
 160     // Now generate random angles, rotate and setup
 161     // the daughters
 162
 163     std::vector<EvtVector4R> v = initDaughters(_x);
 164
 165     size_t N = p->getNDaug();
 166     if(v.size() != N) {
 167
 168       report(INFO,"EvtGen") << "Number of daughters " << N << std::endl;
 169       report(INFO,"EvtGen") << "Momentum vector size " << v.size() << std::endl;
 170       assert(0);
 171     }
 172
 173     for(size_t i=0;i<N;i++){
 174
 175       p->getDaug(i)->init(getDaugs()[i],v[i]);
 176     }
 177   }
 178
 179   virtual EvtAmpFactory<T>* createFactory(const EvtMultiChannelParser& parser) = 0;
 180   virtual std::vector<EvtVector4R> initDaughters(const T& p) const = 0;
 181
 182   // provide access to the decay point and to the amplitude of any decay point.
 183   // this is used by EvtBtoKD3P:
 184   const T & x() const {return _x;}
 185   EvtComplex amplNonCP(const T & x) {return _fact->getAmp()->evaluate(x);}
 186   EvtPdfSum<T>* getPC() {return _fact->getPC();}
 187
 188 protected:
 189   double _probMax;          // Maximum probability
 190   int _nScan;               // Number of points for max prob DP scan
 191   T _x;                     // Decay point
 192
 193   EvtAmpFactory<T>*  _fact; // factory
 194 };
 195
 196
 197 #endif
 198
 199
 200
 201