Fix bug in building local list of valid files.
[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtSVVHelCPMix.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: EvtSVVHelCPMix.cc
12 //
13 // Description: Routine to decay scalar -> 2 vectors
14 //              by specifying the helicity amplitudes, taking appropriate
15 //              weak phases into account to get mixing and CP violation through
16 //              interference. Based on EvtSVVHelAmp. Particularly appropriate for
17 //              Bs->J/Psi+Phi
18 //
19 // Modification history:
20 //
21 //    RYD       November 24, 1996       EvtSVVHelAmp Module
22 //    CATMORE   March 2004              Modified to EvtSVVHelCPMix
23 //
24 // Model takes the following as user-specified arguments:
25 //      deltaM, averageM - mass diference and average of light and heavy mass eigenstates (real scalars)
26 //      gamma, deltagamma - average width and width difference of the l and h eigenstates (real scalars)
27 //      delta1, delta2 - strong phases (real scalars)
28 //      direct weak phase (real scalar) (for Bs->JPsiPhi this will be zero)     
29 //      weak mixing phase (real scalar) (this is equal to 2*arg(Vts Vtb) for Bs->JPsiPhi)
30 //      Magnitudes of helicity amplitudes as in SVV_HELAMP
31 // See Phys Rev D 34 p1404 - p1417 and chapters 5 and 7 of Physics Reports 370 p537-680 for more details
32 //------------------------------------------------------------------------
33 // 
34 #ifdef WIN32 
35   #pragma warning( disable : 4786 ) 
36   // Disable anoying warning about symbol size 
37 #endif 
38 #include <iostream>
39 #include <fstream>
40 #include <stdlib.h>
41 #include <ctype.h>
42 #include "EvtGenBase/EvtParticle.hh"
43 #include "EvtGenBase/EvtGenKine.hh"
44 #include "EvtGenBase/EvtPDL.hh"
45 #include "EvtGenBase/EvtVector4C.hh"
46 #include "EvtGenBase/EvtTensor4C.hh"
47 #include "EvtGenBase/EvtVector3C.hh"
48 #include "EvtGenBase/EvtVector3R.hh"
49 #include "EvtGenBase/EvtTensor3C.hh"
50 #include "EvtGenBase/EvtReport.hh"
51 #include "EvtGenModels/EvtSVVHelCPMix.hh"
52 #include "EvtGenBase/EvtId.hh"
53 #include <string>
54
55 EvtSVVHelCPMix::~EvtSVVHelCPMix() {}
56
57 std::string EvtSVVHelCPMix::getName(){
58
59   return "SVVHELCPMIX";     
60
61 }
62
63
64 EvtDecayBase* EvtSVVHelCPMix::clone(){
65
66   return new EvtSVVHelCPMix;
67
68 }
69
70 void EvtSVVHelCPMix::init(){
71
72   // check that there are 12 arguments
73   checkNArg(12);
74   checkNDaug(2);
75
76   checkSpinParent(EvtSpinType::SCALAR);
77
78   checkSpinDaughter(0,EvtSpinType::VECTOR);
79   checkSpinDaughter(1,EvtSpinType::VECTOR);
80
81   hp = EvtComplex(getArg(0)*cos(getArg(1)),getArg(0)*sin(getArg(1)));
82   h0 = EvtComplex(getArg(2)*cos(getArg(3)),getArg(2)*sin(getArg(3)));
83   hm = EvtComplex(getArg(4)*cos(getArg(5)),getArg(4)*sin(getArg(5)));
84   averageM = getArg(6);
85   deltaM = getArg(7);
86   gamma = getArg(8);
87   deltagamma = getArg(9);
88   weakmixingphase = EvtComplex(cos(getArg(10)),sin(getArg(10)));
89   weakdirectphase = EvtComplex(cos(getArg(11)),sin(getArg(11)));
90 }
91
92
93 void EvtSVVHelCPMix::initProbMax(){
94
95   setProbMax(getArg(0)*getArg(0)+getArg(2)*getArg(2)+getArg(4)*getArg(4));
96   
97 }
98
99
100 void EvtSVVHelCPMix::decay( EvtParticle *p){
101
102   EvtParticle* parent = p;
103   EvtAmp& amp = _amp2;
104   EvtId n_v1 = getDaug(0);
105   EvtId n_v2 = getDaug(1);
106         
107   //  Routine to decay a vector into a vector and scalar.  Started
108   //  by ryd on Oct 17, 1996.
109   // Modified by J.Catmore to take account of CP-violation and mixing
110     
111   int tndaug = 2;
112   EvtId tdaug[2];
113   EvtId Bs=EvtPDL::getId("B_s0");
114   EvtId antiBs=EvtPDL::getId("anti-B_s0");
115   tdaug[0] = n_v1;
116   tdaug[1] = n_v2;
117
118 // Phase space and kinematics
119
120   parent->initializePhaseSpace(tndaug,tdaug);
121
122   EvtParticle *v1,*v2;
123   v1 = parent->getDaug(0);
124   v2 = parent->getDaug(1);
125
126   EvtVector4R momv1 = v1->getP4();
127
128   EvtVector3R v1dir(momv1.get(1),momv1.get(2),momv1.get(3));
129   v1dir=v1dir/v1dir.d3mag();
130
131   // Definition of quantities used in construction of complex amplitudes:
132
133   EvtTensor3C M;  // Tensor as defined in EvtGen manual, equ 117
134   EvtComplex a,b,c; // Helicity amplitudes; EvtGen manual eqns 126-128, also see Phys Lett B 369 p144-150 eqn 15
135   //EvtComplex deltamu = EvtComplex(deltaM, -0.5*deltagamma); // See Phys Rev D 34 p1404
136
137   // conversion from times in mm/c to natural units [GeV]^-1
138   double t = ((parent->getLifetime())/2.998e11)*6.58e-25; 
139
140   // The following two quantities defined in Phys Rev D 34 p1404
141   EvtComplex fplus =  EvtComplex(cos(averageM*t),-1.*sin(averageM*t))*exp(-(gamma/2.0)*t)*
142     (cos(0.5*deltaM*t)*cosh(0.25*deltagamma*t)+EvtComplex(0.0,sin(0.5*deltaM*t)*sinh(0.25*deltagamma*t)));
143   EvtComplex fminus = EvtComplex(cos(averageM*t), -1.*sin(averageM*t))*exp(-(gamma/2.0)*t)*EvtComplex(0.0,1.0)*
144     (sin(0.5*deltaM*t)*cosh(0.25*deltagamma*t)-EvtComplex(0.0,1.0)*sinh(0.25*deltagamma*t)*cos(0.5*deltaM*t)); 
145
146 // See EvtGen manual pp 106-107
147
148   a=-0.5*(hp+hm);
149   b=EvtComplex(0.0,0.5)*(hp-hm);
150   c=(h0+0.5*(hp+hm));
151   
152   M=a*EvtTensor3C::id()+
153     b*EvtGenFunctions::eps(v1dir)+
154     c*EvtGenFunctions::directProd(v1dir,v1dir);
155                 
156   EvtVector3C t0=M.cont1(v1->eps(0).vec().conj());
157   EvtVector3C t1=M.cont1(v1->eps(1).vec().conj());
158   EvtVector3C t2=M.cont1(v1->eps(2).vec().conj());
159
160   EvtVector3C eps0=v2->eps(0).vec().conj();
161   EvtVector3C eps1=v2->eps(1).vec().conj();
162   EvtVector3C eps2=v2->eps(2).vec().conj();
163
164 // We need two sets of equations, one for mesons which were in the Bs state at t=0, and another
165 // for those which were in the antiBs state. Each equation consists of a sum of amplitudes - mod-squaring gives the interference terms.
166
167         EvtComplex amplSum00, amplSum01, amplSum02;
168         EvtComplex amplSum10, amplSum11, amplSum12;
169         EvtComplex amplSum20, amplSum21, amplSum22;
170
171   // First the Bs state:
172   
173   if (parent->getId()==Bs) {
174         
175         amplSum00 = (fplus*weakdirectphase*t0*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps0);
176         amplSum01 = (fplus*weakdirectphase*t0*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps1);
177         amplSum02 = (fplus*weakdirectphase*t0*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t0*eps2);
178
179         amplSum10 = (fplus*weakdirectphase*t1*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps0);
180         amplSum11 = (fplus*weakdirectphase*t1*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps1);
181         amplSum12 = (fplus*weakdirectphase*t1*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t1*eps2);
182
183         amplSum20 = (fplus*weakdirectphase*t2*eps0) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps0);
184         amplSum21 = (fplus*weakdirectphase*t2*eps1) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps1);
185         amplSum22 = (fplus*weakdirectphase*t2*eps2) + (fminus*(1.0/weakdirectphase)*weakmixingphase*t2*eps2);
186   }
187
188   // Now the anti-Bs state:
189
190   if (parent->getId()==antiBs) {
191   
192         amplSum00 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps0) + (fplus*(1.0/weakdirectphase)*t0*eps0);
193         amplSum01 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps1) + (fplus*(1.0/weakdirectphase)*t0*eps1);
194         amplSum02 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t0*eps2) + (fplus*(1.0/weakdirectphase)*t0*eps2);
195
196         amplSum10 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps0) + (fplus*(1.0/weakdirectphase)*t1*eps0);
197         amplSum11 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps1) + (fplus*(1.0/weakdirectphase)*t1*eps1);
198         amplSum12 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t1*eps2) + (fplus*(1.0/weakdirectphase)*t1*eps2);
199
200         amplSum20 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps0) + (fplus*(1.0/weakdirectphase)*t2*eps0);
201         amplSum21 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps1) + (fplus*(1.0/weakdirectphase)*t2*eps1);
202         amplSum22 = (fminus*weakdirectphase*(1.0/weakmixingphase)*t2*eps2) + (fplus*(1.0/weakdirectphase)*t2*eps2);
203
204   }
205
206   // Now set the amplitude
207
208   amp.vertex(0,0,amplSum00);
209   report(INFO,"EvtGen") << "00: " << amplSum00 << std::endl;  
210   amp.vertex(0,1,amplSum01);
211   report(INFO,"EvtGen") << "01: " << amplSum01 << std::endl;
212   amp.vertex(0,2,amplSum02);
213   report(INFO,"EvtGen") << "02: " << amplSum02 << std::endl;
214
215   amp.vertex(1,0,amplSum10);
216   report(INFO,"EvtGen") << "10: " << amplSum10 << std::endl;
217   amp.vertex(1,1,amplSum11);
218   report(INFO,"EvtGen") << "11: " << amplSum11 << std::endl;
219   amp.vertex(1,2,amplSum12);
220   report(INFO,"EvtGen") << "12: " << amplSum12 << std::endl;
221
222   amp.vertex(2,0,amplSum20);
223   report(INFO,"EvtGen") << "20: " << amplSum20 << std::endl;
224   amp.vertex(2,1,amplSum21);
225   report(INFO,"EvtGen") << "21: " << amplSum21 << std::endl;
226   amp.vertex(2,2,amplSum22);
227   report(INFO,"EvtGen") << "22: " << amplSum22 << std::endl;
228
229   return ;
230
231 }
232
233 std::string EvtSVVHelCPMix::getParamName(int i) {
234   switch(i) {
235   case 0:
236     return "plusHelAmp";
237   case 1:
238     return "plusHelAmpPhase";
239   case 2:
240     return "zeroHelAmp";
241   case 3:
242     return "zeroHelAmpPhase";
243   case 4:
244     return "minusHelAmp";
245   case 5:
246     return "minusHelAmpPhase";
247   case 6:
248     return "averageM";
249   case 7:
250     return "deltaM";
251   case 8:
252     return "gamma";
253   case 9:
254     return "deltaGamma";
255   case 10:
256     return "weakMixPhase";
257   case 11:
258     return "weakDirectPhase";
259   default:
260     return "";
261   }
262 }
263
264 std::string EvtSVVHelCPMix::getParamDefault(int i) {
265   switch(i) {
266   case 0:
267     return "1.0";
268   case 1:
269     return "0.0";
270   case 2:
271     return "1.0";
272   case 3:
273     return "0.0";
274   case 4:
275     return "1.0";
276   case 5:
277     return "0.0";
278   default:
279     return "";
280   }
281 }