Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtVSSBMixCPT.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) 2002 INFN-Pisa
10//
11// Module: EvtVSSBMixCPT.cc
12//
13// Description:
14// Routine to decay vector-> scalar scalar with coherent BB-like mixing
15// including CPT effects
16// Based on VSSBMIX
17//
18// Modification history:
19//
20// F. Sandrelli, Fernando M-V March 03, 2002
21//
22//------------------------------------------------------------------------
23//
24#include "EvtGenBase/EvtPatches.hh"
25#include <stdlib.h>
26#include "EvtGenBase/EvtConst.hh"
27#include "EvtGenBase/EvtParticle.hh"
28#include "EvtGenBase/EvtGenKine.hh"
29#include "EvtGenBase/EvtPDL.hh"
30#include "EvtGenBase/EvtReport.hh"
31#include "EvtGenBase/EvtVector4C.hh"
32#include "EvtGenModels/EvtVSSBMixCPT.hh"
33#include "EvtGenBase/EvtId.hh"
34#include <string>
35#include "EvtGenBase/EvtRandom.hh"
36using std::endl;
37
38EvtVSSBMixCPT::~EvtVSSBMixCPT() {}
39
40std::string EvtVSSBMixCPT::getName(){
41 return "VSS_BMIX";
42}
43
44
45EvtDecayBase* EvtVSSBMixCPT::clone(){
46 return new EvtVSSBMixCPT;
47}
48
49void EvtVSSBMixCPT::init(){
50
51 if ( getNArg()>4) checkNArg(14,12,8);
52
53 if (getNArg()<1) {
54 report(ERROR,"EvtGen") << "EvtVSSBMix generator expected "
55 << " at least 1 argument (deltam) but found:"<<getNArg()<<endl;
56 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
57 ::abort();
58 }
59 // check that we are asked to produced exactly 2 daughters
60 //4 are allowed if they are aliased..
61 checkNDaug(2,4);
62
63 if ( getNDaug()==4) {
64 if ( getDaug(0)!=getDaug(2)||getDaug(1)!=getDaug(3)){
65 report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator allows "
66 << " 4 daughters only if 1=3 and 2=4"
67 << " (but 3 and 4 are aliased "<<endl;
68 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
69 ::abort();
70 }
71 }
72 // check that we are asked to decay a vector particle into a pair
73 // of scalar particles
74
75 checkSpinParent(EvtSpinType::VECTOR);
76
77 checkSpinDaughter(0,EvtSpinType::SCALAR);
78 checkSpinDaughter(1,EvtSpinType::SCALAR);
79
80 // check that our daughter particles are charge conjugates of each other
81 if(!(EvtPDL::chargeConj(getDaug(0)) == getDaug(1))) {
82 report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
83 << "to be charge conjugate." << endl
84 << " Found " << EvtPDL::name(getDaug(0)).c_str() << " and "
85 << EvtPDL::name(getDaug(1)).c_str() << endl;
86 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
87 ::abort();
88 }
89 // check that both daughter particles have the same lifetime
90 if(EvtPDL::getctau(getDaug(0)) != EvtPDL::getctau(getDaug(1))) {
91 report(ERROR,"EvtGen") << "EvtVSSBMixCPT generator expected daughters "
92 << "to have the same lifetime." << endl
93 << " Found ctau = "
94 << EvtPDL::getctau(getDaug(0)) << " mm and "
95 << EvtPDL::getctau(getDaug(1)) << " mm" << endl;
96 report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
97 ::abort();
98 }
99 // precompute quantities that will be used to generate events
100 // and print out a summary of parameters for this decay
101
102 // mixing frequency in hbar/mm
103 _freq= getArg(0)/EvtConst::c;
104
105 // deltaG
106 double gamma= 1/EvtPDL::getctau(getDaug(0)); // gamma/c (1/mm)
107 _dGamma=0.0;
108 double dgog=0.0;
109 if ( getNArg() > 1 ) {
110 dgog=getArg(1);
111 _dGamma=dgog*gamma;
112 }
113 // q/p
114 _qoverp = EvtComplex(1.0,0.0);
115 if ( getNArg() > 2){
116 _qoverp = EvtComplex(getArg(2),0.0);
117 }
118 if ( getNArg() > 3) {
119 _qoverp = getArg(2)*EvtComplex(cos(getArg(3)),sin(getArg(3)));
120 }
121 _poverq=1.0/_qoverp;
122
123 // decay amplitudes
124 _A_f=EvtComplex(1.0,0.0);
125 _Abar_f=EvtComplex(0.0,0.0);
126 _A_fbar=_Abar_f; // CPT conservation
127 _Abar_fbar=_A_f; // CPT conservation
128 if ( getNArg() > 4){
129 _A_f=getArg(4)*EvtComplex(cos(getArg(5)),sin(getArg(5))); // this allows for DCSD
130 _Abar_f=getArg(6)*EvtComplex(cos(getArg(7)),sin(getArg(7))); // this allows for DCSD
131 if ( getNArg() > 8 ){
132 // CPT violation in decay
133 _A_fbar=getArg(8)*EvtComplex(cos(getArg(9)),sin(getArg(9)));
134 _Abar_fbar=getArg(10)*EvtComplex(cos(getArg(11)),sin(getArg(11)));
135 } else {
136 // CPT conservation in decay
137 _A_fbar=_Abar_f;
138 _Abar_fbar=_A_f;
139 }
140 }
141
142 // CPT violation in mixing
143 _z = EvtComplex(0.0,0.0);
144 if ( getNArg() > 12 ){
145 _z = EvtComplex(getArg(12),getArg(13));
146 }
147
148
149 // some printout
150 double tau= 1e12*EvtPDL::getctau(getDaug(0))/EvtConst::c; // in ps
151 double dm= 1e-12*getArg(0); // B0/anti-B0 mass difference in hbar/ps
152 double x= dm*tau;
153 double y= dgog*0.5; //y=dgamma/(2*gamma)
154 double qop2 = abs(_qoverp*_qoverp);
155 _chib0_b0bar=qop2*(x*x+y*y)/(qop2*(x*x+y*y)+2+x*x-y*y); // does not include CPT in mixing
156 _chib0bar_b0=(1/qop2)*(x*x+y*y)/((1/qop2)*(x*x+y*y)+2+x*x-y*y); // does not include CPT in mixing
157
158 if ( verbose() ) {
159 report(INFO,"EvtGen") << "VSS_BMIXCPT will generate mixing and CPT/CP effects in mixing:"
160 << endl << endl
161 << " " << EvtPDL::name(getParentId()).c_str() << " --> "
162 << EvtPDL::name(getDaug(0)).c_str() << " + "
163 << EvtPDL::name(getDaug(1)).c_str() << endl << endl
164 << "using parameters:" << endl << endl
165 << " delta(m) = " << dm << " hbar/ps" << endl
166 << " _freq = " << _freq << " hbar/mm" << endl
167 << " dgog = " << dgog <<endl
168 << " dGamma = " << _dGamma <<" hbar/mm" <<endl
169 << " q/p = " << _qoverp << endl
170 << " z = " << _z << endl
171 << " tau = " << tau << " ps" << endl
172 << " x = " << x << endl
173 << " chi(B0->B0bar) = " << _chib0_b0bar << endl
174 << " chi(B0bar->B0) = " << _chib0bar_b0 << endl
175 << " Af = " << _A_f << endl
176 << " Abarf = " << _Abar_f << endl
177 << " Afbar = " << _A_fbar << endl
178 << " Abarfbar = " << _Abar_fbar << endl
179 << endl;
180 }
181}
182
183void EvtVSSBMixCPT::initProbMax(){
184 // this value is ok for reasonable values of all the parameters
185 setProbMax(4.0);
186}
187
188void EvtVSSBMixCPT::decay( EvtParticle *p ){
189
190 static EvtId B0=EvtPDL::getId("B0");
191 static EvtId B0B=EvtPDL::getId("anti-B0");
192
193 // generate a final state according to phase space
194
195 double rndm= EvtRandom::random();
196
197 if ( getNDaug()==4) {
198 EvtId tempDaug[2];
199
200 if ( rndm < 0.5 ) { tempDaug[0]=getDaug(0); tempDaug[1]=getDaug(3); }
201 else{ tempDaug[0]=getDaug(2); tempDaug[1]=getDaug(1); }
202
203 p->initializePhaseSpace(2,tempDaug);
204 }
205 else{ //nominal case.
206 p->initializePhaseSpace(2,getDaugs());
207 }
208
209 EvtParticle *s1,*s2;
210
211 s1 = p->getDaug(0);
212 s2 = p->getDaug(1);
213 //delete any daughters - if there are daughters, they
214 //are from the initialization and will be redone later
215 if ( s1->getNDaug() > 0 ) { s1->deleteDaughters();}
216 if ( s2->getNDaug() > 0 ) { s2->deleteDaughters();}
217
218 EvtVector4R p1= s1->getP4();
219 EvtVector4R p2= s2->getP4();
220
221 // throw a random number to decide if this final state should be mixed
222 rndm= EvtRandom::random();
223 int mixed= (rndm < 0.5) ? 1 : 0;
224
225 // if this decay is mixed, choose one of the 2 possible final states
226 // with equal probability (re-using the same random number)
227 if(mixed==1) {
228 EvtId mixedId= (rndm < 0.25) ? getDaug(0) : getDaug(1);
229 EvtId mixedId2= mixedId;
230 if (getNDaug()==4&&rndm<0.25) mixedId2=getDaug(2);
231 if (getNDaug()==4&&rndm>0.25) mixedId2=getDaug(3);
232 s1->init(mixedId, p1);
233 s2->init(mixedId2, p2);
234 }
235
236
237 // if this decay is unmixed, choose one of the 2 possible final states
238 // with equal probability (re-using the same random number)
239 if(mixed==0) {
240 EvtId unmixedId = (rndm < 0.75) ? getDaug(0) : getDaug(1);
241 EvtId unmixedId2= (rndm < 0.75) ? getDaug(1) : getDaug(0);
242 if (getNDaug()==4&&rndm<0.75) unmixedId2=getDaug(3);
243 if (getNDaug()==4&&rndm>0.75) unmixedId2=getDaug(2);
244 s1->init(unmixedId, p1);
245 s2->init(unmixedId2, p2);
246 }
247
248 // choose a decay time for each final state particle using the
249 // lifetime (which must be the same for both particles) in pdt.table
250 // and calculate the lifetime difference for this event
251 s1->setLifetime();
252 s2->setLifetime();
253 double dct= s1->getLifetime() - s2->getLifetime(); // in mm
254
255 // Convention: _dGamma=GammaLight-GammaHeavy
256 EvtComplex exp1(-0.25*_dGamma*dct,0.5*_freq*dct);
257
258 /*
259 //Find the flavor of the B that decayed first.
260 EvtId firstDec = (dct > 0 ) ? s2->getId() : s1->getId();
261
262 //This tags the flavor of the other particle at that time.
263 EvtId stateAtDeltaTeq0 = ( firstDec==B0 ) ? B0B : B0;
264 */
265 EvtId stateAtDeltaTeq0 = (s2->getId()==B0) ? B0B : B0;
266
267 // calculate the oscillation amplitude, based on wether this event is mixed or not
268 EvtComplex osc_amp;
269
270 //define some useful functions: (see BAD #188 eq. 39 for ref.)
271 EvtComplex gp=0.5*(exp(-1.0*exp1)+exp(exp1));
272 EvtComplex gm=0.5*(exp(-1.0*exp1)-exp(exp1));
273 EvtComplex sqz=sqrt(abs(1-_z*_z))*exp(EvtComplex(0,arg(1-_z*_z)/2));
274
275 EvtComplex BB=gp+_z*gm; // <B0|B0(t)>
276 EvtComplex barBB=-sqz*_qoverp*gm; // <B0bar|B0(t)>
277 EvtComplex BbarB=-sqz*_poverq*gm; // <B0|B0bar(t)>
278 EvtComplex barBbarB=gp-_z*gm; // <B0bar|B0bar(t)>
279
280 //
281 if ( !mixed&&stateAtDeltaTeq0==B0 ) {
282 osc_amp= BB*_A_f+barBB*_Abar_f;
283 }
284 if ( !mixed&&stateAtDeltaTeq0==B0B ) {
285 osc_amp= barBbarB*_Abar_fbar+BbarB*_A_fbar;
286 }
287
288 if ( mixed&&stateAtDeltaTeq0==B0 ) {
289 osc_amp=barBB*_Abar_fbar+BB*_A_fbar;
290 }
291 if ( mixed&&stateAtDeltaTeq0==B0B ) {
292 osc_amp=BbarB*_A_f+barBbarB*_Abar_f;
293 }
294
295 // store the amplitudes for each parent spin basis state
296 double norm=1.0/p1.d3mag();
297 vertex(0,norm*osc_amp*p1*(p->eps(0)));
298 vertex(1,norm*osc_amp*p1*(p->eps(1)));
299 vertex(2,norm*osc_amp*p1*(p->eps(2)));
300
301 return ;
302}
303
0ca57c2f 304std::string EvtVSSBMixCPT::getParamName(int i) {
305 switch(i) {
306 case 0:
307 return "deltaM";
308 case 1:
309 return "deltaGammaOverGamma";
310 case 2:
311 return "qOverP";
312 case 3:
313 return "qOverPPhase";
314 case 4:
315 return "Af";
316 case 5:
317 return "AfPhase";
318 case 6:
319 return "Abarf";
320 case 7:
321 return "AbarfPhase";
322 case 8:
323 return "Afbar";
324 case 9:
325 return "AfbarPhase";
326 case 10:
327 return "Abarfbar";
328 case 11:
329 return "AbarfbarPhase";
330 case 12:
331 return "Z";
332 case 13:
333 return "ZPhase";
334 default:
335 return "";
336 }
337}
da0e9ce3 338
0ca57c2f 339std::string EvtVSSBMixCPT::getParamDefault(int i) {
340 switch(i) {
341 case 3:
342 return "0.0";
343 case 4:
344 return "1.0";
345 case 5:
346 return "0.0";
347 case 6:
348 return "1.0";
349 case 7:
350 return "0.0";
351 default:
352 return "";
353 }
354}