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