Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtSVSCPiso.cpp
1 //--------------------------------------------------------------------------
2 //
3 // Environment:
4 //      This software is part of the EvtGen package developed jointly
5 //      for the BaBar and CLEO collaborations.  If you use all or part
6 //      of it, please give an appropriate acknowledgement.
7 //
8 // Copyright Information: See EvtGen/COPYRIGHT
9 //      Copyright (C) 1998      Caltech, UCSB
10 //
11 // Module: EvtSVSCPiso.cc
12 //
13 // Description: Routine to decay scalar -> vectors scalar
14 //              with CP violation and isospin amplitudes.
15 //              More specifically, it is indended to handle
16 //              decays like B->rho pi and B->a1 pi.
17 //
18 // Modification history:
19 //
20 //    RYD/NK    Febuary 16, 1998          Module created
21 //
22 //------------------------------------------------------------------------
23 // 
24 #include "EvtGenBase/EvtPatches.hh"
25 #include <stdlib.h>
26 #include "EvtGenBase/EvtParticle.hh"
27 #include "EvtGenBase/EvtRandom.hh"
28 #include "EvtGenBase/EvtGenKine.hh"
29 #include "EvtGenBase/EvtCPUtil.hh"
30 #include "EvtGenBase/EvtPDL.hh"
31 #include "EvtGenBase/EvtReport.hh"
32 #include "EvtGenBase/EvtVector4C.hh"
33 #include "EvtGenModels/EvtSVSCPiso.hh"
34 #include "EvtGenBase/EvtId.hh"
35 #include <string>
36 #include "EvtGenBase/EvtConst.hh"
37
38 EvtSVSCPiso::~EvtSVSCPiso() {}
39
40 std::string EvtSVSCPiso::getName(){
41
42   return "SVS_CP_ISO";     
43
44 }
45
46
47 EvtDecayBase* EvtSVSCPiso::clone(){
48
49   return new EvtSVSCPiso;
50
51 }
52
53 void EvtSVSCPiso::init(){
54
55   // check that there are 27 arguments
56   checkNArg(27);
57   checkNDaug(2);
58
59   checkSpinParent(EvtSpinType::SCALAR);
60
61   checkSpinDaughter(0,EvtSpinType::VECTOR);
62   checkSpinDaughter(1,EvtSpinType::SCALAR);
63
64 }
65
66
67 void EvtSVSCPiso::initProbMax(){
68
69 //this might need some revision..
70
71 if ((EvtPDL::chg3(getDaug(0)) > 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
72   setProbMax(2.0*(getArg(3)*getArg(3) + 4.0*getArg(23)*getArg(23)));
73 }
74
75 if ((EvtPDL::chg3(getDaug(0)) < 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
76   setProbMax(2.0*(getArg(5)*getArg(5) + 4.0*getArg(25)*getArg(25)));
77 }
78
79 if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) > 0)) {
80   setProbMax(2.0*(getArg(7)*getArg(7) + 4.0*getArg(23)*getArg(23)));
81 }
82
83 if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) < 0)) {
84   setProbMax(2.0*(getArg(9)*getArg(9) + 4.0*getArg(25)*getArg(25)));
85 }
86
87 if ((EvtPDL::chg3(getDaug(0)) > 0) && (EvtPDL::chg3(getDaug(1)) < 0)) {
88   setProbMax(2.0*(getArg(11)*getArg(11) + getArg(23)*getArg(23) + 
89                   getArg(19)*getArg(19) + getArg(13)*getArg(13) + 
90                   getArg(25)*getArg(25) + getArg(21)*getArg(21)));
91 }
92
93 if ((EvtPDL::chg3(getDaug(0)) < 0) && (EvtPDL::chg3(getDaug(1)) > 0)) {
94   setProbMax(2.0*(getArg(15)*getArg(15) + getArg(23)*getArg(23) + 
95                   getArg(19)*getArg(19) + getArg(17)*getArg(17) + 
96                   getArg(25)*getArg(25) + getArg(21)*getArg(21)));
97 }
98
99 if ((EvtPDL::chg3(getDaug(0)) == 0) && (EvtPDL::chg3(getDaug(1)) == 0)) {
100    setProbMax(2.0*(getArg(7)*getArg(7) + getArg(3)*getArg(3) + getArg(11)*getArg(11) + 
101                   getArg(15)*getArg(15) + 4.0*getArg(19)*getArg(19) + getArg(9)*getArg(9)+
102                    getArg(5)*getArg(5) + getArg(13)*getArg(13) + getArg(17)*getArg(17) + 
103                    4.0*getArg(21)*getArg(21)));
104 }
105
106 }
107
108
109 void EvtSVSCPiso::decay( EvtParticle *p){
110
111   //added by Lange Jan4,2000
112   static EvtId B0=EvtPDL::getId("B0");
113   static EvtId B0B=EvtPDL::getId("anti-B0");
114
115   double t;
116   EvtId other_b;
117   int charged(0);
118
119   int first_time=0;
120   int flip=0;
121   EvtId ds[2];
122
123
124 //randomly generate the tag (B0 or B0B) 
125
126    double tag = EvtRandom::Flat(0.0,1.0);
127    if (tag < 0.5) {
128  
129      EvtCPUtil::getInstance()->OtherB(p,t,other_b,1.0);
130      other_b = B0;
131    }
132    else {
133     
134      EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.0);
135      other_b = B0B;
136    }
137
138   if (p->getNDaug()==0) first_time=1;
139
140   if (first_time){
141     if (EvtRandom::Flat(0.0,1.0)<getArg(2)) flip=1;
142   }
143   else{
144     if (getDaug(0)!=p->getDaug(0)->getId()) flip=1;
145   }
146
147   if (!flip) {
148     ds[0]=getDaug(0);
149     ds[1]=getDaug(1);
150   }
151   else{
152     ds[0]=EvtPDL::chargeConj(getDaug(0));
153     ds[1]=EvtPDL::chargeConj(getDaug(1));
154   }
155
156   p->initializePhaseSpace(getNDaug(),ds);
157
158   EvtParticle *v,*s;
159   v=p->getDaug(0);
160   s=p->getDaug(1);
161
162    EvtComplex amp;
163
164    EvtComplex A_f,Abar_f;
165    EvtComplex A_fbar,Abar_fbar;
166    EvtComplex Apm, Apm_bar, Amp, Amp_bar;
167
168    EvtComplex Tp0, Tp0_bar, T0p, T0p_bar,Tpm, Tpm_bar, Tmp, Tmp_bar;
169    EvtComplex P1, P1_bar, P0, P0_bar;
170
171    Tp0 = EvtComplex(getArg(3)*cos(getArg(4)),getArg(3)*sin(getArg(4)));
172    Tp0_bar = EvtComplex(getArg(5)*cos(getArg(6)),getArg(5)*sin(getArg(6)));
173    T0p = EvtComplex(getArg(7)*cos(getArg(8)),getArg(7)*sin(getArg(8)));
174    T0p_bar = EvtComplex(getArg(9)*cos(getArg(10)),getArg(9)*sin(getArg(10)));
175    Tpm = EvtComplex(getArg(11)*cos(getArg(12)),getArg(11)*sin(getArg(12)));
176    Tpm_bar = EvtComplex(getArg(13)*cos(getArg(14)),getArg(13)*sin(getArg(14)));
177    Tmp = EvtComplex(getArg(15)*cos(getArg(16)),getArg(15)*sin(getArg(16)));
178    Tmp_bar = EvtComplex(getArg(17)*cos(getArg(18)),getArg(17)*sin(getArg(18)));
179    P0 = EvtComplex(getArg(19)*cos(getArg(20)),getArg(19)*sin(getArg(20)));
180    P0_bar = EvtComplex(getArg(21)*cos(getArg(22)),getArg(21)*sin(getArg(22)));
181    P1 = EvtComplex(getArg(23)*cos(getArg(24)),getArg(23)*sin(getArg(24)));
182    P1_bar = EvtComplex(getArg(25)*cos(getArg(26)),getArg(25)*sin(getArg(26)));
183
184
185 //***********************charged modes****************************
186
187  if ((EvtPDL::chg3(getDaug(0)) > 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
188
189 //V+ S0, so T+0 + 2 P1
190    
191    charged = 1;
192    A_f = Tp0 + 2.0*P1;
193  }
194
195  if ((EvtPDL::chg3(getDaug(0)) < 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
196
197 //V- S0, so T+0_bar + 2P1_bar
198    
199    charged = 1;
200    A_f = Tp0_bar + 2.0*P1_bar; 
201  }
202  
203  if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) > 0)) {
204
205 //V0 S+, so T0+ - 2 P1
206
207    charged = 1;
208    A_f = T0p - 2.0*P1;
209  }
210
211  if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) < 0)) {
212
213 //V0 S-, so T0+_bar - 2 P1_bar
214
215    charged = 1;
216    A_f = T0p_bar - 2.0*P1_bar;
217  }
218
219
220 //***********************neutral modes***************************
221
222
223 //V+ S-, so Af = T+- + P1 + P0
224 Apm = Tpm + P1 + P0;
225 Apm_bar = Tpm_bar + P1_bar + P0_bar;
226
227 //V- S+, so Af = T-+ - P1 + P0
228 Amp = Tmp - P1 + P0;
229 Amp_bar = Tmp_bar - P1_bar + P0;
230
231
232  if ((EvtPDL::chg3(getDaug(0)) > 0 ) && (EvtPDL::chg3(getDaug(1)) < 0)) {
233
234 //V+ S-
235    charged = 0;
236    A_f = Apm; 
237    Abar_f = Apm_bar;
238    A_fbar = Amp;
239    Abar_fbar = Amp_bar;
240
241  }
242
243  if ((EvtPDL::chg3(getDaug(0)) < 0 ) && (EvtPDL::chg3(getDaug(1)) > 0)) {
244
245 //V- S+
246    charged = 0;
247    A_f = Amp; 
248    Abar_f = Amp_bar;
249    A_fbar = Apm;
250    Abar_fbar = Apm_bar;
251
252  }
253
254  if ((EvtPDL::chg3(getDaug(0)) == 0 ) && (EvtPDL::chg3(getDaug(1)) == 0)) {
255
256 //V0 S0
257    charged = 0;
258    A_f = T0p + Tp0 - Tpm - Tmp - 2.0*P0 ; 
259    Abar_f = T0p_bar + Tp0_bar - Tpm_bar - Tmp_bar - 2.0*P0_bar;
260    A_fbar = A_f;
261    Abar_fbar = Abar_f;
262
263  }
264
265 if (charged==0) {
266    
267    if (!flip) {
268      if (other_b==B0B){
269
270        amp=A_f*cos(getArg(1)*t/(2*EvtConst::c))+
271          EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
272          EvtComplex(0.0,1.0)*Abar_f*sin(getArg(1)*t/(2*EvtConst::c));
273      }
274      if (other_b==B0){
275             
276        amp=A_f*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
277          EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+       
278          Abar_f*cos(getArg(1)*t/(2*EvtConst::c));
279      }
280    }
281    else{
282      if (other_b==B0B){
283
284        amp=A_fbar*cos(getArg(1)*t/(2*EvtConst::c))+
285          EvtComplex(cos(-2.0*getArg(0)),sin(-2.0*getArg(0)))*
286          EvtComplex(0.0,1.0)*Abar_fbar*sin(getArg(1)*t/(2*EvtConst::c));
287      }
288      if (other_b==B0){
289
290        amp=A_fbar*EvtComplex(cos(2.0*getArg(0)),sin(2.0*getArg(0)))*
291          EvtComplex(0.0,1.0)*sin(getArg(1)*t/(2*EvtConst::c))+       
292          Abar_fbar*cos(getArg(1)*t/(2*EvtConst::c));
293      }
294    }
295
296 }
297 else amp = A_f;
298
299   EvtVector4R p4_parent;
300
301   p4_parent=v->getP4()+s->getP4();
302
303   double norm=1.0/v->getP4().d3mag();
304
305    vertex(0,amp*norm*p4_parent*(v->epsParent(0)));
306    vertex(1,amp*norm*p4_parent*(v->epsParent(1)));
307    vertex(2,amp*norm*p4_parent*(v->epsParent(2)));
308
309   return ;
310 }
311