Line data Source code
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 0 : EvtIntervalDecayAmp()
45 0 : : _probMax(0.), _nScan(0), _fact(0)
46 0 : {}
47 :
48 : EvtIntervalDecayAmp(const EvtIntervalDecayAmp<T>& other)
49 : : _probMax(other._probMax), _nScan(other._nScan),
50 : COPY_PTR(_fact)
51 : {}
52 :
53 : virtual ~EvtIntervalDecayAmp()
54 0 : {
55 0 : delete _fact;
56 0 : }
57 :
58 :
59 : // Initialize model
60 :
61 : virtual void init()
62 : {
63 : // Collect model parameters and parse them
64 :
65 0 : vector<std::string> args;
66 : int i;
67 0 : for(i=0;i<getNArg();i++) args.push_back(getArgStr(i));
68 0 : EvtMultiChannelParser parser;
69 0 : parser.parse(args);
70 :
71 : // Create factory and interval
72 :
73 0 : if(VERBOSE) report(Severity::Info,"EvtGen") << "Create factory and interval" << std::endl;
74 0 : _fact = createFactory(parser);
75 :
76 : // Maximum PDF value over the Dalitz plot can be specified, or a scan
77 : // can be performed.
78 :
79 0 : _probMax = parser.pdfMax();
80 0 : _nScan = parser.nScan();
81 0 : if(VERBOSE) report(Severity::Info,"EvtGen") << "Pdf maximum " << _probMax << std::endl;
82 0 : if(VERBOSE) report(Severity::Info,"EvtGen") << "Scan number " << _nScan << std::endl;
83 0 : }
84 :
85 :
86 : virtual void initProbMax()
87 : {
88 0 : if(0 == _nScan) {
89 :
90 0 : if(_probMax > 0) setProbMax(_probMax);
91 0 : else assert(0);
92 0 : }
93 : else {
94 :
95 : double factor = 1.2; // increase maximum probability by 20%
96 0 : EvtAmpPdf<T> pdf(*_fact->getAmp());
97 0 : EvtPdfSum<T>* pc = _fact->getPC();
98 0 : EvtPdfDiv<T> pdfdiv(pdf,*pc);
99 0 : printf("Sampling %d points to find maximum\n",_nScan);
100 0 : EvtPdfMax<T> x = pdfdiv.findMax(*pc,_nScan);
101 0 : _probMax = factor * x.value();
102 0 : printf("Found maximum %f\n",x.value());
103 0 : printf("Increase to %f\n",_probMax);
104 0 : setProbMax(_probMax);
105 0 : }
106 0 : }
107 :
108 : virtual void decay(EvtParticle *p)
109 : {
110 : // Set things up in most general way
111 :
112 0 : static EvtId B0=EvtPDL::getId("B0");
113 0 : static EvtId B0B=EvtPDL::getId("anti-B0");
114 0 : double t;
115 0 : EvtId other_b;
116 0 : EvtComplex ampl(0.,0.);
117 :
118 : // Sample using pole-compensator pdf
119 :
120 0 : EvtPdfSum<T>* pc = getPC();
121 0 : _x = pc->randomPoint();
122 :
123 0 : if(_fact->isCPModel()) {
124 :
125 : // Time-dependent Dalitz plot changes
126 : // Dec 2005 (ddujmic@slac.stanford.edu)
127 :
128 0 : EvtComplex A = _fact->getAmp()->evaluate(_x);
129 0 : EvtComplex Abar = _fact->getAmpConj()->evaluate(_x);
130 :
131 0 : EvtCPUtil::getInstance()->OtherB(p,t,other_b);
132 :
133 0 : double dm = _fact->dm();
134 0 : double mixAmpli = _fact->mixAmpli();
135 0 : double mixPhase = _fact->mixPhase();
136 0 : EvtComplex qoverp( cos(mixPhase)*mixAmpli, sin(mixPhase)*mixAmpli);
137 0 : EvtComplex poverq( cos(mixPhase)/mixAmpli, -sin(mixPhase)/mixAmpli);
138 :
139 :
140 0 : if (other_b==B0B) ampl = A*cos(dm*t/(2*EvtConst::c)) +
141 0 : EvtComplex(0.,1.)*Abar*sin(dm*t/(2*EvtConst::c))*qoverp;
142 0 : if (other_b==B0) ampl = Abar*cos(dm*t/(2*EvtConst::c)) +
143 0 : EvtComplex(0.,1.)*A*sin(dm*t/(2*EvtConst::c))*poverq;
144 :
145 :
146 0 : }
147 : else {
148 :
149 0 : ampl = amplNonCP(_x);
150 : }
151 :
152 : // Pole-compensate
153 :
154 0 : double comp = sqrt(pc->evaluate(_x));
155 0 : assert(comp > 0);
156 0 : vertex(ampl/comp);
157 :
158 : // Now generate random angles, rotate and setup
159 : // the daughters
160 :
161 0 : std::vector<EvtVector4R> v = initDaughters(_x);
162 :
163 0 : size_t N = p->getNDaug();
164 0 : if(v.size() != N) {
165 :
166 0 : report(Severity::Info,"EvtGen") << "Number of daughters " << N << std::endl;
167 0 : report(Severity::Info,"EvtGen") << "Momentum vector size " << v.size() << std::endl;
168 0 : assert(0);
169 : }
170 :
171 0 : for(size_t i=0;i<N;i++){
172 :
173 0 : p->getDaug(i)->init(getDaugs()[i],v[i]);
174 : }
175 0 : }
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 0 : EvtComplex amplNonCP(const T & x) {return _fact->getAmp()->evaluate(x);}
184 0 : 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 :
|
