]>
Commit | Line | Data |
---|---|---|
da0e9ce3 | 1 | //----------------------------------------------------------------------- |
2 | // File and Version Information: | |
0ca57c2f | 3 | // $Id: EvtIntervalDecayAmp.hh,v 1.4 2009-03-16 16:39:16 robbep Exp $ |
da0e9ce3 | 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 | ||
0ca57c2f | 112 | static EvtId B0=EvtPDL::getId("B0"); |
113 | static EvtId B0B=EvtPDL::getId("anti-B0"); | |
da0e9ce3 | 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 | ||
0ca57c2f | 131 | EvtCPUtil::getInstance()->OtherB(p,t,other_b); |
da0e9ce3 | 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 |