Converting TEvtGen to native cmake
[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtGoityRoberts.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: EvtGoityRoberts.cc
12 //
13 // Description: Routine to decay vector-> scalar scalar
14 //
15 // Modification history:
16 //
17 //    RYD     November 24, 1996        Module created
18 //
19 //------------------------------------------------------------------------
20 // 
21 #include "EvtGenBase/EvtPatches.hh"
22 #include <stdlib.h>
23 #include "EvtGenBase/EvtParticle.hh"
24 #include "EvtGenBase/EvtGenKine.hh"
25 #include "EvtGenBase/EvtPDL.hh"
26 #include "EvtGenBase/EvtReport.hh"
27 #include "EvtGenModels/EvtGoityRoberts.hh"
28 #include "EvtGenBase/EvtTensor4C.hh"
29 #include "EvtGenBase/EvtDiracSpinor.hh"
30 #include <string>
31 #include "EvtGenBase/EvtVector4C.hh"
32
33 EvtGoityRoberts::~EvtGoityRoberts() {}
34
35 std::string EvtGoityRoberts::getName(){
36
37   return "GOITY_ROBERTS";     
38
39 }
40
41
42 EvtDecayBase* EvtGoityRoberts::clone(){
43
44   return new EvtGoityRoberts;
45
46 }
47
48 void EvtGoityRoberts::init(){
49
50   // check that there are 0 arguments
51   checkNArg(0);
52   checkNDaug(4);
53
54   checkSpinParent(EvtSpinType::SCALAR);
55   checkSpinDaughter(1,EvtSpinType::SCALAR);
56   checkSpinDaughter(2,EvtSpinType::DIRAC);
57   checkSpinDaughter(3,EvtSpinType::NEUTRINO);
58
59 }
60
61
62 void EvtGoityRoberts::initProbMax() {
63
64    setProbMax( 3000.0);
65 }      
66
67 void EvtGoityRoberts::decay( EvtParticle *p){
68
69   //added by Lange Jan4,2000
70   static EvtId DST0=EvtPDL::getId("D*0");
71   static EvtId DSTB=EvtPDL::getId("anti-D*0");
72   static EvtId DSTP=EvtPDL::getId("D*+");
73   static EvtId DSTM=EvtPDL::getId("D*-");
74   static EvtId D0=EvtPDL::getId("D0");
75   static EvtId D0B=EvtPDL::getId("anti-D0");
76   static EvtId DP=EvtPDL::getId("D+");
77   static EvtId DM=EvtPDL::getId("D-");
78
79
80
81   EvtId meson=getDaug(0);
82
83   if (meson==DST0||meson==DSTP||meson==DSTM||meson==DSTB) {
84     DecayBDstarpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
85   }
86   else{
87     if (meson==D0||meson==DP||meson==DM||meson==D0B) {
88       DecayBDpilnuGR(p,getDaug(0),getDaug(2),getDaug(3));
89     }
90     else{
91       report(ERROR,"EvtGen") << "Wrong daugther in EvtGoityRoberts!\n";
92     }
93   }
94   return ;
95 }
96
97 void EvtGoityRoberts::DecayBDstarpilnuGR(EvtParticle *pb,EvtId ndstar,
98                                          EvtId nlep, EvtId /*nnu*/)
99 {
100
101   pb->initializePhaseSpace(getNDaug(),getDaugs());
102
103   //added by Lange Jan4,2000
104   static EvtId EM=EvtPDL::getId("e-");
105   static EvtId EP=EvtPDL::getId("e+");
106   static EvtId MUM=EvtPDL::getId("mu-");
107   static EvtId MUP=EvtPDL::getId("mu+");
108
109   EvtParticle *dstar, *pion, *lepton, *neutrino;
110   
111   // pb->makeDaughters(getNDaug(),getDaugs());
112   dstar=pb->getDaug(0);
113   pion=pb->getDaug(1);
114   lepton=pb->getDaug(2);
115   neutrino=pb->getDaug(3);
116
117   EvtVector4C l1, l2, et0, et1, et2;
118   
119   EvtVector4R v,vp,p4_pi;
120   double w;
121   
122   v.set(1.0,0.0,0.0,0.0);       //4-velocity of B meson
123   vp=(1.0/dstar->getP4().mass())*dstar->getP4();  //4-velocity of D*
124   p4_pi=pion->getP4();
125
126   w=v*vp;                       //four velocity transfere.
127
128   EvtTensor4C omega;
129
130   double mb=EvtPDL::getMeanMass(pb->getId());     //B mass
131   double md=EvtPDL::getMeanMass(ndstar);   //D* mass
132
133   EvtComplex dmb(0.0460,-0.5*0.00001);   // B*-B mass splitting ?
134   EvtComplex dmd(0.1421,-0.5*0.00006);
135                              // The last two sets of numbers should
136                              // be correctly calculated from the
137                              // dstar and pion charges.
138   double g = 0.5;         // EvtAmplitude proportional to these coupling constants
139   double alpha3 =  0.690; // See table I in G&R's paper
140   double alpha1 = -1.430;
141   double alpha2 = -0.140;
142   double f0 = 0.093;      // The pion decay constants set to 93 MeV
143
144   EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2  
145   EvtComplex dmt1(0.392,-0.5*1.040);
146   EvtComplex dmt2(0.709,-0.5*0.405);
147                    
148   double betas=0.285;      // magic number for meson wave function ground state
149   double betap=0.280;      // magic number for meson wave function state "1"
150   double betad=0.260;      // magic number for meson wave function state "2"
151   double betasp=betas*betas+betap*betap;
152   double betasd=betas*betas+betad*betad;
153
154   double lambdabar=0.750;  //M(0-,1-) - mQ From Goity&Roberts's code
155
156 // Isgur&Wise fct
157   double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
158   double xi1= -1.0*sqrt(2.0/3.0)*(
159        lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
160        exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
161   double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
162                pow((2*betas*betap/(betasp)),2.5)*
163                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
164   double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
165                pow((2*betas*betad/(betasd)),3.5)*
166                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
167
168   //report(INFO,"EvtGen") <<"rho's:"<<rho1<<rho2<<endl;
169
170   EvtComplex h1nr,h2nr,h3nr,f1nr,f2nr;
171   EvtComplex f3nr,f4nr,f5nr,f6nr,knr,g1nr,g2nr,g3nr,g4nr,g5nr;
172   EvtComplex h1r,h2r,h3r,f1r,f2r,f3r,f4r,f5r,f6r,kr,g1r,g2r,g3r,g4r,g5r;
173   EvtComplex h1,h2,h3,f1,f2,f3,f4,f5,f6,k,g1,g2,g3,g4,g5;
174
175   // Non-resonance part
176   h1nr = -g*xi*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmb));
177   h2nr = -g*xi/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmb));
178   h3nr = -(g*xi/(f0*md))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
179                          -EvtComplex((1.0+w)/(p4_pi*vp),0.0));
180
181   f1nr = -(g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb) -
182          1.0/(EvtComplex(p4_pi*vp,0.0)+dmd));
183   f2nr = f1nr*mb/md;
184   f3nr = EvtComplex(0.0);
185   f4nr = EvtComplex(0.0);
186   f5nr = (g*xi/(2*f0*mb*md))*(EvtComplex(1.0,0.0)
187                  +(p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmb));
188   f6nr = (g*xi/(2*f0*mb))*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb)
189                  -EvtComplex(1.0/(p4_pi*vp),0.0));
190
191   knr = (g*xi/(2*f0))*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmb) +
192                  EvtComplex((p4_pi*(v-w*vp))/(p4_pi*vp),0.0));
193   
194   g1nr = EvtComplex(0.0);
195   g2nr = EvtComplex(0.0);
196   g3nr = EvtComplex(0.0);
197   g4nr = (g*xi)/(f0*md*EvtComplex(p4_pi*vp));
198   g5nr = EvtComplex(0.0);
199
200   // Resonance part (D** removed by hand - alainb)
201   h1r = -alpha1*rho1*(p4_pi*v)/(f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
202         alpha2*rho2*(p4_pi*(v+2.0*w*v-vp))
203         /(3*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
204         alpha3*xi1*(p4_pi*v)/(f0*mb*md*EvtComplex(p4_pi*v,0.0)+dmt3); 
205   h2r = -alpha2*(1+w)*rho2/(3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
206         alpha3*xi1/(f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
207   h3r = alpha2*rho2*(1+w)/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
208         alpha3*xi1/(f0*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
209
210   f1r = -alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) -
211         alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
212   f2r = f1r*mb/md;
213   f3r = EvtComplex(0.0);
214   f4r = EvtComplex(0.0);
215   f5r = alpha1*rho1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt1)) +
216         alpha2*rho2*(p4_pi*(vp-v/3.0-2.0/3.0*w*v))/
217         (2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
218         alpha3*xi1*(p4_pi*v)/(2*f0*mb*md*(EvtComplex(p4_pi*v,0.0)+dmt3));
219   f6r = alpha2*rho2*(w-1.0)/(6*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
220         alpha3*xi1/(2*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt3));
221
222   kr = -alpha1*rho1*(w-1.0)*(p4_pi*v)/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt1)) -
223        alpha2*rho2*(w-1.0)*(p4_pi*(vp-w*v))
224        /(3*f0*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
225        alpha3*xi1*(p4_pi*(vp-w*v))/(2*f0*(EvtComplex(p4_pi*v,0.0)+dmt3));
226   
227   g1r = EvtComplex(0.0);
228   g2r = EvtComplex(0.0);
229   g3r = -g2r;
230   g4r = 2.0*alpha2*rho2/(3*f0*md*(EvtComplex(p4_pi*v,0.0)+dmt2));
231   g5r = EvtComplex(0.0);
232
233   //Sum
234   h1 = h1nr + h1r;
235   h2 = h2nr + h2r;
236   h3 = h3nr + h3r;
237
238   f1 = f1nr + f1r;
239   f2 = f2nr + f2r;
240   f3 = f3nr + f3r;
241   f4 = f4nr + f4r;
242   f5 = f5nr + f5r;
243   f6 = f6nr + f6r;
244
245   k = knr+kr;
246   
247   g1 = g1nr + g1r;
248   g2 = g2nr + g2r;
249   g3 = g3nr + g3r;
250   g4 = g4nr + g4r;
251   g5 = g5nr + g5r;
252
253   EvtTensor4C g_metric;
254   g_metric.setdiag(1.0,-1.0,-1.0,-1.0);
255
256   if (nlep==EM||nlep==MUM){ 
257     omega=EvtComplex(0.0,0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
258                              h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
259                              h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
260         f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
261                        f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
262         f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
263         EvtComplex(0.0,0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
264                               (g1*p4_pi+g2*mb*v))+
265         EvtComplex(0.0,0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
266                              dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
267
268    l1=EvtLeptonVACurrent(lepton->spParent(0),neutrino->spParentNeutrino());
269    l2=EvtLeptonVACurrent(lepton->spParent(1),neutrino->spParentNeutrino());
270   }
271   else{
272     if (nlep==EP||nlep==MUP){ 
273       omega=EvtComplex(0.0,-0.5)*dual(h1*mb*md*EvtGenFunctions::directProd(v,vp)+
274                              h2*mb*EvtGenFunctions::directProd(v,p4_pi)+
275                                       h3*md*EvtGenFunctions::directProd(vp,p4_pi))+
276         f1*mb*EvtGenFunctions::directProd(v,p4_pi)+f2*md*EvtGenFunctions::directProd(vp,p4_pi)+
277                        f3*EvtGenFunctions::directProd(p4_pi,p4_pi)+f4*mb*mb*EvtGenFunctions::directProd(v,v)+
278         f5*mb*md*EvtGenFunctions::directProd(vp,v)+f6*mb*EvtGenFunctions::directProd(p4_pi,v)+k*g_metric+
279         EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd(dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v),
280                               (g1*p4_pi+g2*mb*v))+
281         EvtComplex(0.0,-0.5)*EvtGenFunctions::directProd((g3*mb*v+g4*md*vp+g5*p4_pi),
282                              dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v));
283
284    l1=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(0));
285    l2=EvtLeptonVACurrent(neutrino->spParentNeutrino(),lepton->spParent(1));
286     }
287     else{
288    report(DEBUG,"EvtGen") << "42387dfs878w wrong lepton number\n";
289     }
290  }
291
292   et0=omega.cont2( dstar->epsParent(0).conj() );
293   et1=omega.cont2( dstar->epsParent(1).conj() );
294   et2=omega.cont2( dstar->epsParent(2).conj() );
295
296   vertex(0,0,l1.cont(et0));
297   vertex(0,1,l2.cont(et0));
298
299   vertex(1,0,l1.cont(et1));
300   vertex(1,1,l2.cont(et1));
301
302   vertex(2,0,l1.cont(et2));
303   vertex(2,1,l2.cont(et2));
304
305   return;
306
307 }
308
309 void EvtGoityRoberts::DecayBDpilnuGR(EvtParticle *pb,EvtId nd,
310                                      EvtId nlep, EvtId /*nnu*/)
311
312 {
313   //added by Lange Jan4,2000
314   static EvtId EM=EvtPDL::getId("e-");
315   static EvtId EP=EvtPDL::getId("e+");
316   static EvtId MUM=EvtPDL::getId("mu-");
317   static EvtId MUP=EvtPDL::getId("mu+");
318
319   EvtParticle *d, *pion, *lepton, *neutrino;
320
321   pb->initializePhaseSpace(getNDaug(),getDaugs());
322   d=pb->getDaug(0);
323   pion=pb->getDaug(1);
324   lepton=pb->getDaug(2);
325   neutrino=pb->getDaug(3);
326
327   EvtVector4C l1, l2, et0, et1, et2;
328  
329   EvtVector4R v,vp,p4_pi;
330   double w;
331   
332   v.set(1.0,0.0,0.0,0.0);       //4-velocity of B meson
333   vp=(1.0/d->getP4().mass())*d->getP4();  //4-velocity of D
334   p4_pi=pion->getP4();                  //4-momentum of pion
335   w=v*vp;                       //four velocity transfer.
336   
337   double mb=EvtPDL::getMeanMass(pb->getId());     //B mass
338   double md=EvtPDL::getMeanMass(nd);   //D* mass
339   EvtComplex dmb(0.0460,-0.5*0.00001);   //B mass splitting ?
340                       //The last two numbers should be
341                       //correctly calculated from the
342                       //dstar and pion particle number.
343
344   double g = 0.5;         // Amplitude proportional to these coupling constants
345   double alpha3 =  0.690; // See table I in G&R's paper
346   double alpha1 = -1.430;
347   double alpha2 = -0.140;
348   double f0=0.093;        // The pion decay constant set to 93 MeV
349
350   EvtComplex dmt3(0.563,-0.5*0.191); // Mass splitting = dmt - iGamma/2  
351   EvtComplex dmt1(0.392,-0.5*1.040);
352   EvtComplex dmt2(0.709,-0.5*0.405);
353                    
354   double betas=0.285;      // magic number for meson wave function ground state
355   double betap=0.280;      // magic number for meson wave function state "1"
356   double betad=0.260;      // magic number for meson wave function state "2"
357   double betasp=betas*betas+betap*betap;
358   double betasd=betas*betas+betad*betad;
359
360   double lambdabar=0.750;  //M(0-,1-) - mQ From Goity&Roberts's code
361
362   // Isgur&Wise fct
363   double xi = exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
364   double xi1= -1.0*sqrt(2.0/3.0)*(lambdabar*lambdabar*(w*w-1.0)/(4*betas*betas))*
365               exp(lambdabar*lambdabar*(1.0-w*w)/(4*betas*betas));
366   double rho1= sqrt(1.0/2.0)*(lambdabar/betas)*
367                pow((2*betas*betap/(betasp)),2.5)*
368                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasp));
369   double rho2= sqrt(1.0/8.0)*(lambdabar*lambdabar/(betas*betas))*
370                pow((2*betas*betad/(betasd)),3.5)*
371                exp(lambdabar*lambdabar*(1.0-w*w)/(2*betasd));
372
373   EvtComplex h,a1,a2,a3;
374   EvtComplex hnr,a1nr,a2nr,a3nr;
375   EvtComplex hr,a1r,a2r,a3r;
376
377 // Non-resonance part (D* and D** removed by hand - alainb)
378   hnr = g*xi*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb*md);
379   a1nr= -1.0*g*xi*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0);
380   a2nr= g*xi*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmb))/(2*f0*mb);
381   a3nr=EvtComplex(0.0,0.0);
382
383 // Resonance part (D** remove by hand - alainb)
384   hr = alpha2*rho2*(w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0*mb*md) +
385        alpha3*xi1*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb*md);
386   a1r= -1.0*alpha2*rho2*(w*w-1)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt2))/(6*f0) -
387        alpha3*xi1*(1+w)*(1.0/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0);
388   a2r= alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*mb) +
389        alpha2*rho2*(0.5*p4_pi*(w*vp-v)+p4_pi*(vp-w*v))/
390                   (3*f0*mb*(EvtComplex(p4_pi*v,0.0)+dmt2)) +
391        alpha3*xi1*((p4_pi*(v+vp))/(EvtComplex(p4_pi*v,0.0)+dmt3))/(2*f0*mb);
392   a3r= -1.0*alpha1*rho1*((p4_pi*v)/(EvtComplex(p4_pi*v,0.0)+dmt1))/(2*f0*md) -
393        alpha2*rho2*((p4_pi*(vp-w*v))/(EvtComplex(p4_pi*v,0.0)+dmt2))/(2*f0*md);
394
395 // Sum
396   h=hnr+hr;
397   a1=a1nr+a1r;
398   a2=a2nr+a2r;
399   a3=a3nr+a3r;
400
401   EvtVector4C omega;
402
403   if ( nlep==EM|| nlep==MUM ) {
404     omega=EvtComplex(0.0,-1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
405                  a1*p4_pi+a2*mb*v+a3*md*vp;
406     l1=EvtLeptonVACurrent(
407              lepton->spParent(0),neutrino->spParentNeutrino());
408     l2=EvtLeptonVACurrent(
409              lepton->spParent(1),neutrino->spParentNeutrino());
410   }
411   else{
412     if ( nlep==EP|| nlep==MUP ) {
413      omega=EvtComplex(0.0,1.0)*h*mb*md*dual(EvtGenFunctions::directProd(vp,p4_pi)).cont2(v)+
414                  a1*p4_pi+a2*mb*v+a3*md*vp;
415      l1=EvtLeptonVACurrent(
416               neutrino->spParentNeutrino(),lepton->spParent(0));
417      l2=EvtLeptonVACurrent(
418               neutrino->spParentNeutrino(),lepton->spParent(1));
419     }
420     else{
421      report(ERROR,"EvtGen") << "42387dfs878w wrong lepton number\n";
422     }
423   }
424
425   vertex(0,l1*omega);
426   vertex(1,l2*omega);
427
428 return;
429
430 }
431