Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtIntervalDecayAmp.hh
CommitLineData
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
39template <class T>
40class EvtIntervalDecayAmp : public EvtDecayAmp {
41
42public:
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
186protected:
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