00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020 #ifndef EVT_INTERVAL_DECAY_AMP
00021 #define EVT_INTERVAL_DECAY_AMP
00022
00023 #define VERBOSE true
00024 #include <iostream>
00025 #include <vector>
00026 #include <string>
00027 #include "EvtGenBase/EvtDecayAmp.hh"
00028 #include "EvtGenBase/EvtParticle.hh"
00029 #include "EvtGenBase/EvtMacros.hh"
00030 #include "EvtGenBase/EvtPdf.hh"
00031 #include "EvtGenBase/EvtAmpFactory.hh"
00032 #include "EvtGenBase/EvtMultiChannelParser.hh"
00033 #include "EvtGenBase/EvtAmpPdf.hh"
00034 #include "EvtGenBase/EvtCPUtil.hh"
00035 #include "EvtGenBase/EvtPDL.hh"
00036 #include "EvtGenBase/EvtCyclic3.hh"
00037
00038 template <class T>
00039 class EvtIntervalDecayAmp : public EvtDecayAmp {
00040
00041 public:
00042
00043 EvtIntervalDecayAmp()
00044 : _probMax(0.), _nScan(0), _fact(0)
00045 {}
00046
00047 EvtIntervalDecayAmp(const EvtIntervalDecayAmp<T>& other)
00048 : _probMax(other._probMax), _nScan(other._nScan),
00049 COPY_PTR(_fact)
00050 {}
00051
00052 virtual ~EvtIntervalDecayAmp()
00053 {
00054 delete _fact;
00055 }
00056
00057
00058
00059
00060 virtual void init()
00061 {
00062
00063
00064 vector<std::string> args;
00065 int i;
00066 for(i=0;i<getNArg();i++) args.push_back(getArgStr(i));
00067 EvtMultiChannelParser parser;
00068 parser.parse(args);
00069
00070
00071
00072 if(VERBOSE) report(INFO,"EvtGen") << "Create factory and interval" << std::endl;
00073 _fact = createFactory(parser);
00074
00075
00076
00077
00078
00079 _probMax = parser.pdfMax();
00080 _nScan = parser.nScan();
00081 if(VERBOSE) report(INFO,"EvtGen") << "Pdf maximum " << _probMax << std::endl;
00082 if(VERBOSE) report(INFO,"EvtGen") << "Scan number " << _nScan << std::endl;
00083 }
00084
00085
00086 virtual void initProbMax()
00087 {
00088 if(0 == _nScan) {
00089
00090 if(_probMax > 0) setProbMax(_probMax);
00091 else assert(0);
00092 }
00093 else {
00094
00095 double factor = 1.2;
00096 EvtAmpPdf<T> pdf(*_fact->getAmp());
00097 EvtPdfSum<T>* pc = _fact->getPC();
00098 EvtPdfDiv<T> pdfdiv(pdf,*pc);
00099 printf("Sampling %d points to find maximum\n",_nScan);
00100 EvtPdfMax<T> x = pdfdiv.findMax(*pc,_nScan);
00101 _probMax = factor * x.value();
00102 printf("Found maximum %f\n",x.value());
00103 printf("Increase to %f\n",_probMax);
00104 setProbMax(_probMax);
00105 }
00106 }
00107
00108 virtual void decay(EvtParticle *p)
00109 {
00110
00111
00112 static EvtId B0=EvtPDL::getId("B0");
00113 static EvtId B0B=EvtPDL::getId("anti-B0");
00114 double t;
00115 EvtId other_b;
00116 EvtComplex ampl(0.,0.);
00117
00118
00119
00120 EvtPdfSum<T>* pc = getPC();
00121 _x = pc->randomPoint();
00122
00123 if(_fact->isCPModel()) {
00124
00125 EvtCPUtil::OtherB(p,t,other_b);
00126 EvtComplex A = _fact->getAmp()->evaluate(_x);
00127 EvtComplex Abar = _fact->getAmpConj()->evaluate(_x);
00128 double dm = _fact->dm();
00129
00130 if (other_b==B0B) ampl=A*cos(dm*t/(2*EvtConst::c))+Abar*sin(dm*t/(2*EvtConst::c));
00131 if (other_b==B0) ampl=A*cos(dm*t/(2*EvtConst::c))-Abar*sin(dm*t/(2*EvtConst::c));
00132 }
00133 else {
00134
00135 ampl = amplNonCP(_x);
00136 }
00137
00138
00139
00140 double comp = sqrt(pc->evaluate(_x));
00141 assert(comp > 0);
00142 vertex(ampl/comp);
00143
00144
00145
00146
00147 std::vector<EvtVector4R> v = initDaughters(_x);
00148
00149 int N = p->getNDaug();
00150 if(v.size() != N) {
00151
00152 report(INFO,"EvtGen") << "Number of daughters " << N << std::endl;
00153 report(INFO,"EvtGen") << "Momentum vector size " << v.size() << std::endl;
00154 assert(0);
00155 }
00156
00157 int i;
00158 for(i=0;i<N;i++){
00159
00160 p->getDaug(i)->init(getDaugs()[i],v[i]);
00161 }
00162 }
00163
00164 virtual EvtAmpFactory<T>* createFactory(const EvtMultiChannelParser& parser) = 0;
00165 virtual std::vector<EvtVector4R> initDaughters(const T& p) const = 0;
00166
00167
00168
00169 const T & x() const {return _x;}
00170 EvtComplex amplNonCP(const T & x) {return _fact->getAmp()->evaluate(x);}
00171 EvtPdfSum<T>* getPC() {return _fact->getPC();}
00172
00173 protected:
00174 double _probMax;
00175 int _nScan;
00176 T _x;
00177
00178 EvtAmpFactory<T>* _fact;
00179 };
00180
00181
00182 #endif
00183
00184
00185
00186