TEvtGen/EvtGen/EvtGenModels/EvtIntervalDecayAmp.hh

   1 //-----------------------------------------------------------------------
   2 // File and Version Information:
   3 //      $Id: EvtIntervalDecayAmp.hh,v 1.4 2009-03-16 16:39:16 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) 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     double t;
 115     EvtId other_b;
 116     EvtComplex ampl(0.,0.);
 117
 118     // Sample using pole-compensator pdf
 119
 120     EvtPdfSum<T>* pc = getPC();
 121     _x = pc->randomPoint();
 122
 123     if(_fact->isCPModel()) {
 124
 125       // Time-dependent Dalitz plot changes
 126       // Dec 2005 (ddujmic@slac.stanford.edu)
 127
 128       EvtComplex A    = _fact->getAmp()->evaluate(_x);
 129       EvtComplex Abar = _fact->getAmpConj()->evaluate(_x);
 130
 131       EvtCPUtil::getInstance()->OtherB(p,t,other_b);
 132
 133       double dm = _fact->dm();
 134       double mixAmpli = _fact->mixAmpli();
 135       double mixPhase = _fact->mixPhase();
 136       EvtComplex qoverp( cos(mixPhase)*mixAmpli,  sin(mixPhase)*mixAmpli);
 137       EvtComplex poverq( cos(mixPhase)/mixAmpli, -sin(mixPhase)/mixAmpli);
 138
 139
 140       if (other_b==B0B) ampl = A*cos(dm*t/(2*EvtConst::c))  +
 141                           EvtComplex(0.,1.)*Abar*sin(dm*t/(2*EvtConst::c))*qoverp;
 142       if (other_b==B0)  ampl = Abar*cos(dm*t/(2*EvtConst::c))  +
 143                           EvtComplex(0.,1.)*A*sin(dm*t/(2*EvtConst::c))*poverq;
 144
 145
 146     }
 147     else {
 148
 149       ampl = amplNonCP(_x);
 150     }
 151
 152     // Pole-compensate
 153
 154     double comp = sqrt(pc->evaluate(_x));
 155     assert(comp > 0);
 156     vertex(ampl/comp);
 157
 158     // Now generate random angles, rotate and setup
 159     // the daughters
 160
 161     std::vector<EvtVector4R> v = initDaughters(_x);
 162
 163     size_t N = p->getNDaug();
 164     if(v.size() != N) {
 165
 166       report(INFO,"EvtGen") << "Number of daughters " << N << std::endl;
 167       report(INFO,"EvtGen") << "Momentum vector size " << v.size() << std::endl;
 168       assert(0);
 169     }
 170
 171     for(size_t i=0;i<N;i++){
 172
 173       p->getDaug(i)->init(getDaugs()[i],v[i]);
 174     }
 175   }
 176
 177   virtual EvtAmpFactory<T>* createFactory(const EvtMultiChannelParser& parser) = 0;
 178   virtual std::vector<EvtVector4R> initDaughters(const T& p) const = 0;
 179
 180   // provide access to the decay point and to the amplitude of any decay point.
 181   // this is used by EvtBtoKD3P:
 182   const T & x() const {return _x;}
 183   EvtComplex amplNonCP(const T & x) {return _fact->getAmp()->evaluate(x);}
 184   EvtPdfSum<T>* getPC() {return _fact->getPC();}
 185
 186 protected:
 187   double _probMax;          // Maximum probability
 188   int _nScan;               // Number of points for max prob DP scan
 189   T _x;                     // Decay point
 190
 191   EvtAmpFactory<T>*  _fact; // factory
 192 };
 193
 194
 195 #endif
 196
 197
 198
 199