Fix for definitions for CINT
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtSSD_DirectCP.cpp
1 // $Id: EvtSSD_DirectCP.cpp,v 1.2 2009-03-16 16:24:05 robbep Exp $
2 // Generation of direct CP violation in hadronic environment
3 // Patrick Robbe, LHCb,  08 Nov 2006
4 // 
5 #include <stdlib.h>
6 #include "EvtGenBase/EvtParticle.hh"
7 #include "EvtGenBase/EvtRandom.hh"
8 #include "EvtGenBase/EvtGenKine.hh"
9 #include "EvtGenBase/EvtCPUtil.hh"
10 #include "EvtGenBase/EvtPDL.hh"
11 #include "EvtGenBase/EvtReport.hh"
12 #include "EvtGenBase/EvtVector4C.hh"
13 #include "EvtGenBase/EvtTensor4C.hh"
14 #include "EvtGenModels/EvtSSD_DirectCP.hh"
15 #include <string>
16 #include "EvtGenBase/EvtConst.hh"
17
18 EvtSSD_DirectCP::~EvtSSD_DirectCP() {}
19
20 std::string EvtSSD_DirectCP::getName( ){
21
22   return "SSD_DirectCP" ;
23
24 }
25
26
27 EvtDecayBase* EvtSSD_DirectCP::clone(){
28
29   return new EvtSSD_DirectCP;
30
31 }
32
33 void EvtSSD_DirectCP::init(){
34
35   // check that there is 1 argument and 2-body decay
36
37   checkNArg(1);
38   checkNDaug(2);
39
40   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
41   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
42   
43   if ( (!(d1type == EvtSpinType::SCALAR || d2type == EvtSpinType::SCALAR))||
44        (!((d2type==EvtSpinType::SCALAR)||(d2type==EvtSpinType::VECTOR)||
45           (d2type==EvtSpinType::TENSOR)))||
46        (!((d1type==EvtSpinType::SCALAR)||(d1type==EvtSpinType::VECTOR)||
47           (d1type==EvtSpinType::TENSOR)))
48        ) {
49     report(ERROR,"EvtGen") << "EvtSSD_DirectCP generator expected "
50                            << "one of the daugters to be a scalar, "
51                            << "the other either scalar, vector, or tensor, "
52                            << "found:"
53                            << EvtPDL::name(getDaug(0)).c_str()
54                            <<" and "
55                            <<EvtPDL::name(getDaug(1)).c_str()<<std::endl;
56     report(ERROR,"EvtGen") << "Will terminate execution!"<<std::endl;
57     ::abort();
58   }
59   
60   _acp = getArg( 0 ) ; // A_CP defined as A_CP = (BR(fbar)-BR(f))/(BR(fbar)+BR(f))
61
62 }
63
64 void EvtSSD_DirectCP::initProbMax() {
65   double theProbMax = 1. ;
66
67   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
68   EvtSpinType::spintype d1type=EvtPDL::getSpinType(getDaug(0));
69   if (d1type==EvtSpinType::TENSOR||d2type==EvtSpinType::TENSOR) theProbMax*=10;
70
71   setProbMax(theProbMax);
72 }
73
74 void EvtSSD_DirectCP::decay( EvtParticle *p) {
75
76   bool flip = false ;
77   EvtId daugs[2];
78   
79   // decide it is B or Bbar:
80   if ( EvtRandom::Flat(0.,1.) < ( ( 1. - _acp ) / 2. ) ) {
81     // it is a B
82     if ( EvtPDL::getStdHep( getParentId() ) < 0 ) flip = true ;
83   } else {
84     // it is a Bbar
85     if ( EvtPDL::getStdHep( getParentId() ) > 0 ) flip = true ;
86   }
87   
88   if ( flip ) {
89     if ( ( isB0Mixed( p ) ) || ( isBsMixed( p ) ) ) {
90       p->getParent()
91         ->setId( EvtPDL::chargeConj( p->getParent()->getId() ) ) ;
92       p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
93     }
94     else {
95       p->setId( EvtPDL::chargeConj( p->getId() ) ) ;
96     }
97   }  
98
99   if (!flip) {
100     daugs[0]=getDaug(0);
101     daugs[1]=getDaug(1);
102   }
103   else{
104     daugs[0]=EvtPDL::chargeConj(getDaug(0));
105     daugs[1]=EvtPDL::chargeConj(getDaug(1));
106   }
107
108   EvtParticle *d;
109   p->initializePhaseSpace(2, daugs);
110
111   EvtVector4R p4_parent=p->getP4Restframe();
112   double m_parent=p4_parent.mass();
113
114   EvtSpinType::spintype d2type=EvtPDL::getSpinType(getDaug(1));
115
116   EvtVector4R momv;
117   EvtVector4R moms;
118
119   if (d2type==EvtSpinType::SCALAR){
120     d2type=EvtPDL::getSpinType(getDaug(0));
121     d= p->getDaug(0);
122     momv = d->getP4();
123     moms = p->getDaug(1)->getP4();
124   }
125   else{
126     d= p->getDaug(1);
127     momv = d->getP4();
128     moms = p->getDaug(0)->getP4();
129   }
130
131   if (d2type==EvtSpinType::SCALAR) {
132     vertex(1.);
133   }
134   
135   if (d2type==EvtSpinType::VECTOR) {
136     
137     double norm=momv.mass()/(momv.d3mag()*p->mass());
138     
139     vertex(0,norm*p4_parent*(d->epsParent(0)));
140     vertex(1,norm*p4_parent*(d->epsParent(1)));
141     vertex(2,norm*p4_parent*(d->epsParent(2)));
142   
143   }
144
145   if (d2type==EvtSpinType::TENSOR) {
146
147     double norm=
148       d->mass()*d->mass()/(m_parent*d->getP4().d3mag()*d->getP4().d3mag());
149  
150    
151    vertex(0,norm*d->epsTensorParent(0).cont1(p4_parent)*p4_parent);
152    vertex(1,norm*d->epsTensorParent(1).cont1(p4_parent)*p4_parent);
153    vertex(2,norm*d->epsTensorParent(2).cont1(p4_parent)*p4_parent);
154    vertex(3,norm*d->epsTensorParent(3).cont1(p4_parent)*p4_parent);
155    vertex(4,norm*d->epsTensorParent(4).cont1(p4_parent)*p4_parent);  
156   }
157 }
158
159 bool EvtSSD_DirectCP::isB0Mixed ( EvtParticle * p ) {
160   if ( ! ( p->getParent() ) ) return false ;
161
162   static EvtId B0 =EvtPDL::getId("B0");
163   static EvtId B0B=EvtPDL::getId("anti-B0");
164
165   if ( ( p->getId() != B0 ) && ( p->getId() != B0B ) ) return false ;
166
167   if ( ( p->getParent()->getId() == B0 ) ||
168        ( p->getParent()->getId() == B0B ) ) return true ;
169
170   return false ;
171 }
172
173 bool EvtSSD_DirectCP::isBsMixed ( EvtParticle * p ) {
174   if ( ! ( p->getParent() ) ) return false ;
175
176   static EvtId BS0=EvtPDL::getId("B_s0");
177   static EvtId BSB=EvtPDL::getId("anti-B_s0");
178
179   if ( ( p->getId() != BS0 ) && ( p->getId() != BSB ) ) return false ;
180
181   if ( ( p->getParent()->getId() == BS0 ) ||
182        ( p->getParent()->getId() == BSB ) ) return true ;
183
184   return false ;
185 }
186
187 std::string EvtSSD_DirectCP::getParamName(int i) {
188   switch(i) {
189   case 0:
190     return "ACP";
191   default:
192     return "";
193   }
194 }