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