Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtLambdaB2LambdaV.cpp
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 EvtLambdaB2LambdaV::EvtLambdaB2LambdaV()
22 {
23   //set facility name
24   fname="EvtGen.EvtLambdaB2LambdaV";  
25 }
26
27
28 //------------------------------------------------------------------------
29 // Destructor
30 //------------------------------------------------------------------------
31 EvtLambdaB2LambdaV::~EvtLambdaB2LambdaV()
32 {}
33
34
35 //------------------------------------------------------------------------
36 // Method 'getName'
37 //------------------------------------------------------------------------
38 std::string EvtLambdaB2LambdaV::getName()
39 {
40   return  "LAMBDAB2LAMBDAV";
41 }
42
43
44 //------------------------------------------------------------------------
45 // Method 'clone'
46 //------------------------------------------------------------------------
47 EvtDecayBase* EvtLambdaB2LambdaV::clone()
48 {
49   return new EvtLambdaB2LambdaV;
50
51 }
52
53
54 //------------------------------------------------------------------------
55 // Method 'initProbMax'
56 //------------------------------------------------------------------------
57 void 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 //------------------------------------------------------------------------
69 void 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 //------------------------------------------------------------------------
190 void 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 //------------------------------------------------------------------------
328 double 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 //------------------------------------------------------------------------
339 double 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 //------------------------------------------------------------------------
367 double 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 //------------------------------------------------------------------------
437 EvtLambda2PPiForLambdaB2LambdaV::EvtLambda2PPiForLambdaB2LambdaV()
438
439   //set facility name
440   fname="EvtGen.EvtLambda2PPiForLambdaB2LambdaV";  
441 }
442
443
444 //------------------------------------------------------------------------
445 // Destructor
446 //------------------------------------------------------------------------
447 EvtLambda2PPiForLambdaB2LambdaV::~EvtLambda2PPiForLambdaB2LambdaV()
448 {
449 }
450
451
452 //------------------------------------------------------------------------
453 // Method 'getName'
454 //------------------------------------------------------------------------
455 std::string EvtLambda2PPiForLambdaB2LambdaV::getName()
456 {
457   return "LAMBDA2PPIFORLAMBDAB2LAMBDAV";
458 }
459
460
461 //------------------------------------------------------------------------
462 // Method 'clone'
463 //------------------------------------------------------------------------
464 EvtDecayBase* EvtLambda2PPiForLambdaB2LambdaV::clone()
465 {
466   return new EvtLambda2PPiForLambdaB2LambdaV;
467 }
468
469
470 //------------------------------------------------------------------------
471 // Method 'initProbMax'
472 //------------------------------------------------------------------------
473 void 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 //------------------------------------------------------------------------
496 void 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 //------------------------------------------------------------------------
617 void 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
738 report(INFO,fname.c_str())<<" Lambdab  px: "<<px_M<<std::endl;
739 report(INFO,fname.c_str())<<" Lambdab  py: "<<py_M<<std::endl;
740 report(INFO,fname.c_str())<<" Lambdab  pz: "<<pz_M<<std::endl;
741 report(INFO,fname.c_str())<<" Lambdab  E: "<<E_M<<std::endl;
742 report(INFO,fname.c_str())<<" Lambdab2  px:  "<<q_lambdab2.get(1)<<std::endl;
743 report(INFO,fname.c_str())<<" Lambdab2  py:  "<<q_lambdab2.get(2)<<std::endl;
744 report(INFO,fname.c_str())<<" Lambdab2  pz:  "<<q_lambdab2.get(3)<<std::endl;
745 report(INFO,fname.c_str())<<" Lambdab2  E:   "<<q_lambdab2.get(0)<<std::endl;
746 report(INFO,fname.c_str())<<" Lambdab3  px:  "<<q_lambdab3.get(1)<<std::endl;
747 report(INFO,fname.c_str())<<" Lambdab3  py:  "<<q_lambdab3.get(2)<<std::endl;
748 report(INFO,fname.c_str())<<" Lambdab3  pz:  "<<q_lambdab3.get(3)<<std::endl;
749 report(INFO,fname.c_str())<<" Lambdab3  E:   "<<q_lambdab3.get(0)<<std::endl;
750 report(INFO,fname.c_str())<<" Lambda 0  px:  "<<P_lambda.get(1)<<std::endl;
751 report(INFO,fname.c_str())<<" Lambda 0  py:  "<<P_lambda.get(2)<<std::endl;
752 report(INFO,fname.c_str())<<" Lambda 0  pz:  "<<P_lambda.get(3)<<std::endl;
753 report(INFO,fname.c_str())<<" Lambda 0  E:   "<<P_lambda.get(0)<<std::endl;
754 report(INFO,fname.c_str())<<" Lambda 1  px:  "<<q_lambda1.get(1)<<std::endl;
755 report(INFO,fname.c_str())<<" Lambda 1  py:  "<<q_lambda1.get(2)<<std::endl;
756 report(INFO,fname.c_str())<<" Lambda 1  pz:  "<<q_lambda1.get(3)<<std::endl;
757 report(INFO,fname.c_str())<<" Lambda 1  E:   "<<q_lambda1.get(0)<<std::endl;
758 report(INFO,fname.c_str())<<" Lambda 2  px:  "<<q_lambda2.get(1)<<std::endl;
759 report(INFO,fname.c_str())<<" Lambda 2  py:  "<<q_lambda2.get(2)<<std::endl;
760 report(INFO,fname.c_str())<<" Lambda 2  pz:  "<<q_lambda2.get(3)<<std::endl;
761 report(INFO,fname.c_str())<<" Lambda 2  E:   "<<q_lambda2.get(0)<<std::endl;
762
763 report(INFO,fname.c_str())<<" Lambda  px: "<<px<<std::endl;
764 report(INFO,fname.c_str())<<" Lambda  py: "<<py<<std::endl;
765 report(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  
772 report(INFO,fname.c_str())<<" pion 3 px:  "<<q_pion3.get(1)<<std::endl;
773 report(INFO,fname.c_str())<<" pion 3 px:  "<<q_pion3.get(1)<<std::endl;
774 report(INFO,fname.c_str())<<" pion 3 py:  "<<q_pion3.get(2)<<std::endl;
775 report(INFO,fname.c_str())<<" pion 3 pz:  "<<q_pion3.get(3)<<std::endl;
776 report(INFO,fname.c_str())<<" pion 3 E:   "<<q_pion3.get(0)<<std::endl;
777
778 report(INFO,fname.c_str())<<" pion 5 px:  "<<q_pion5.get(1)<<std::endl;
779 report(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
790 report(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  
795 report(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
801 report(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 //------------------------------------------------------------------------
848 EvtV2VpVmForLambdaB2LambdaV::EvtV2VpVmForLambdaB2LambdaV()
849 {
850   //set facility name
851   fname="EvtGen.EvtV2VpVmForLambdaB2LambdaV";  
852 }
853
854
855 //------------------------------------------------------------------------
856 // Destructor
857 //------------------------------------------------------------------------
858 EvtV2VpVmForLambdaB2LambdaV::~EvtV2VpVmForLambdaB2LambdaV()
859 {}
860
861
862 //------------------------------------------------------------------------
863 // Method 'getName'
864 //------------------------------------------------------------------------
865 std::string EvtV2VpVmForLambdaB2LambdaV::getName()
866 {
867   return "V2VPVMFORLAMBDAB2LAMBDAV";
868 }
869
870 //------------------------------------------------------------------------
871 // Method 'clone'
872 //------------------------------------------------------------------------
873 EvtDecayBase* EvtV2VpVmForLambdaB2LambdaV::clone()
874 {
875   return new EvtV2VpVmForLambdaB2LambdaV;
876 }
877
878
879 //------------------------------------------------------------------------
880 // Method 'initProbMax'
881 //------------------------------------------------------------------------
882 void 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 //------------------------------------------------------------------------
905 void 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 //------------------------------------------------------------------------
1017 void 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;
1147 report(INFO,fname.c_str())<<" Lambdab  py: "<<py_M<<std::endl;
1148 report(INFO,fname.c_str())<<" Lambdab  pz: "<<pz_M<<std::endl;
1149 report(INFO,fname.c_str())<<" Lambdab  E: "<<E_M<<std::endl;
1150 report(INFO,fname.c_str())<<" Lambdab2  px:  "<<q_lambdab2.get(1)<<std::endl;
1151 report(INFO,fname.c_str())<<" Lambdab2  py:  "<<q_lambdab2.get(2)<<std::endl;
1152 report(INFO,fname.c_str())<<" Lambdab2  pz:  "<<q_lambdab2.get(3)<<std::endl;
1153 report(INFO,fname.c_str())<<" Lambdab2  E:   "<<q_lambdab2.get(0)<<std::endl;
1154 report(INFO,fname.c_str())<<" Lambdab3  px:  "<<q_lambdab3.get(1)<<std::endl;
1155 report(INFO,fname.c_str())<<" Lambdab3  py:  "<<q_lambdab3.get(2)<<std::endl;
1156 report(INFO,fname.c_str())<<" Lambdab3  pz:  "<<q_lambdab3.get(3)<<std::endl;
1157 report(INFO,fname.c_str())<<" Lambdab3  E:   "<<q_lambdab3.get(0)<<std::endl;
1158 report(INFO,fname.c_str())<<" V 0  px:  "<<P_V.get(1)<<std::endl;
1159 report(INFO,fname.c_str())<<" V 0  py:  "<<P_V.get(2)<<std::endl;
1160 report(INFO,fname.c_str())<<" V 0  pz:  "<<P_V.get(3)<<std::endl;
1161 report(INFO,fname.c_str())<<" V 0  E:   "<<P_V.get(0)<<std::endl;
1162 report(INFO,fname.c_str())<<" V 1  px:  "<<q_V1.get(1)<<std::endl;
1163 report(INFO,fname.c_str())<<" V 1  py:  "<<q_V1.get(2)<<std::endl;
1164 report(INFO,fname.c_str())<<" V 1  pz:  "<<q_V1.get(3)<<std::endl;
1165 report(INFO,fname.c_str())<<" V 1  E:   "<<q_V1.get(0)<<std::endl;
1166 report(INFO,fname.c_str())<<" V 2  px:  "<<q_V2.get(1)<<std::endl;
1167 report(INFO,fname.c_str())<<" V 2  py:  "<<q_V2.get(2)<<std::endl;
1168 report(INFO,fname.c_str())<<" V 2  pz:  "<<q_V2.get(3)<<std::endl;
1169 report(INFO,fname.c_str())<<" V 2  E:   "<<q_V2.get(0)<<std::endl;
1170 report(INFO,fname.c_str())<<" V  px: "<<px<<std::endl;
1171 report(INFO,fname.c_str())<<" V  py: "<<py<<std::endl;
1172 report(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; 
1177 report(INFO,fname.c_str())<<" Vm 3 px:  "<<q_Vm3.get(1)<<std::endl;
1178 report(INFO,fname.c_str())<<" Vm 3 px:  "<<q_Vm3.get(1)<<std::endl;
1179 report(INFO,fname.c_str())<<" Vm 3 py:  "<<q_Vm3.get(2)<<std::endl;
1180 report(INFO,fname.c_str())<<" Vm 3 pz:  "<<q_Vm3.get(3)<<std::endl;
1181 report(INFO,fname.c_str())<<" Vm 3 E:   "<<q_Vm3.get(0)<<std::endl;
1182 report(INFO,fname.c_str())<<" Vm 5 px:  "<<q_Vm5.get(1)<<std::endl;
1183 report(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;
1190 report(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 }