Line data Source code
1 : #include "EvtGenModels/EvtLambdaB2LambdaV.hh"
2 : #include "EvtGenBase/EvtRandom.hh"
3 : #include "EvtGenBase/EvtPatches.hh"
4 : #include <stdlib.h>
5 : #include <fstream>
6 : #include <stdio.h>
7 : #include <string>
8 : #include "EvtGenBase/EvtGenKine.hh"
9 : #include "EvtGenBase/EvtParticle.hh"
10 : #include "EvtGenBase/EvtPDL.hh"
11 : #include "EvtGenBase/EvtReport.hh"
12 :
13 : using std::fstream ;
14 : //************************************************************************
15 : //* *
16 : //* Class EvtLambdaB2LambdaV *
17 : //* *
18 : //************************************************************************
19 : //DECLARE_ALGORITHM_FACTORY( EvtLambdaB2LambdaV );
20 :
21 0 : EvtLambdaB2LambdaV::EvtLambdaB2LambdaV()
22 0 : {
23 : //set facility name
24 0 : fname="EvtGen.EvtLambdaB2LambdaV";
25 0 : }
26 :
27 :
28 : //------------------------------------------------------------------------
29 : // Destructor
30 : //------------------------------------------------------------------------
31 : EvtLambdaB2LambdaV::~EvtLambdaB2LambdaV()
32 0 : {}
33 :
34 :
35 : //------------------------------------------------------------------------
36 : // Method 'getName'
37 : //------------------------------------------------------------------------
38 : std::string EvtLambdaB2LambdaV::getName()
39 : {
40 0 : return "LAMBDAB2LAMBDAV";
41 : }
42 :
43 :
44 : //------------------------------------------------------------------------
45 : // Method 'clone'
46 : //------------------------------------------------------------------------
47 : EvtDecayBase* EvtLambdaB2LambdaV::clone()
48 : {
49 0 : return new EvtLambdaB2LambdaV;
50 :
51 0 : }
52 :
53 :
54 : //------------------------------------------------------------------------
55 : // Method 'initProbMax'
56 : //------------------------------------------------------------------------
57 : void EvtLambdaB2LambdaV::initProbMax()
58 : {
59 : //maximum (case where C=0)
60 0 : double Max = 1+fabs(A*B);
61 0 : report(Severity::Debug,fname.c_str())<<" PDF max value : "<<Max<<std::endl;
62 0 : setProbMax(Max);
63 0 : }
64 :
65 :
66 : //------------------------------------------------------------------------
67 : // Method 'init'
68 : //------------------------------------------------------------------------
69 : void EvtLambdaB2LambdaV::init()
70 : {
71 : bool antiparticle=false;
72 :
73 : //introduction
74 0 : report(Severity::Debug,fname.c_str())<< "*************************************************"<<std::endl;
75 0 : report(Severity::Debug,fname.c_str())<< "* Event Model Class : EvtLambdaB2LambdaV *"<<std::endl;
76 0 : report(Severity::Debug,fname.c_str())<< "*************************************************"<<std::endl;
77 :
78 : //check the number of arguments
79 0 : checkNArg(2);
80 :
81 :
82 : //check the number of daughters
83 0 : checkNDaug(2);
84 :
85 0 : const EvtId Id_mother = getParentId();
86 0 : const EvtId Id_daug1 = getDaug(0);
87 0 : const EvtId Id_daug2 = getDaug(1);
88 :
89 :
90 : //lambdab identification
91 0 : if (Id_mother==EvtPDL::getId("Lambda_b0"))
92 : {
93 : antiparticle=false;
94 0 : }
95 0 : else if (Id_mother==EvtPDL::getId("anti-Lambda_b0"))
96 : {
97 : antiparticle=true;
98 : }
99 : else
100 : {
101 0 : report(Severity::Error,fname.c_str())<<" Mother is not a Lambda_b0 or an anti-Lambda_b0, but a "
102 0 : <<EvtPDL::name(Id_mother)<<std::endl;
103 0 : abort();
104 : }
105 :
106 : //lambda
107 0 : if ( !(Id_daug1==EvtPDL::getId("Lambda0") && !antiparticle ) && !(Id_daug1==EvtPDL::getId("anti-Lambda0") && antiparticle) )
108 : {
109 0 : if (!antiparticle)
110 : {
111 0 : report(Severity::Error,fname.c_str()) << " Daughter1 is not a Lambda0, but a "
112 0 : << EvtPDL::name(Id_daug1)<<std::endl;
113 0 : }
114 : else
115 0 : { report(Severity::Error,fname.c_str()) << " Daughter1 is not an anti-Lambda0, but a "
116 0 : << EvtPDL::name(Id_daug1)<<std::endl;
117 : }
118 0 : abort();
119 : }
120 :
121 :
122 : //identification meson V
123 0 : if (getArg(1)==1) Vtype=VID::JPSI;
124 0 : else if (getArg(1)==2) Vtype=VID::RHO;
125 0 : else if (getArg(1)==3) Vtype=VID::OMEGA;
126 0 : else if (getArg(1)==4) Vtype=VID::RHO_OMEGA_MIXING;
127 : else
128 : {
129 0 : report(Severity::Error,fname.c_str()) << " Vtype " <<getArg(1)<<" is unknown"<<std::endl;
130 0 : abort();
131 : }
132 :
133 :
134 : //check vector meson V
135 0 : if (Id_daug2==EvtPDL::getId("J/psi") && Vtype==VID::JPSI)
136 : {
137 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : Lambda_b0 -> Lambda J/psi"<<std::endl;
138 0 : else report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda J/psi"<<std::endl;
139 : }
140 0 : else if (Id_daug2==EvtPDL::getId("rho0") && Vtype==VID::RHO )
141 : {
142 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : Lambda_b0 -> Lambda rho0"<<std::endl;
143 0 : else report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda rho0"<<std::endl;
144 : }
145 0 : else if (Id_daug2==EvtPDL::getId("omega") && Vtype==VID::OMEGA)
146 : {
147 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : Lambda_b0 -> Lambda omega"<<std::endl;
148 0 : else report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : anti-Lambda_b0 -> anti-Lambda omega"<<std::endl;
149 : }
150 0 : else if ((Id_daug2==EvtPDL::getId("omega") || Id_daug2==EvtPDL::getId("rho0") )&& Vtype==VID::RHO_OMEGA_MIXING)
151 : {
152 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : "
153 0 : <<"Lambda_b0 -> Lambda rho-omega-mixing"<<std::endl;
154 0 : else report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : "
155 0 : <<"anti-Lambda_b0 -> anti-Lambda rho-omega-mixing"<<std::endl;
156 : }
157 :
158 : else
159 : {
160 0 : report(Severity::Error,fname.c_str())<<" Daughter2 is not a J/psi, phi or rho0 but a "
161 0 : <<EvtPDL::name(Id_daug2)<<std::endl;
162 0 : abort();
163 : }
164 :
165 : //fix dynamics parameters
166 0 : B = (double) getArg(0);
167 0 : C = EvtComplex((sqrt(2.)/2.),(sqrt(2.)/2.));
168 0 : switch(Vtype)
169 : {
170 0 : case VID::JPSI : A = 0.490; break;
171 : case VID::RHO :
172 : case VID::OMEGA :
173 0 : case VID::RHO_OMEGA_MIXING : A = 0.194; break;
174 0 : default : A = 0; break;
175 : }
176 0 : report(Severity::Debug,fname.c_str())<<" LambdaB decay parameters : "<<std::endl;
177 0 : report(Severity::Debug,fname.c_str())<<" - lambda asymmetry A = "<<A<<std::endl;
178 0 : report(Severity::Debug,fname.c_str())<<" - lambdab polarisation B = "<<B<<std::endl;
179 0 : report(Severity::Debug,fname.c_str())<<" - lambdab density matrix rho+- C = "<<C<<std::endl;
180 :
181 :
182 :
183 :
184 0 : }
185 :
186 :
187 : //------------------------------------------------------------------------
188 : // Method 'decay'
189 : //------------------------------------------------------------------------
190 : void EvtLambdaB2LambdaV::decay( EvtParticle *lambdab)
191 : {
192 0 : lambdab->makeDaughters(getNDaug(),getDaugs());
193 :
194 : //get lambda and V particles
195 0 : EvtParticle * lambda = lambdab->getDaug(0);
196 0 : EvtParticle * V = lambdab->getDaug(1);
197 :
198 : //get resonance masses
199 : // - LAMBDAB -> mass given by the preceding class
200 : // - LAMBDA -> nominal mass
201 : // - V -> getVmass method
202 0 : double MASS_LAMBDAB = lambdab->mass();
203 0 : double MASS_LAMBDA = EvtPDL::getMeanMass(EvtPDL::getId("Lambda0"));
204 0 : double MASS_V = getVMass(MASS_LAMBDAB,MASS_LAMBDA);
205 :
206 : //generate random angles
207 0 : double phi = EvtRandom::Flat(0,2*EvtConst::pi);
208 0 : double theta = acos( EvtRandom::Flat(-1,+1));
209 0 : report(Severity::Debug,fname.c_str())<<" Angular angles : theta = "<<theta
210 0 : <<" ; phi = "<<phi<<std::endl;
211 : //computate resonance quadrivectors
212 0 : double E_lambda = (MASS_LAMBDAB*MASS_LAMBDAB + MASS_LAMBDA*MASS_LAMBDA - MASS_V*MASS_V)
213 0 : /(2*MASS_LAMBDAB);
214 0 : double E_V = (MASS_LAMBDAB*MASS_LAMBDAB + MASS_V*MASS_V - MASS_LAMBDA*MASS_LAMBDA)
215 0 : /(2*MASS_LAMBDAB);
216 0 : double P = sqrt(E_lambda*E_lambda-lambda->mass()*lambda->mass());
217 :
218 :
219 :
220 :
221 0 : EvtVector4R P_lambdab=lambdab->getP4();
222 :
223 0 : double px = P_lambdab.get(1);
224 0 : double py = P_lambdab.get(2);
225 0 : double pz = P_lambdab.get(3);
226 0 : double E = P_lambdab.get(0);
227 0 : report(Severity::Info,fname.c_str())<<"E of lambdab: "<< P_lambdab.get(0)<<std::endl;
228 0 : report(Severity::Info,fname.c_str())<<"E of lambdab: "<< E<<std::endl;
229 :
230 :
231 0 : EvtVector4R q_lambdab2 (E,
232 0 : ((1/(sqrt(pow(px,2)+pow(py,2))))*((px*(px))+(py*(py)))),
233 0 : ((1/(sqrt(pow(px,2)+pow(py,2))))*(-((py)*(px))+(px*(py)))),
234 : (pz));
235 :
236 0 : EvtVector4R q_lambdab3 (E,
237 0 : q_lambdab2.get(3),
238 0 : q_lambdab2.get(1),
239 0 : q_lambdab2.get(2));
240 :
241 :
242 0 : EvtVector4R q_lambda0 (E_lambda,
243 0 : P*sin(theta)*cos(phi),
244 0 : P*sin(theta)*sin(phi),
245 0 : P*cos(theta) );
246 :
247 0 : EvtVector4R q_V0 (E_V,
248 0 : -P*sin(theta)*cos(phi),
249 0 : -P*sin(theta)*sin(phi),
250 0 : -P*cos(theta) );
251 :
252 :
253 0 : EvtVector4R q_lambda1 (E_lambda,
254 0 : q_lambda0.get(2),
255 0 : q_lambda0.get(3),
256 0 : q_lambda0.get(1) );
257 :
258 0 : EvtVector4R q_V1 (E_V,
259 0 : q_V0.get(2),
260 0 : q_V0.get(3),
261 0 : q_V0.get(1) );
262 :
263 0 : EvtVector4R q_lambda (E_lambda,
264 0 : ((1/(sqrt(pow(px,2)+pow(py,2))))*((px*(q_lambda1.get(1))) - (py*(q_lambda1.get(2))))),
265 0 : ((1/(sqrt(pow(px,2)+pow(py,2))))*((py*(q_lambda1.get(1))) + (px*(q_lambda1.get(2))))),
266 0 : (q_lambda1.get(3)));
267 :
268 :
269 0 : EvtVector4R q_V (E_V,
270 0 : ((1/(sqrt(pow(px,2)+pow(py,2))))*((px*(q_V1.get(1))) - (py*(q_V1.get(2))))),
271 0 : ((1/(sqrt(pow(px,2)+pow(py,2))))*((py*(q_V1.get(1))) + (px*(q_V1.get(2))))),
272 0 : (q_V1.get(3)));
273 :
274 :
275 :
276 :
277 :
278 0 : lambda->getP4();
279 0 : V->getP4();
280 0 : report(Severity::Info,fname.c_str())<<" LambdaB px: "<<px<<std::endl;
281 0 : report(Severity::Info,fname.c_str())<<" LambdaB py: "<<py<<std::endl;
282 0 : report(Severity::Info,fname.c_str())<<" LambdaB pz: "<<pz<<std::endl;
283 0 : report(Severity::Info,fname.c_str())<<" LambdaB E: "<<E<<std::endl;
284 :
285 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 E: "<<q_lambdab3.get(0)<<std::endl;
286 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 px: "<<q_lambda0.get(1)<<std::endl;
287 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 py: "<<q_lambda0.get(2)<<std::endl;
288 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 pz: "<<q_lambda0.get(3)<<std::endl;
289 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 E: "<<q_lambda0.get(0)<<std::endl;
290 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 px: "<<q_lambda1.get(1)<<std::endl;
291 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 py: "<<q_lambda1.get(2)<<std::endl;
292 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 pz: "<<q_lambda1.get(3)<<std::endl;
293 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 E: "<<q_lambda1.get(0)<<std::endl;
294 0 : report(Severity::Info,fname.c_str())<<" Lambda px: "<<q_lambda.get(1)<<std::endl;
295 0 : report(Severity::Info,fname.c_str())<<" Lambda py: "<<q_lambda.get(2)<<std::endl;
296 0 : report(Severity::Info,fname.c_str())<<" Lambda pz: "<<q_lambda.get(3)<<std::endl;
297 0 : report(Severity::Info,fname.c_str())<<" Lambda E: "<<q_lambda0.get(3)<<std::endl;
298 0 : report(Severity::Info,fname.c_str())<<" V 0 px: "<<q_V0.get(1)<<std::endl;
299 0 : report(Severity::Info,fname.c_str())<<" V 0 py: "<<q_V0.get(2)<<std::endl;
300 0 : report(Severity::Info,fname.c_str())<<" V 0 pz: "<<q_V0.get(3)<<std::endl;
301 0 : report(Severity::Info,fname.c_str())<<" V 0 E: "<<q_V0.get(0)<<std::endl;
302 0 : report(Severity::Info,fname.c_str())<<" V 1 px: "<<q_V1.get(1)<<std::endl;
303 0 : report(Severity::Info,fname.c_str())<<" V 1 py: "<<q_V1.get(2)<<std::endl;
304 0 : report(Severity::Info,fname.c_str())<<" V 1 pz: "<<q_V1.get(3)<<std::endl;
305 0 : report(Severity::Info,fname.c_str())<<" V 1 E: "<<q_V1.get(0)<<std::endl;
306 0 : report(Severity::Info,fname.c_str())<<" V px: "<<q_V.get(1)<<std::endl;
307 0 : report(Severity::Info,fname.c_str())<<" V py: "<<q_V.get(2)<<std::endl;
308 0 : report(Severity::Info,fname.c_str())<<" V pz: "<<q_V.get(3)<<std::endl;
309 0 : report(Severity::Info,fname.c_str())<<" V E: "<<q_V0.get(3)<<std::endl;
310 : //set quadrivectors to particles
311 0 : lambda ->init(getDaugs()[0],q_lambda);
312 0 : V ->init(getDaugs()[1],q_V );
313 :
314 : //computate pdf
315 0 : double pdf = 1 + A*B*cos(theta) + 2*A*real(C*EvtComplex(cos(phi),sin(phi)))*sin(theta);
316 :
317 0 : report(Severity::Debug,fname.c_str())<<" LambdaB decay pdf value : "<<pdf<<std::endl;
318 : //set probability
319 0 : setProb(pdf);
320 :
321 : return;
322 0 : }
323 :
324 :
325 : //------------------------------------------------------------------------
326 : // Method 'BreitWignerRelPDF'
327 : //------------------------------------------------------------------------
328 : double EvtLambdaB2LambdaV::BreitWignerRelPDF(double m,double _m0, double _g0)
329 : {
330 0 : double a = (_m0 * _g0) * (_m0 * _g0);
331 0 : double b = (m*m - _m0*_m0)*(m*m - _m0*_m0);
332 0 : return a/(b+a);
333 : }
334 :
335 :
336 : //------------------------------------------------------------------------
337 : // Method 'RhoOmegaMixingPDF'
338 : //------------------------------------------------------------------------
339 : double EvtLambdaB2LambdaV::RhoOmegaMixingPDF(double m, double _mr, double _gr, double _mo, double _go)
340 : {
341 0 : double a = m*m - _mr*_mr;
342 0 : double b = m*m - _mo*_mo;
343 0 : double c = _gr*_mr;
344 0 : double d = _go*_mo;
345 : double re_pi = -3500e-6; //GeV^2
346 : double im_pi = -300e-6; //GeV^2
347 : double va_pi = re_pi+im_pi;
348 :
349 : //computate pdf value
350 0 : double f = 1/(a*a+c*c) * ( 1+
351 0 : (va_pi*va_pi+2*b*re_pi+2*d*im_pi)/(b*b+d*d));
352 :
353 : //computate pdf max value
354 : a = 0;
355 0 : b = _mr*_mr - _mo*_mo;
356 :
357 0 : double maxi = 1/(a*a+c*c) * ( 1+
358 0 : (va_pi*va_pi+2*b*re_pi+2*d*im_pi)/(b*b+d*d));
359 :
360 0 : return f/maxi;
361 : }
362 :
363 :
364 : //------------------------------------------------------------------------
365 : // Method 'getVMass'
366 : //------------------------------------------------------------------------
367 : double EvtLambdaB2LambdaV::getVMass(double MASS_LAMBDAB, double MASS_LAMBDA)
368 : {
369 : //JPSI case
370 0 : if (Vtype==VID::JPSI)
371 : {
372 0 : return EvtPDL::getMass(EvtPDL::getId("J/psi"));
373 : }
374 :
375 : //RHO,OMEGA,MIXING case
376 : else
377 : {
378 : //parameters
379 0 : double MASS_RHO = EvtPDL::getMeanMass(EvtPDL::getId("rho0"));
380 0 : double MASS_OMEGA = EvtPDL::getMeanMass(EvtPDL::getId("omega"));
381 0 : double WIDTH_RHO = EvtPDL::getWidth(EvtPDL::getId("rho0"));
382 0 : double WIDTH_OMEGA = EvtPDL::getWidth(EvtPDL::getId("omega"));
383 0 : double MASS_PION = EvtPDL::getMeanMass(EvtPDL::getId("pi-"));
384 0 : double _max = MASS_LAMBDAB - MASS_LAMBDA;
385 0 : double _min = 2*MASS_PION;
386 :
387 : double mass=0; bool test=false; int ntimes=10000;
388 0 : do
389 : {
390 0 : double y = EvtRandom::Flat(0,1);
391 :
392 : //generate mass
393 0 : mass = (_max-_min) * EvtRandom::Flat(0,1) + _min;
394 :
395 : //pdf calculate
396 : double f=0;
397 0 : switch(Vtype)
398 : {
399 0 : case VID::RHO : f = BreitWignerRelPDF(mass,MASS_RHO,WIDTH_RHO); break;
400 0 : case VID::OMEGA : f = BreitWignerRelPDF(mass,MASS_OMEGA,WIDTH_OMEGA); break;
401 0 : case VID::RHO_OMEGA_MIXING : f = RhoOmegaMixingPDF(mass,MASS_RHO,WIDTH_RHO,
402 0 : MASS_OMEGA,WIDTH_OMEGA); break;
403 0 : default : f = 1; break;
404 : }
405 :
406 0 : ntimes--;
407 0 : if (y<f) test=true;
408 0 : }while(ntimes && !test);
409 :
410 : //looping 10000 times
411 0 : if (ntimes==0)
412 : {
413 0 : report(Severity::Info,fname.c_str()) << "Tried accept/reject:10000"
414 0 : <<" times, and rejected all the times!"<<std::endl;
415 0 : report(Severity::Info,fname.c_str()) << "Is therefore accepting the last event!"<<std::endl;
416 0 : }
417 :
418 : //return V particle mass
419 : return mass;
420 : }
421 0 : }
422 :
423 :
424 :
425 :
426 :
427 : //************************************************************************
428 : //* *
429 : //* Class EvtLambda2PPiForLambdaB2LambdaV *
430 : //* *
431 : //************************************************************************
432 :
433 :
434 : //------------------------------------------------------------------------
435 : // Constructor
436 : //------------------------------------------------------------------------
437 0 : EvtLambda2PPiForLambdaB2LambdaV::EvtLambda2PPiForLambdaB2LambdaV()
438 0 : {
439 : //set facility name
440 0 : fname="EvtGen.EvtLambda2PPiForLambdaB2LambdaV";
441 0 : }
442 :
443 :
444 : //------------------------------------------------------------------------
445 : // Destructor
446 : //------------------------------------------------------------------------
447 : EvtLambda2PPiForLambdaB2LambdaV::~EvtLambda2PPiForLambdaB2LambdaV()
448 0 : {
449 0 : }
450 :
451 :
452 : //------------------------------------------------------------------------
453 : // Method 'getName'
454 : //------------------------------------------------------------------------
455 : std::string EvtLambda2PPiForLambdaB2LambdaV::getName()
456 : {
457 0 : return "LAMBDA2PPIFORLAMBDAB2LAMBDAV";
458 : }
459 :
460 :
461 : //------------------------------------------------------------------------
462 : // Method 'clone'
463 : //------------------------------------------------------------------------
464 : EvtDecayBase* EvtLambda2PPiForLambdaB2LambdaV::clone()
465 : {
466 0 : return new EvtLambda2PPiForLambdaB2LambdaV;
467 0 : }
468 :
469 :
470 : //------------------------------------------------------------------------
471 : // Method 'initProbMax'
472 : //------------------------------------------------------------------------
473 : void EvtLambda2PPiForLambdaB2LambdaV::initProbMax()
474 : {
475 : //maximum (case where D is real)
476 : double Max=0;
477 0 : if (A==0) Max=1;
478 0 : else if (C==0 || real(D)==0) Max=1+fabs(A*B);
479 0 : else if (B==0) Max=1+ EvtConst::pi/2.0*fabs(C*A*real(D));
480 : else
481 : {
482 0 : double theta_max = atan(- EvtConst::pi/2.0*C*real(D)/B);
483 0 : double max1 = 1 + fabs(A*B*cos(theta_max)
484 0 : - EvtConst::pi/2.0*C*A*real(D)*sin(theta_max));
485 0 : double max2 = 1 + fabs(A*B);
486 0 : if (max1>max2) Max=max1; else Max=max2;
487 : }
488 0 : report(Severity::Debug,fname.c_str())<<" PDF max value : "<<Max<<std::endl;
489 0 : setProbMax(Max);
490 0 : }
491 :
492 :
493 : //------------------------------------------------------------------------
494 : // Method 'init'
495 : //------------------------------------------------------------------------
496 : void EvtLambda2PPiForLambdaB2LambdaV::init()
497 : {
498 : bool antiparticle=false;
499 :
500 : //introduction
501 0 : report(Severity::Debug,fname.c_str())<<" ***********************************************************"<<std::endl;
502 0 : report(Severity::Debug,fname.c_str())<<" * Event Model Class : EvtLambda2PPiForLambdaB2LambdaV *"<<std::endl;
503 0 : report(Severity::Debug,fname.c_str())<<" ***********************************************************"<<std::endl;
504 :
505 : //check the number of arguments
506 0 : checkNArg(2);
507 :
508 : //check the number of daughters
509 0 : checkNDaug(2);
510 :
511 0 : const EvtId Id_mother = getParentId();
512 0 : const EvtId Id_daug1 = getDaug(0);
513 0 : const EvtId Id_daug2 = getDaug(1);
514 :
515 : //lambda identification
516 0 : if (Id_mother==EvtPDL::getId("Lambda0"))
517 : {
518 : antiparticle=false;
519 0 : }
520 0 : else if (Id_mother==EvtPDL::getId("anti-Lambda0"))
521 : {
522 : antiparticle=true;
523 : }
524 : else
525 : {
526 0 : report(Severity::Error,fname.c_str())<<" Mother is not a Lambda0 or an anti-Lambda0, but a "
527 0 : <<EvtPDL::name(Id_mother)<<std::endl;
528 0 : abort();
529 : }
530 :
531 : //proton identification
532 0 : if ( !(Id_daug1==EvtPDL::getId("p+") && !antiparticle ) && !(Id_daug1==EvtPDL::getId("anti-p-") && antiparticle) )
533 : {
534 0 : if (!antiparticle)
535 : {
536 0 : report(Severity::Error,fname.c_str()) << " Daughter1 is not a p+, but a "
537 0 : << EvtPDL::name(Id_daug1)<<std::endl;
538 0 : }
539 : else
540 0 : { report(Severity::Error,fname.c_str()) << " Daughter1 is not an anti-p-, but a "
541 0 : << EvtPDL::name(Id_daug1)<<std::endl;
542 : }
543 0 : abort();
544 : }
545 :
546 : //pion identification
547 0 : if ( !(Id_daug2==EvtPDL::getId("pi-") && !antiparticle ) && !(Id_daug2==EvtPDL::getId("pi+") && antiparticle) )
548 : {
549 0 : if (!antiparticle)
550 : {
551 0 : report(Severity::Error,fname.c_str()) << " Daughter2 is not a p-, but a "
552 0 : << EvtPDL::name(Id_daug1)<<std::endl;
553 0 : }
554 : else
555 0 : { report(Severity::Error,fname.c_str()) << " Daughter2 is not an p+, but a "
556 0 : << EvtPDL::name(Id_daug1)<<std::endl;
557 : }
558 0 : abort();
559 : }
560 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : Lambda0 -> p+ pi-"<<std::endl;
561 0 : else report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : Anti-Lambda0 -> anti-p- pi+"<<std::endl;
562 :
563 : //identification meson V
564 0 : if (getArg(1)==1)
565 : {
566 0 : Vtype=VID::JPSI;
567 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" From : Lambda_b0 -> Lambda J/psi"<<std::endl;
568 0 : else report(Severity::Debug,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda J/psi"<<std::endl;
569 : }
570 0 : else if (getArg(1)==2)
571 : {
572 0 : Vtype=VID::RHO;
573 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" From : Lambda_b0 -> Lambda rho0"<<std::endl;
574 0 : else report(Severity::Debug,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda rho0"<<std::endl;
575 : }
576 0 : else if (getArg(1)==3)
577 : {
578 0 : Vtype=VID::OMEGA;
579 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" From : Lambda_b0 -> Lambda omega"<<std::endl;
580 0 : else report(Severity::Debug,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda omega"<<std::endl;
581 : }
582 0 : else if (getArg(1)==4)
583 : {
584 0 : Vtype=VID::RHO_OMEGA_MIXING;
585 : }
586 : else
587 : {
588 0 : report(Severity::Error,fname.c_str()) << " Vtype " <<getArg(1)<<" is unknown"<<std::endl;
589 0 : if (!antiparticle) report(Severity::Debug,fname.c_str())<<" From : Lambda_b0 -> Lambda rho-omega-mixing"<<std::endl;
590 0 : else report(Severity::Debug,fname.c_str())<<" From : anti-Lambda_b0 -> anti-Lambda rho-omega-mixing"<<std::endl; abort();
591 : }
592 :
593 : //constants
594 0 : A = 0.642;
595 0 : C = (double) getArg(0);
596 0 : switch(Vtype)
597 : {
598 0 : case VID::JPSI : B = -0.167; D = EvtComplex(0.25,0.0); break;
599 : case VID::RHO :
600 : case VID::OMEGA :
601 0 : case VID::RHO_OMEGA_MIXING : B = -0.21; D = EvtComplex(0.31,0.0); break;
602 0 : default : B = 0; D = EvtComplex(0,0); break;
603 : }
604 :
605 :
606 0 : report(Severity::Debug,fname.c_str())<<" Lambda decay parameters : "<<std::endl;
607 0 : report(Severity::Debug,fname.c_str())<<" - proton asymmetry A = "<<A<<std::endl;
608 0 : report(Severity::Debug,fname.c_str())<<" - lambda polarisation B = "<<B<<std::endl;
609 0 : report(Severity::Debug,fname.c_str())<<" - lambdaB polarisation C = "<<C<<std::endl;
610 0 : report(Severity::Debug,fname.c_str())<<" - lambda density matrix rho+- D = "<<D<<std::endl;
611 0 : }
612 :
613 :
614 : //------------------------------------------------------------------------
615 : // Method 'decay'
616 : //------------------------------------------------------------------------
617 : void EvtLambda2PPiForLambdaB2LambdaV::decay( EvtParticle *lambda )
618 : {
619 0 : lambda->makeDaughters(getNDaug(),getDaugs());
620 :
621 : //get proton and pion particles
622 0 : EvtParticle * proton = lambda->getDaug(0);
623 0 : EvtParticle * pion = lambda->getDaug(1);
624 :
625 : //get resonance masses
626 : // - LAMBDA -> mass given by EvtLambdaB2LambdaV class
627 : // - PROTON & PION -> nominal mass
628 0 : double MASS_LAMBDA = lambda->mass();
629 0 : double MASS_PROTON = EvtPDL::getMeanMass(EvtPDL::getId("p+"));
630 0 : double MASS_PION = EvtPDL::getMeanMass(EvtPDL::getId("pi-"));
631 :
632 : //generate random angles
633 0 : double phi = EvtRandom::Flat(0,2*EvtConst::pi);
634 0 : double theta = acos( EvtRandom::Flat(-1,+1));
635 0 : report(Severity::Debug,fname.c_str())<<" Angular angles : theta = "<<theta<<" ; phi = "<<phi<<std::endl;
636 :
637 : //computate resonance quadrivectors
638 0 : double E_proton = (MASS_LAMBDA*MASS_LAMBDA + MASS_PROTON*MASS_PROTON - MASS_PION*MASS_PION)
639 0 : /(2*MASS_LAMBDA);
640 0 : double E_pion = (MASS_LAMBDA*MASS_LAMBDA + MASS_PION*MASS_PION - MASS_PROTON*MASS_PROTON)
641 0 : /(2*MASS_LAMBDA);
642 0 : double P = sqrt(E_proton*E_proton-proton->mass()*proton->mass());
643 :
644 0 : EvtVector4R P_lambda=lambda->getP4();
645 0 : EvtParticle *Mother_lambda=lambda->getParent();
646 0 : EvtVector4R lambdab=Mother_lambda->getP4();
647 :
648 :
649 :
650 0 : double E_lambda =P_lambda.get(0);
651 0 : double px_M =lambdab.get(1);
652 0 : double py_M =lambdab.get(2);
653 0 : double pz_M =lambdab.get(3);
654 0 : double E_M =lambdab.get(0);
655 :
656 0 : EvtVector4R q_lambdab2 (E_M,
657 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(px_M))+(py_M*(py_M)))),
658 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-((py_M)*(px_M))+(px_M*(py_M)))),
659 : (pz_M));
660 :
661 0 : EvtVector4R q_lambdab3 (E_M,
662 0 : q_lambdab2.get(3),
663 0 : q_lambdab2.get(1),
664 0 : q_lambdab2.get(2));
665 :
666 :
667 :
668 0 : EvtVector4R q_lambda1 (E_lambda,
669 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(P_lambda.get(1))) + (py_M*(P_lambda.get(2))))),
670 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-(py_M*(P_lambda.get(1))) + (px_M*(P_lambda.get(2))))),
671 0 : P_lambda.get(3));
672 :
673 0 : EvtVector4R q_lambda2 (E_lambda,
674 0 : q_lambda1.get(3),
675 0 : q_lambda1.get(1),
676 0 : q_lambda1.get(2));
677 :
678 :
679 :
680 :
681 :
682 0 : double px=q_lambda2.get(1);
683 0 : double py=q_lambda2.get(2);
684 0 : double pz=q_lambda2.get(3);
685 :
686 :
687 :
688 :
689 0 : EvtVector4R q_lambda4 (q_lambda2.get(0),
690 0 : ((1/(sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2) + pow(q_lambda2.get(3),2))))* (1/(sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2))))*((q_lambda2.get(1))*(q_lambda2.get(1))*(q_lambda2.get(3))+((q_lambda2.get(2))*(q_lambda2.get(2))*(q_lambda2.get(3))) - ((q_lambda2.get(3))*(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2))))),
691 0 : ((((q_lambda2.get(2))*(q_lambda2.get(1)))-((q_lambda2.get(1))*(q_lambda2.get(2))))/(sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2)))),
692 0 : (((1/sqrt(pow(q_lambda2.get(1),2) + pow(q_lambda2.get(2),2) + pow(q_lambda2.get(3),2)))*( ((q_lambda2.get(1))*(q_lambda2.get(1))) +((q_lambda2.get(2))*(q_lambda2.get(2))) + ((q_lambda2.get(3))*(q_lambda2.get(3)))))) );
693 :
694 :
695 :
696 :
697 0 : EvtVector4R q_proton1 (E_proton,
698 0 : P*sin(theta)*cos(phi),
699 0 : P*sin(theta)*sin(phi),
700 0 : P*cos(theta));
701 0 : EvtVector4R q_pion1 (E_pion,
702 0 : -P*sin(theta)*cos(phi),
703 0 : -P*sin(theta)*sin(phi),
704 0 : -P*cos(theta));
705 :
706 :
707 0 : EvtVector4R q_proton3 (E_proton,
708 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_proton1.get(1))*(px)*(pz) - ((q_proton1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_proton1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
709 0 : (((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))* (1/(sqrt(pow(px,2) + pow(py,2))))*(((q_proton1.get(1)))*(py)*(pz) + ((q_proton1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_proton1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
710 0 : (((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))*((-(q_proton1.get(1)))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_proton1.get(3))*(pz)))));
711 :
712 0 : EvtVector4R q_pion3 (E_pion,
713 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_pion1.get(1))*(px)*(pz) - ((q_pion1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_pion1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
714 0 : (((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_pion1.get(1))*(py)*(pz) + ((q_pion1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) + (((q_pion1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
715 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))*((-(q_pion1.get(1)))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_pion1.get(3))*(pz)))));
716 :
717 0 : EvtVector4R q_proton5 (q_proton3.get(0),
718 0 : (q_proton3.get(2)),
719 0 : (q_proton3.get(3)),
720 0 : (q_proton3.get(1)));
721 :
722 0 : EvtVector4R q_pion5 (q_pion3.get(0),
723 0 : (q_pion3.get(2)),
724 0 : (q_pion3.get(3)),
725 0 : (q_pion3.get(1)));
726 :
727 0 : EvtVector4R q_proton (q_proton5.get(0),
728 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_proton5.get(1)))-(py_M*(q_proton5.get(2))))),
729 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_proton5.get(1)))+(px_M*(q_proton5.get(2))))),
730 0 : (q_proton5.get(3)));
731 :
732 :
733 0 : EvtVector4R q_pion (q_pion5.get(0),
734 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_pion5.get(1)))-(py_M*(q_pion5.get(2))))),
735 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_pion5.get(1)))+(px_M*(q_pion5.get(2))))),
736 0 : (q_pion5.get(3)));
737 :
738 0 : report(Severity::Info,fname.c_str())<<" Lambdab px: "<<px_M<<std::endl;
739 0 : report(Severity::Info,fname.c_str())<<" Lambdab py: "<<py_M<<std::endl;
740 0 : report(Severity::Info,fname.c_str())<<" Lambdab pz: "<<pz_M<<std::endl;
741 0 : report(Severity::Info,fname.c_str())<<" Lambdab E: "<<E_M<<std::endl;
742 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 px: "<<q_lambdab2.get(1)<<std::endl;
743 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 py: "<<q_lambdab2.get(2)<<std::endl;
744 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 pz: "<<q_lambdab2.get(3)<<std::endl;
745 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 E: "<<q_lambdab2.get(0)<<std::endl;
746 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 px: "<<q_lambdab3.get(1)<<std::endl;
747 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 py: "<<q_lambdab3.get(2)<<std::endl;
748 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 pz: "<<q_lambdab3.get(3)<<std::endl;
749 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 E: "<<q_lambdab3.get(0)<<std::endl;
750 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 px: "<<P_lambda.get(1)<<std::endl;
751 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 py: "<<P_lambda.get(2)<<std::endl;
752 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 pz: "<<P_lambda.get(3)<<std::endl;
753 0 : report(Severity::Info,fname.c_str())<<" Lambda 0 E: "<<P_lambda.get(0)<<std::endl;
754 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 px: "<<q_lambda1.get(1)<<std::endl;
755 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 py: "<<q_lambda1.get(2)<<std::endl;
756 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 pz: "<<q_lambda1.get(3)<<std::endl;
757 0 : report(Severity::Info,fname.c_str())<<" Lambda 1 E: "<<q_lambda1.get(0)<<std::endl;
758 0 : report(Severity::Info,fname.c_str())<<" Lambda 2 px: "<<q_lambda2.get(1)<<std::endl;
759 0 : report(Severity::Info,fname.c_str())<<" Lambda 2 py: "<<q_lambda2.get(2)<<std::endl;
760 0 : report(Severity::Info,fname.c_str())<<" Lambda 2 pz: "<<q_lambda2.get(3)<<std::endl;
761 0 : report(Severity::Info,fname.c_str())<<" Lambda 2 E: "<<q_lambda2.get(0)<<std::endl;
762 :
763 0 : report(Severity::Info,fname.c_str())<<" Lambda px: "<<px<<std::endl;
764 0 : report(Severity::Info,fname.c_str())<<" Lambda py: "<<py<<std::endl;
765 0 : report(Severity::Info,fname.c_str())<<" Lambda pz: "<<pz<<std::endl;
766 :
767 0 : report(Severity::Info,fname.c_str())<<" pion 1 px: "<<q_pion1.get(1)<<std::endl;
768 0 : report(Severity::Info,fname.c_str())<<" pion 1 py: "<<q_pion1.get(2)<<std::endl;
769 0 : report(Severity::Info,fname.c_str())<<" pion 1 pz: "<<q_pion1.get(3)<<std::endl;
770 0 : report(Severity::Info,fname.c_str())<<" pion 1 E: "<<q_pion1.get(0)<<std::endl;
771 :
772 0 : report(Severity::Info,fname.c_str())<<" pion 3 px: "<<q_pion3.get(1)<<std::endl;
773 0 : report(Severity::Info,fname.c_str())<<" pion 3 px: "<<q_pion3.get(1)<<std::endl;
774 0 : report(Severity::Info,fname.c_str())<<" pion 3 py: "<<q_pion3.get(2)<<std::endl;
775 0 : report(Severity::Info,fname.c_str())<<" pion 3 pz: "<<q_pion3.get(3)<<std::endl;
776 0 : report(Severity::Info,fname.c_str())<<" pion 3 E: "<<q_pion3.get(0)<<std::endl;
777 :
778 0 : report(Severity::Info,fname.c_str())<<" pion 5 px: "<<q_pion5.get(1)<<std::endl;
779 0 : report(Severity::Info,fname.c_str())<<" pion 5 py: "<<q_pion5.get(2)<<std::endl;
780 0 : report(Severity::Info,fname.c_str())<<" pion 5 pz: "<<q_pion5.get(3)<<std::endl;
781 0 : report(Severity::Info,fname.c_str())<<" pion 5 E: "<<q_pion5.get(0)<<std::endl;
782 :
783 :
784 :
785 0 : report(Severity::Info,fname.c_str())<<" proton 1 px: "<<q_proton1.get(1)<<std::endl;
786 0 : report(Severity::Info,fname.c_str())<<" proton 1 py: "<<q_proton1.get(2)<<std::endl;
787 0 : report(Severity::Info,fname.c_str())<<" proton 1 pz: "<<q_proton1.get(3)<<std::endl;
788 0 : report(Severity::Info,fname.c_str())<<" proton 1 E: "<<q_proton1.get(0)<<std::endl;
789 :
790 0 : report(Severity::Info,fname.c_str())<<" proton 3 px: "<<q_proton3.get(1)<<std::endl;
791 0 : report(Severity::Info,fname.c_str())<<" proton 3 py: "<<q_proton3.get(2)<<std::endl;
792 0 : report(Severity::Info,fname.c_str())<<" proton 3 pz: "<<q_proton3.get(3)<<std::endl;
793 0 : report(Severity::Info,fname.c_str())<<" proton 3 E: "<<q_proton3.get(0)<<std::endl;
794 :
795 0 : report(Severity::Info,fname.c_str())<<" proton 5 px: "<<q_proton5.get(1)<<std::endl;
796 0 : report(Severity::Info,fname.c_str())<<" proton 5 py: "<<q_proton5.get(2)<<std::endl;
797 0 : report(Severity::Info,fname.c_str())<<" proton 5 pz: "<<q_proton5.get(3)<<std::endl;
798 0 : report(Severity::Info,fname.c_str())<<" proton 5 E: "<<q_proton5.get(0)<<std::endl;
799 :
800 :
801 0 : report(Severity::Info,fname.c_str())<<" proton px: "<<q_proton.get(1)<<std::endl;
802 0 : report(Severity::Info,fname.c_str())<<" proton py: "<<q_proton.get(2)<<std::endl;
803 0 : report(Severity::Info,fname.c_str())<<"proton pz: "<<q_proton.get(3)<<std::endl;
804 0 : report(Severity::Info,fname.c_str())<<" pion px: "<<q_pion.get(1)<<std::endl;
805 0 : report(Severity::Info,fname.c_str())<<" pion py: "<<q_pion.get(2)<<std::endl;
806 0 : report(Severity::Info,fname.c_str())<<" pion pz: "<<q_pion.get(3)<<std::endl;
807 :
808 :
809 :
810 :
811 :
812 : ;
813 :
814 :
815 :
816 :
817 :
818 :
819 : ///////////*******NEW********//////////////////////
820 :
821 : //set quadrivectors to particles
822 0 : proton->init(getDaugs()[0],q_proton);
823 0 : pion ->init(getDaugs()[1],q_pion );
824 :
825 : //computate pdf
826 : //double pdf = 1 + A*B*cos(theta) - EvtConst::pi/2.0*C*A*real(D*EvtComplex(cos(phi),sin(phi)))*sin(theta);
827 0 : double pdf = 1 + A*B*cos(theta) + 2*A*real(D*EvtComplex(cos(phi),sin(phi)))*sin(theta);
828 0 : report(Severity::Debug,fname.c_str())<<" Lambda decay pdf value : "<<pdf<<std::endl;
829 : //set probability
830 0 : setProb(pdf);
831 :
832 : return;
833 0 : }
834 :
835 :
836 :
837 :
838 : //************************************************************************
839 : //* *
840 : //* Class EvtLambda2PPiForLambdaB2LambdaV *
841 : //* *
842 : //************************************************************************
843 :
844 :
845 : //------------------------------------------------------------------------
846 : // Constructor
847 : //------------------------------------------------------------------------
848 0 : EvtV2VpVmForLambdaB2LambdaV::EvtV2VpVmForLambdaB2LambdaV()
849 0 : {
850 : //set facility name
851 0 : fname="EvtGen.EvtV2VpVmForLambdaB2LambdaV";
852 0 : }
853 :
854 :
855 : //------------------------------------------------------------------------
856 : // Destructor
857 : //------------------------------------------------------------------------
858 : EvtV2VpVmForLambdaB2LambdaV::~EvtV2VpVmForLambdaB2LambdaV()
859 0 : {}
860 :
861 :
862 : //------------------------------------------------------------------------
863 : // Method 'getName'
864 : //------------------------------------------------------------------------
865 : std::string EvtV2VpVmForLambdaB2LambdaV::getName()
866 : {
867 0 : return "V2VPVMFORLAMBDAB2LAMBDAV";
868 : }
869 :
870 : //------------------------------------------------------------------------
871 : // Method 'clone'
872 : //------------------------------------------------------------------------
873 : EvtDecayBase* EvtV2VpVmForLambdaB2LambdaV::clone()
874 : {
875 0 : return new EvtV2VpVmForLambdaB2LambdaV;
876 0 : }
877 :
878 :
879 : //------------------------------------------------------------------------
880 : // Method 'initProbMax'
881 : //------------------------------------------------------------------------
882 : void EvtV2VpVmForLambdaB2LambdaV::initProbMax()
883 : {
884 : //maximum
885 : double Max = 0;
886 0 : if (Vtype==VID::JPSI)
887 : {
888 0 : if ((1-3*A)>0) Max=2*(1-A);
889 0 : else Max=1+A;
890 : }
891 : else
892 : {
893 0 : if ((3*A-1)>=0) Max=2*A;
894 0 : else Max=1-A;
895 : }
896 :
897 0 : report(Severity::Debug,fname.c_str())<<" PDF max value : "<<Max<<std::endl;
898 0 : setProbMax(Max);
899 0 : }
900 :
901 :
902 : //------------------------------------------------------------------------
903 : // Method 'init'
904 : //------------------------------------------------------------------------
905 : void EvtV2VpVmForLambdaB2LambdaV::init()
906 : {
907 : //introduction
908 0 : report(Severity::Debug,fname.c_str())<<" ***********************************************************"<<std::endl;
909 0 : report(Severity::Debug,fname.c_str())<<" * Event Model Class : EvtV2VpVmForLambdaB2LambdaV *"<<std::endl;
910 0 : report(Severity::Debug,fname.c_str())<<" ***********************************************************"<<std::endl;
911 :
912 : //check the number of arguments
913 0 : checkNArg(2);
914 :
915 : //check the number of daughters
916 0 : checkNDaug(2);
917 :
918 0 : const EvtId Id_mother = getParentId();
919 0 : const EvtId Id_daug1 = getDaug(0);
920 0 : const EvtId Id_daug2 = getDaug(1);
921 :
922 : //identification meson V
923 0 : if (getArg(1)==1) Vtype=VID::JPSI;
924 0 : else if (getArg(1)==2) Vtype=VID::RHO;
925 0 : else if (getArg(1)==3) Vtype=VID::OMEGA;
926 0 : else if (getArg(1)==4) Vtype=VID::RHO_OMEGA_MIXING;
927 : else
928 : {
929 0 : report(Severity::Error,fname.c_str()) << " Vtype " <<getArg(1)<<" is unknown"<<std::endl;
930 0 : abort();
931 : }
932 :
933 : //vector meson V
934 0 : if (Id_mother==EvtPDL::getId("J/psi") && Vtype==VID::JPSI)
935 : {
936 : }
937 0 : else if (Id_mother==EvtPDL::getId("omega") && Vtype==VID::OMEGA)
938 : {
939 : }
940 0 : else if (Id_mother==EvtPDL::getId("rho0") && Vtype==VID::RHO)
941 : {
942 : }
943 0 : else if ((Id_mother==EvtPDL::getId("rho0") || Id_mother==EvtPDL::getId("omega")) && Vtype==VID::RHO_OMEGA_MIXING)
944 : {
945 : }
946 : else
947 : {
948 0 : report(Severity::Error,fname.c_str())<<" Mother is not a J/psi, phi or rho0 but a "
949 0 : <<EvtPDL::name(Id_mother)<<std::endl;
950 0 : abort();
951 : }
952 :
953 : //daughters for each V possibility
954 0 : switch(Vtype)
955 : {
956 : case VID::JPSI :
957 0 : if (Id_daug1!=EvtPDL::getId("mu+"))
958 : {
959 0 : report(Severity::Error,fname.c_str()) << " Daughter1 is not a mu+, but a "
960 0 : << EvtPDL::name(Id_daug1)<<std::endl;
961 0 : abort();
962 : }
963 0 : if (Id_daug2!=EvtPDL::getId("mu-"))
964 : {
965 0 : report(Severity::Error,fname.c_str()) << " Daughter2 is not a mu-, but a "
966 0 : << EvtPDL::name(Id_daug2)<<std::endl;
967 0 : abort();
968 : }
969 0 : report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : J/psi -> mu+ mu-"<<std::endl;
970 0 : break;
971 :
972 : case VID::RHO :
973 : case VID::OMEGA :
974 : case VID::RHO_OMEGA_MIXING :
975 0 : if (Id_daug1!=EvtPDL::getId("pi+"))
976 : {
977 0 : report(Severity::Error,fname.c_str()) << " Daughter1 is not a pi+, but a "
978 0 : << EvtPDL::name(Id_daug1)<<std::endl;
979 0 : abort();
980 : }
981 0 : if (Id_daug2!=EvtPDL::getId("pi-"))
982 : {
983 0 : report(Severity::Error,fname.c_str()) << " Daughter2 is not a pi-, but a "
984 0 : << EvtPDL::name(Id_daug2)<<std::endl;
985 0 : abort();
986 : }
987 0 : if (Vtype==VID::RHO) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : rho0 -> pi+ pi-"<<std::endl;
988 0 : if (Vtype==VID::OMEGA) report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : omega -> pi+ pi-"<<std::endl;
989 0 : if (Vtype==VID::RHO_OMEGA_MIXING)
990 0 : report(Severity::Debug,fname.c_str())<<" Decay mode successfully initialized : rho-omega mixing -> pi+ pi-"<<std::endl; break;
991 :
992 : default :
993 0 : report(Severity::Error,fname.c_str()) << "No decay mode chosen ! "<<std::endl;
994 0 : abort();
995 : break;
996 : }
997 :
998 : //fix dynamics parameters
999 0 : switch(Vtype)
1000 : {
1001 0 : case VID::JPSI : A = 0.66; break;
1002 : case VID::RHO :
1003 : case VID::OMEGA :
1004 0 : case VID::RHO_OMEGA_MIXING : A = 0.79; break;
1005 0 : default : A = 0; break;
1006 : }
1007 :
1008 0 : report(Severity::Debug,fname.c_str())<<" V decay parameters : "<<std::endl;
1009 0 : report(Severity::Debug,fname.c_str())<<" - V density matrix rho00 A = "<<A<<std::endl;
1010 :
1011 :
1012 0 : }
1013 :
1014 : //------------------------------------------------------------------------
1015 : // Method 'decay'
1016 : //------------------------------------------------------------------------
1017 : void EvtV2VpVmForLambdaB2LambdaV::decay( EvtParticle *V )
1018 : {
1019 0 : V->makeDaughters(getNDaug(),getDaugs());
1020 :
1021 : //get Vp and Vm particles
1022 0 : EvtParticle * Vp = V->getDaug(0);
1023 0 : EvtParticle * Vm = V->getDaug(1);
1024 :
1025 : //get resonance masses
1026 : // - V -> mass given by EvtLambdaB2LambdaV class
1027 : // - Vp & Vm -> nominal mass
1028 0 : double MASS_V = V->mass();
1029 : double MASS_VM = 0;
1030 0 : switch(Vtype)
1031 : {
1032 0 : case VID::JPSI : MASS_VM=EvtPDL::getMeanMass(EvtPDL::getId("mu-")); break;
1033 : case VID::RHO :
1034 : case VID::OMEGA :
1035 0 : case VID::RHO_OMEGA_MIXING : MASS_VM=EvtPDL::getMeanMass(EvtPDL::getId("pi-")); break;
1036 0 : default : MASS_VM=0; break;
1037 : }
1038 : double MASS_VP = MASS_VM;
1039 :
1040 : //generate random angles
1041 0 : double phi = EvtRandom::Flat(0,2*EvtConst::pi);
1042 0 : double theta = acos( EvtRandom::Flat(-1,+1));
1043 0 : report(Severity::Debug,fname.c_str())<<" Angular angles : theta = "<<theta<<" ; phi = "<<phi<<std::endl;
1044 :
1045 : //computate resonance quadrivectors
1046 0 : double E_Vp = (MASS_V*MASS_V + MASS_VP*MASS_VP - MASS_VM*MASS_VM)
1047 0 : /(2*MASS_V);
1048 : double E_Vm = (MASS_V*MASS_V + MASS_VM*MASS_VM - MASS_VP*MASS_VP)
1049 : /(2*MASS_V);
1050 0 : double P = sqrt(E_Vp*E_Vp-Vp->mass()*Vp->mass());
1051 :
1052 0 : EvtVector4R P_V=V->getP4();
1053 0 : EvtParticle *Mother_V=V->getParent();
1054 0 : EvtVector4R lambdab=Mother_V->getP4();
1055 :
1056 :
1057 0 : double E_V=(P_V.get(0));
1058 0 : double px_M=lambdab.get(1);
1059 0 : double py_M=lambdab.get(2);
1060 0 : double pz_M=lambdab.get(3);
1061 0 : double E_M=lambdab.get(0);
1062 :
1063 0 : EvtVector4R q_lambdab2 (E_M,
1064 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(px_M))+(py_M*(py_M)))),
1065 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-((py_M)*(px_M))+(px_M*(py_M)))),
1066 : (pz_M));
1067 :
1068 0 : EvtVector4R q_lambdab3 (E_M,
1069 0 : q_lambdab2.get(3),
1070 0 : q_lambdab2.get(1),
1071 0 : q_lambdab2.get(2));
1072 :
1073 :
1074 0 : EvtVector4R q_V1 (E_V,
1075 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(P_V.get(1))) + (py_M*(P_V.get(2))))),
1076 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*(-(py_M*(P_V.get(1))) + (px_M*(P_V.get(2))))),
1077 0 : P_V.get(3));
1078 :
1079 0 : EvtVector4R q_V2 (E_V,
1080 0 : q_V1.get(3),
1081 0 : q_V1.get(1),
1082 0 : q_V1.get(2));
1083 :
1084 :
1085 :
1086 0 : double px= -(q_V2.get(1));
1087 0 : double py=-(q_V2.get(2));
1088 0 : double pz=-(q_V2.get(3));
1089 :
1090 :
1091 :
1092 :
1093 0 : EvtVector4R q_V4 (q_V2.get(0),
1094 0 : ((1/(sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2) + pow(q_V2.get(3),2))))* (1/(sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2))))*((q_V2.get(1))*(q_V2.get(1))*(q_V2.get(3))+((q_V2.get(2))*(q_V2.get(2))*(q_V2.get(3))) - ((q_V2.get(3))*(pow(q_V2.get(1),2) + pow(q_V2.get(2),2))))),
1095 0 : ((((q_V2.get(2))*(q_V2.get(1)))-((q_V2.get(1))*(q_V2.get(2))))/(sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2)))),
1096 0 : (((1/sqrt(pow(q_V2.get(1),2) + pow(q_V2.get(2),2) + pow(q_V2.get(3),2)))*( ((q_V2.get(1))*(q_V2.get(1))) +((q_V2.get(2))*(q_V2.get(2))) + ((q_V2.get(3))*(q_V2.get(3)))))) );
1097 :
1098 :
1099 :
1100 0 : EvtVector4R q_Vp1 (E_Vp,
1101 0 : P*sin(theta)*cos(phi),
1102 0 : P*sin(theta)*sin(phi),
1103 0 : P*cos(theta));
1104 0 : EvtVector4R q_Vm1 (E_Vm,
1105 0 : -P*sin(theta)*cos(phi),
1106 0 : -P*sin(theta)*sin(phi),
1107 0 : -P*cos(theta));
1108 :
1109 0 : EvtVector4R q_Vp3 (q_Vp1.get(0),
1110 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_Vp1.get(1))*(px)*(pz)+((q_Vp1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vp1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
1111 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*(((q_Vp1.get(1)))*(py)*(pz) - ((q_Vp1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vp1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
1112 0 : ((-(1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))*((q_Vp1.get(1))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_Vp1.get(3))*(pz)))));
1113 :
1114 0 : EvtVector4R q_Vm3 (q_Vm1.get(0),
1115 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*((q_Vm1.get(1))*(px)*(pz)+((q_Vm1.get(2))*(py)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vm1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(px)))),
1116 0 : ((1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))* (1/(sqrt(pow(px,2) + pow(py,2))))*(((q_Vm1.get(1)))*(py)*(pz) - ((q_Vm1.get(2))*(px)*((sqrt(pow(px,2) + pow(py,2) + pow(pz,2))))) - (((q_Vm1.get(3)))*(sqrt(pow(px,2) + pow(py,2)))*(py)))) ,
1117 0 : ((-(1/(sqrt(pow(px,2) + pow(py,2) + pow(pz,2)))))*((q_Vm1.get(1))*((sqrt(pow(px,2) + pow(py,2)))) + ((q_Vm1.get(3))*(pz)))));
1118 :
1119 :
1120 :
1121 :
1122 :
1123 :
1124 0 : EvtVector4R q_Vp5 (q_Vp3.get(0),
1125 0 : (q_Vp3.get(2)),
1126 0 : (q_Vp3.get(3)),
1127 0 : (q_Vp3.get(1)));
1128 :
1129 0 : EvtVector4R q_Vm5 (q_Vm3.get(0),
1130 0 : (q_Vm3.get(2)),
1131 0 : (q_Vm3.get(3)),
1132 0 : (q_Vm3.get(1)));
1133 :
1134 :
1135 0 : EvtVector4R q_Vp (q_Vp5.get(0),
1136 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_Vp5.get(1)))-(py_M*(q_Vp5.get(2))))),
1137 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_Vp5.get(1)))+(px_M*(q_Vp5.get(2))))),
1138 0 : (q_Vp5.get(3)));
1139 :
1140 :
1141 0 : EvtVector4R q_Vm (q_Vm5.get(0),
1142 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((px_M*(q_Vm5.get(1)))-(py_M*(q_Vm5.get(2))))),
1143 0 : ((1/(sqrt(pow(px_M,2)+pow(py_M,2))))*((py_M*(q_Vm5.get(1)))+(px_M*(q_Vm5.get(2))))),
1144 0 : (q_Vm5.get(3)));
1145 :
1146 0 : report(Severity::Info,fname.c_str())<<" Lambdab px: "<<px_M<<std::endl;
1147 0 : report(Severity::Info,fname.c_str())<<" Lambdab py: "<<py_M<<std::endl;
1148 0 : report(Severity::Info,fname.c_str())<<" Lambdab pz: "<<pz_M<<std::endl;
1149 0 : report(Severity::Info,fname.c_str())<<" Lambdab E: "<<E_M<<std::endl;
1150 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 px: "<<q_lambdab2.get(1)<<std::endl;
1151 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 py: "<<q_lambdab2.get(2)<<std::endl;
1152 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 pz: "<<q_lambdab2.get(3)<<std::endl;
1153 0 : report(Severity::Info,fname.c_str())<<" Lambdab2 E: "<<q_lambdab2.get(0)<<std::endl;
1154 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 px: "<<q_lambdab3.get(1)<<std::endl;
1155 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 py: "<<q_lambdab3.get(2)<<std::endl;
1156 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 pz: "<<q_lambdab3.get(3)<<std::endl;
1157 0 : report(Severity::Info,fname.c_str())<<" Lambdab3 E: "<<q_lambdab3.get(0)<<std::endl;
1158 0 : report(Severity::Info,fname.c_str())<<" V 0 px: "<<P_V.get(1)<<std::endl;
1159 0 : report(Severity::Info,fname.c_str())<<" V 0 py: "<<P_V.get(2)<<std::endl;
1160 0 : report(Severity::Info,fname.c_str())<<" V 0 pz: "<<P_V.get(3)<<std::endl;
1161 0 : report(Severity::Info,fname.c_str())<<" V 0 E: "<<P_V.get(0)<<std::endl;
1162 0 : report(Severity::Info,fname.c_str())<<" V 1 px: "<<q_V1.get(1)<<std::endl;
1163 0 : report(Severity::Info,fname.c_str())<<" V 1 py: "<<q_V1.get(2)<<std::endl;
1164 0 : report(Severity::Info,fname.c_str())<<" V 1 pz: "<<q_V1.get(3)<<std::endl;
1165 0 : report(Severity::Info,fname.c_str())<<" V 1 E: "<<q_V1.get(0)<<std::endl;
1166 0 : report(Severity::Info,fname.c_str())<<" V 2 px: "<<q_V2.get(1)<<std::endl;
1167 0 : report(Severity::Info,fname.c_str())<<" V 2 py: "<<q_V2.get(2)<<std::endl;
1168 0 : report(Severity::Info,fname.c_str())<<" V 2 pz: "<<q_V2.get(3)<<std::endl;
1169 0 : report(Severity::Info,fname.c_str())<<" V 2 E: "<<q_V2.get(0)<<std::endl;
1170 0 : report(Severity::Info,fname.c_str())<<" V px: "<<px<<std::endl;
1171 0 : report(Severity::Info,fname.c_str())<<" V py: "<<py<<std::endl;
1172 0 : report(Severity::Info,fname.c_str())<<" V pz: "<<pz<<std::endl;
1173 0 : report(Severity::Info,fname.c_str())<<" Vm 1 px: "<<q_Vm1.get(1)<<std::endl;
1174 0 : report(Severity::Info,fname.c_str())<<" Vm 1 py: "<<q_Vm1.get(2)<<std::endl;
1175 0 : report(Severity::Info,fname.c_str())<<" Vm 1 pz: "<<q_Vm1.get(3)<<std::endl;
1176 0 : report(Severity::Info,fname.c_str())<<" Vm 1 E: "<<q_Vm1.get(0)<<std::endl;
1177 0 : report(Severity::Info,fname.c_str())<<" Vm 3 px: "<<q_Vm3.get(1)<<std::endl;
1178 0 : report(Severity::Info,fname.c_str())<<" Vm 3 px: "<<q_Vm3.get(1)<<std::endl;
1179 0 : report(Severity::Info,fname.c_str())<<" Vm 3 py: "<<q_Vm3.get(2)<<std::endl;
1180 0 : report(Severity::Info,fname.c_str())<<" Vm 3 pz: "<<q_Vm3.get(3)<<std::endl;
1181 0 : report(Severity::Info,fname.c_str())<<" Vm 3 E: "<<q_Vm3.get(0)<<std::endl;
1182 0 : report(Severity::Info,fname.c_str())<<" Vm 5 px: "<<q_Vm5.get(1)<<std::endl;
1183 0 : report(Severity::Info,fname.c_str())<<" Vm 5 py: "<<q_Vm5.get(2)<<std::endl;
1184 0 : report(Severity::Info,fname.c_str())<<" Vm 5 pz: "<<q_Vm5.get(3)<<std::endl;
1185 0 : report(Severity::Info,fname.c_str())<<" Vm 5 E: "<<q_Vm5.get(0)<<std::endl;
1186 0 : report(Severity::Info,fname.c_str())<<" Vp 1 px: "<<q_Vp1.get(1)<<std::endl;
1187 0 : report(Severity::Info,fname.c_str())<<" Vp 1 py: "<<q_Vp1.get(2)<<std::endl;
1188 0 : report(Severity::Info,fname.c_str())<<" Vp 1 pz: "<<q_Vp1.get(3)<<std::endl;
1189 0 : report(Severity::Info,fname.c_str())<<" Vp 1 E: "<<q_Vp1.get(0)<<std::endl;
1190 0 : report(Severity::Info,fname.c_str())<<" Vp 3 px: "<<q_Vp3.get(1)<<std::endl;
1191 0 : report(Severity::Info,fname.c_str())<<" Vp 3 py: "<<q_Vp3.get(2)<<std::endl;
1192 0 : report(Severity::Info,fname.c_str())<<" Vp 3 pz: "<<q_Vp3.get(3)<<std::endl;
1193 0 : report(Severity::Info,fname.c_str())<<" Vp 3 E: "<<q_Vp3.get(0)<<std::endl;
1194 0 : report(Severity::Info,fname.c_str())<<" Vp 5 px: "<<q_Vp5.get(1)<<std::endl;
1195 0 : report(Severity::Info,fname.c_str())<<" Vp 5 py: "<<q_Vp5.get(2)<<std::endl;
1196 0 : report(Severity::Info,fname.c_str())<<" Vp 5 pz: "<<q_Vp5.get(3)<<std::endl;
1197 0 : report(Severity::Info,fname.c_str())<<" Vp 5 E: "<<q_Vp5.get(0)<<std::endl;
1198 0 : report(Severity::Info,fname.c_str())<<" Vp px: "<<q_Vp.get(1)<<std::endl;
1199 0 : report(Severity::Info,fname.c_str())<<" Vp py: "<<q_Vp.get(2)<<std::endl;
1200 0 : report(Severity::Info,fname.c_str())<<"Vp pz: "<<q_Vp.get(3)<<std::endl;
1201 0 : report(Severity::Info,fname.c_str())<<" Vm px: "<<q_Vm.get(1)<<std::endl;
1202 0 : report(Severity::Info,fname.c_str())<<" Vm py: "<<q_Vm.get(2)<<std::endl;
1203 0 : report(Severity::Info,fname.c_str())<<" Vm pz: "<<q_Vm.get(3)<<std::endl;
1204 :
1205 : //set quadrivectors to particles
1206 0 : Vp->init(getDaugs()[0],q_Vp);
1207 0 : Vm->init(getDaugs()[1],q_Vm);
1208 :
1209 : //computate pdf
1210 : double pdf = 0;
1211 0 : if (Vtype==VID::JPSI)
1212 : {
1213 : //leptonic case
1214 0 : pdf = (1-3*A)*cos(theta)*cos(theta) + (1+A);
1215 0 : }
1216 : else
1217 : {
1218 : //hadronic case
1219 0 : pdf = (3*A-1)*cos(theta)*cos(theta) + (1-A);
1220 :
1221 : }
1222 0 : report(Severity::Debug,fname.c_str())<<" V decay pdf value : "<<pdf<<std::endl;
1223 :
1224 : //set probability
1225 0 : setProb(pdf);
1226 : return;
1227 0 : }
|
