Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtSVVNONCPEIGEN.cpp
CommitLineData
da0e9ce3 1//--------------------------------------------------------------------------
2//
3// Environment:
4// This software is part of the EvtGen package developed jointly
5// for the BaBar and CLEO collaborations. If you use all or part
6// of it, please give an appropriate acknowledgement.
7//
8// Copyright Information: See EvtGen/COPYRIGHT
9// Copyright (C) 2001 Royal Holloway, University of London
10//
11// Module: EvtSVVNONCPEIGEN.cc
12//
13// Description: Routine to decay scalar -> vector vector
14// and has CP violation.
15//
16// This model does all the ckm-suppressed decays and mixing for you. It randomly 'overwrites'
17// any reco or tagging state as set in the Y(4S) decay model (VSS_(B)MIX) with its own generated states.
18//
19// As such, the corresponding dec file requires only one decay-mode description, for example:
20// Decay MyB0
21// 1.000 rho+ MyD*- SVV_NONCPEIGEN dm beta gamma 0.322 0.31 0.941 0 0.107 1.42 0.02 0 0.02 0 0.02 0 ;
22// EndDecay
23// and furthermore Y(4S) only needs to decay to B0's (or B0bar's).
24// The decay above should be a CKM-favored mode (eg. B0->D*-rho+ or B0bar->D*+rho-).
25// All ckm-suppressed decays and the mixing are derived from this line in the ::Decay function.
26//
27// There are 15 or 27 arguments. The first three are dm, phase1
28// and phase2. dm is the B0-B0bar mass difference. Phases 1
29// and 2 are the CKM weak phases relevant for the particular mode,
30// eg for B-->DstRho phase1 is beta and phase2 is gamma.
31//
32// The next arguments are the 2 amplitudes (= 12 input parameters)
33// in the order: A_f, Abar_f. In the example above, the 'A_f' amplitude now
34// stands for the ckm-favored decay 'B0->D*-rho+', and 'Abar_f' stands for 'B0bar->D*-rho+'
35//
36// Each amplitude has its 3 helicity states in the order +, 0, -, which are each
37// specified by a magnitude and a strong phase.
38//
39// The last 2 arguments A_fbar and Abar_fbar (=12 input parameters) are not necessary,
40// but can included if one wants to set them differently from A_f, Abar_f.
41//
42// Mind you that Hbar_+- = H_-+ (ignoring the weak phase, which flips sign).
43// It is custumary to select one set of helicity states (eg H_+-) and to adopt these for
44// the CP-conjugate decays as well (ie. depict Hbar_-+ with H_+-), which is the interpretation
45// we use for the input-parameters above.
46// However, the angular decay in EvtGen is just a formula in which helicity amplitudes are 'plugged' in,
47// making no difference between B0 or B0bar decays. In the model below we (thus) account for the +-
48// flipping between B0 and B0bar.
49//
50//
51// Modification history:
52// Ajit Kurup 9 March 2001 Module created (from EvtSVSNONCPEIGEN)
53// Max Baak 01/16/2004 Fix of Helicity amplitude ordering.
54// Decay also works for B0bar decays.
55//------------------------------------------------------------------------
56//
57#include "EvtGenBase/EvtPatches.hh"
58#include <stdlib.h>
59#include "EvtGenBase/EvtParticle.hh"
60#include "EvtGenBase/EvtRandom.hh"
61#include "EvtGenBase/EvtGenKine.hh"
62#include "EvtGenBase/EvtCPUtil.hh"
63#include "EvtGenBase/EvtPDL.hh"
64#include "EvtGenBase/EvtReport.hh"
65#include "EvtGenBase/EvtVector4C.hh"
66#include "EvtGenModels/EvtSVVNONCPEIGEN.hh"
67#include <string>
68#include "EvtGenModels/EvtSVVHelAmp.hh"
69#include "EvtGenBase/EvtConst.hh"
70
71EvtSVVNONCPEIGEN::~EvtSVVNONCPEIGEN() {}
72
73std::string EvtSVVNONCPEIGEN::getName(){
74
75 return "SVV_NONCPEIGEN";
76
77}
78
79
80EvtDecayBase* EvtSVVNONCPEIGEN::clone(){
81
82 return new EvtSVVNONCPEIGEN;
83
84}
85
86void EvtSVVNONCPEIGEN::init(){
87
88 // check that there are 27 arguments
89 checkNArg(27,15);
90 checkNDaug(2);
91
92 checkSpinDaughter(0,EvtSpinType::VECTOR);
93 checkSpinDaughter(1,EvtSpinType::VECTOR);
94
95 // The ordering of A_f is :
96 // A_f[0-2] = A_f
97 // A_f[3-5] = Abar_f
98 // A_f[6-8] = A_fbar
99 // A_f[9-11] = Abar_fbar
100 //
101 // Each of the 4 amplitudes include the 3 different helicity states in
102 // the order +, 0, -. See more about helicity amplitude ordering in ::decay
103
104 int i=0;
105 int j=(getNArg()-3)/2;
106
107 for(i=0; i<j; ++i){
108 _A_f[i] = getArg((2*i)+3) * EvtComplex( cos(getArg((2*i)+4)),sin(getArg((2*i)+4)) );
109 }
110
111 // If only 6 amplitudes are specified, calculate the last 6 from the first 6:
112 if(6 == j){
113 for(i = 0; i < 3; ++i){
114 _A_f[6+i] = _A_f[3+i];
115 _A_f[9+i] = _A_f[i];
116 }
117 }
118}
119
120void EvtSVVNONCPEIGEN::initProbMax() {
121 double probMax = 0;
122 for (int i = 0; i < 12; ++i){
123 double amp = abs(_A_f[i]);
124 probMax += amp * amp;
125 }
126
127 setProbMax(probMax);
128}
129
130void EvtSVVNONCPEIGEN::decay( EvtParticle *p){
131
132 //added by Lange Jan4,2000
133 static EvtId B0=EvtPDL::getId("B0");
134 static EvtId B0B=EvtPDL::getId("anti-B0");
135
136 double t;
137 EvtId other_b;
138 EvtId daugs[2];
139
140
141 // MB: flip selects the final of the decay
142 int flip = ((p->getId() == B0) ? 0 : 1);
143 daugs[0]=getDaug(0);
144 daugs[1]=getDaug(1);
145 p->initializePhaseSpace(2,daugs);
146
0ca57c2f 147 EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5);
da0e9ce3 148
149 EvtComplex amp[3];
150
151 double dmt2 = getArg(0) * t / (2 * EvtConst::c);
152 double phiCKM = (2.0 * getArg(1) + getArg(2)); // 2b+g
153 EvtComplex ePlusIPhi(cos(phiCKM), sin(phiCKM));
154 EvtComplex eMinusIPhi(cos(-phiCKM), sin(-phiCKM));
155
156 // flip == 0 : D*-rho+
157 // flip == 1 : D*+rho-
158
159 if (!flip) {
160 if (other_b==B0B){
161 // At t=0 we have a B0
162 for (int i=0; i<3; ++i) {
163 amp[i] = _A_f[i]*cos(dmt2) + eMinusIPhi*EvtComplex(0.0,sin(dmt2))*_A_f[i+3];
164 }
165 }
166 if (other_b==B0){
167 // At t=0 we have a B0bar
168 for(int i=0; i<3; ++i) {
169 amp[i] = _A_f[i]*ePlusIPhi*EvtComplex(0.0,sin(dmt2)) + _A_f[i+3]*cos(dmt2);
170 }
171 }
172 } else{
173 if (other_b==B0B){
174 // At t=0 we have a B0
175
176 // M.Baak 01/16/2004
177 // Note: \bar{H}+- = H-+
178 // If one wants to use the correct helicities for B0 and B0bar decays but the same formula-notation (as done in EvtSVV_HelAmp),
179 // count the B0bar helicities backwards. (Equivalently, one could flip the chi angle.)
180
181 for(int i=0; i<3; ++i) {
182 amp[i] = _A_f[8-i]*cos(dmt2) + eMinusIPhi*EvtComplex(0.0,sin(dmt2))*_A_f[11-i];
183 }
184 }
185 if (other_b==B0){
186 // At t=0 we have a B0bar
187 for(int i=0; i<3; ++i) {
188 amp[i] = _A_f[8-i] * ePlusIPhi * EvtComplex(0.0,sin(dmt2)) + _A_f[11-i]*cos(dmt2);
189 }
190 }
191 }
192
193 EvtSVVHelAmp::SVVHel(p,_amp2,daugs[0],daugs[1],amp[0],amp[1],amp[2]);
194
195 return ;
196}
197
0ca57c2f 198std::string EvtSVVNONCPEIGEN::getParamName(int i) {
199 switch(i) {
200 case 0:
201 return "deltaM";
202 case 1:
203 return "weakPhase1";
204 case 2:
205 return "weakPhase2";
206 case 3:
207 return "AfPlusHelAmp";
208 case 4:
209 return "AfPlusHelAmpPhase";
210 case 5:
211 return "AfZeroHelAmp";
212 case 6:
213 return "AfZeroHelAmpPhase";
214 case 7:
215 return "AfMinusHelAmp";
216 case 8:
217 return "AfMinusHelAmpPhase";
218 case 9:
219 return "AbarfPlusHelAmp";
220 case 10:
221 return "AbarfPlusHelAmpPhase";
222 case 11:
223 return "AbarfZeroHelAmp";
224 case 12:
225 return "AbarfZeroHelAmpPhase";
226 case 13:
227 return "AbarfMinusHelAmp";
228 case 14:
229 return "AbarfMinusHelAmpPhase";
230 case 15:
231 return "AfbarPlusHelAmp";
232 case 16:
233 return "AfbarPlusHelAmpPhase";
234 case 17:
235 return "AfbarZeroHelAmp";
236 case 18:
237 return "AfbarZeroHelAmpPhase";
238 case 19:
239 return "AfbarMinusHelAmp";
240 case 20:
241 return "AfbarMinusHelAmpPhase";
242 case 21:
243 return "AbarfbarPlusHelAmp";
244 case 22:
245 return "AbarfbarPlusHelAmpPhase";
246 case 23:
247 return "AbarfbarZeroHelAmp";
248 case 24:
249 return "AbarfbarZeroHelAmpPhase";
250 case 25:
251 return "AbarfbarMinusHelAmp";
252 case 26:
253 return "AbarfbarMinusHelAmpPhase";
254 default:
255 return "";
256 }
257}
258
259std::string EvtSVVNONCPEIGEN::getParamDefault(int i) {
260 switch(i) {
261 case 3:
262 return "1.0";
263 case 4:
264 return "0.0";
265 case 5:
266 return "1.0";
267 case 6:
268 return "0.0";
269 case 7:
270 return "1.0";
271 case 8:
272 return "0.0";
273 case 9:
274 return "1.0";
275 case 10:
276 return "0.0";
277 case 11:
278 return "1.0";
279 case 12:
280 return "0.0";
281 case 13:
282 return "1.0";
283 case 14:
284 return "0.0";
285 default:
286 return "";
287 }
288}