[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenBase / EvtSemiLeptonicTensorAmp.cpp

//--------------------------------------------------------------------------
//
// Environment:
//      This software is part of the EvtGen package developed jointly
//      for the BaBar and CLEO collaborations.  If you use all or part
//      of it, please give an appropriate acknowledgement.
//
// Copyright Information: See EvtGen/COPYRIGHT
//      Copyright (C) 1998      Caltech, UCSB
//
// Module: EvtSemiLeptonicTensorAmp.cc
//
// Description: Routine to implement semileptonic decays to pseudo-scalar
//              mesons. 
//
// Modification history:
//
//    DJL       April 17,1998       Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtSemiLeptonicTensorAmp.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtAmp.hh"
#include "EvtGenBase/EvtSemiLeptonicFF.hh"

void EvtSemiLeptonicTensorAmp::CalcAmp( EvtParticle *parent,
					EvtAmp& amp,
					EvtSemiLeptonicFF *FormFactors ) {
  static EvtId EM=EvtPDL::getId("e-");
  static EvtId MUM=EvtPDL::getId("mu-");
  static EvtId TAUM=EvtPDL::getId("tau-");
  static EvtId EP=EvtPDL::getId("e+");
  static EvtId MUP=EvtPDL::getId("mu+");
  static EvtId TAUP=EvtPDL::getId("tau+");

  static EvtId D0=EvtPDL::getId("D0");
  static EvtId D0B=EvtPDL::getId("anti-D0");
  static EvtId DP=EvtPDL::getId("D+");
  static EvtId DM=EvtPDL::getId("D-");
  static EvtId DSM=EvtPDL::getId("D_s-");
  static EvtId DSP=EvtPDL::getId("D_s+");

  //Add the lepton and neutrino 4 momenta to find q2

  EvtVector4R q = parent->getDaug(1)->getP4() 
                    + parent->getDaug(2)->getP4(); 
  double q2 = (q.mass2());

  double hf,kf,bpf,bmf;

  FormFactors->gettensorff(parent->getId(),
                           parent->getDaug(0)->getId(),
                           q2,
			   parent->getDaug(0)->mass(),
                           &hf, 
                           &kf, 
                           &bpf, 
                           &bmf);


  double costhl_flag = 1.0;

  if(parent->getId()==D0||parent->getId()==D0B||
     parent->getId()==DP||parent->getId()==DM) {
    costhl_flag = -1.0;
  }
  if(parent->getId()==DSP||parent->getId()==DSM) {
    costhl_flag = -1.0;
  }
  hf = hf * costhl_flag;

  EvtVector4R p4b;
  p4b.set(parent->mass(),0.0,0.0,0.0);
 
  EvtVector4R p4meson = parent->getDaug(0)->getP4();
 
  EvtVector4C l1,l2;

  EvtId l_num = parent->getDaug(1)->getId();

  EvtVector4C ep_meson_b[5];

  ep_meson_b[0] = ((parent->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
  ep_meson_b[1] = ((parent->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
  ep_meson_b[2] = ((parent->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
  ep_meson_b[3] = ((parent->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
  ep_meson_b[4] = ((parent->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();

  EvtVector4R pp,pm;

  pp=p4b+p4meson;
  pm=p4b-p4meson;

  //lange - October 31,2002 - try to lessen the mass dependence of probmax
  double q2max = p4b.mass2() + p4meson.mass2() - 2.0*p4b.mass()*p4meson.mass();
  double q2maxin=1.0/q2max;

  EvtComplex ep_meson_bb[5];

  ep_meson_bb[0]=ep_meson_b[0]*(p4b);
  ep_meson_bb[1]=ep_meson_b[1]*(p4b);
  ep_meson_bb[2]=ep_meson_b[2]*(p4b);
  ep_meson_bb[3]=ep_meson_b[3]*(p4b);
  ep_meson_bb[4]=ep_meson_b[4]*(p4b);


  EvtVector4C tds0,tds1,tds2,tds3,tds4;

  EvtTensor4C tds;
  if (l_num==EM||l_num==MUM||l_num==TAUM){
    EvtTensor4C tdual=EvtComplex(0.0,hf)*dual(EvtGenFunctions::directProd(pp,pm));
    tds0=tdual.cont2(ep_meson_b[0])
      -kf*ep_meson_b[0]
      -bpf*ep_meson_bb[0]*pp-bmf*ep_meson_bb[0]*pm;
    tds0*=q2maxin;

    tds1=tdual.cont2(ep_meson_b[1])
      -kf*ep_meson_b[1]
      -bpf*ep_meson_bb[1]*pp-bmf*ep_meson_bb[1]*pm;
    tds1*=q2maxin;

    tds2=tdual.cont2(ep_meson_b[2])
      -kf*ep_meson_b[2]
      -bpf*ep_meson_bb[2]*pp-bmf*ep_meson_bb[2]*pm;
    tds2*=q2maxin;

    tds3=tdual.cont2(ep_meson_b[3])
      -kf*ep_meson_b[3]
      -bpf*ep_meson_bb[3]*pp-bmf*ep_meson_bb[3]*pm;
    tds3*=q2maxin;

    tds4=tdual.cont2(ep_meson_b[4])
      -kf*ep_meson_b[4]
      -bpf*ep_meson_bb[4]*pp-bmf*ep_meson_bb[4]*pm;
    tds4*=q2maxin;


    l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
                          parent->getDaug(2)->spParentNeutrino());
    l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
                          parent->getDaug(2)->spParentNeutrino());
  }
  else{
    if (l_num==EP||l_num==MUP||l_num==TAUP){
     EvtTensor4C tdual=EvtComplex(0.0,-hf)*dual(EvtGenFunctions::directProd(pp,pm));
      tds0=tdual.cont2(ep_meson_b[0])
        -kf*ep_meson_b[0]
        -bpf*ep_meson_bb[0]*pp-bmf*ep_meson_bb[0]*pm;
      tds0*=q2maxin;

      tds1=tdual.cont2(ep_meson_b[1])
        -kf*ep_meson_b[1]
        -bpf*ep_meson_bb[1]*pp-bmf*ep_meson_bb[1]*pm;
      tds1*=q2maxin;

      tds2=tdual.cont2(ep_meson_b[2])
        -kf*ep_meson_b[2]
        -bpf*ep_meson_bb[2]*pp-bmf*ep_meson_bb[2]*pm;
      tds2*=q2maxin;

      tds3=tdual.cont2(ep_meson_b[3])
        -kf*ep_meson_b[3]
        -bpf*ep_meson_bb[3]*pp-bmf*ep_meson_bb[3]*pm;
      tds3*=q2maxin;

      tds4=tdual.cont2(ep_meson_b[4])
        -kf*ep_meson_b[4]
        -bpf*ep_meson_bb[4]*pp-bmf*ep_meson_bb[4]*pm;
      tds4*=q2maxin;

      l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
			    parent->getDaug(1)->spParent(0));
      l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
			    parent->getDaug(1)->spParent(1));
    }
    else{
      report(ERROR,"EvtGen") << "dfnb89agngri wrong lepton number\n";
    }
  }
 
  amp.vertex(0,0,l1*tds0);
  amp.vertex(0,1,l2*tds0);

  amp.vertex(1,0,l1*tds1);
  amp.vertex(1,1,l2*tds1);

  amp.vertex(2,0,l1*tds2);
  amp.vertex(2,1,l2*tds2);

  amp.vertex(3,0,l1*tds3);
  amp.vertex(3,1,l2*tds3);

  amp.vertex(4,0,l1*tds4);
  amp.vertex(4,1,l2*tds4);

  return;
 
}
Commit	Line	Data
da0e9ce3	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: EvtSemiLeptonicTensorAmp.cc
	12	//
	13	// Description: Routine to implement semileptonic decays to pseudo-scalar
	14	// mesons.
	15	//
	16	// Modification history:
	17	//
	18	// DJL April 17,1998 Module created
	19	//
	20	//------------------------------------------------------------------------
	21	//
	22	#include "EvtGenBase/EvtPatches.hh"
	23	#include "EvtGenBase/EvtPatches.hh"
	24	#include "EvtGenBase/EvtParticle.hh"
	25	#include "EvtGenBase/EvtGenKine.hh"
	26	#include "EvtGenBase/EvtPDL.hh"
	27	#include "EvtGenBase/EvtReport.hh"
	28	#include "EvtGenBase/EvtTensor4C.hh"
	29	#include "EvtGenBase/EvtVector4C.hh"
	30	#include "EvtGenBase/EvtDiracSpinor.hh"
	31	#include "EvtGenBase/EvtSemiLeptonicTensorAmp.hh"
	32	#include "EvtGenBase/EvtId.hh"
	33	#include "EvtGenBase/EvtAmp.hh"
	34	#include "EvtGenBase/EvtSemiLeptonicFF.hh"
	35
	36	void EvtSemiLeptonicTensorAmp::CalcAmp( EvtParticle *parent,
	37	EvtAmp& amp,
	38	EvtSemiLeptonicFF *FormFactors ) {
	39	static EvtId EM=EvtPDL::getId("e-");
	40	static EvtId MUM=EvtPDL::getId("mu-");
	41	static EvtId TAUM=EvtPDL::getId("tau-");
	42	static EvtId EP=EvtPDL::getId("e+");
	43	static EvtId MUP=EvtPDL::getId("mu+");
	44	static EvtId TAUP=EvtPDL::getId("tau+");
	45
	46	static EvtId D0=EvtPDL::getId("D0");
	47	static EvtId D0B=EvtPDL::getId("anti-D0");
	48	static EvtId DP=EvtPDL::getId("D+");
	49	static EvtId DM=EvtPDL::getId("D-");
	50	static EvtId DSM=EvtPDL::getId("D_s-");
	51	static EvtId DSP=EvtPDL::getId("D_s+");
	52
	53	//Add the lepton and neutrino 4 momenta to find q2
	54
	55	EvtVector4R q = parent->getDaug(1)->getP4()
	56	+ parent->getDaug(2)->getP4();
	57	double q2 = (q.mass2());
	58
	59	double hf,kf,bpf,bmf;
	60
	61	FormFactors->gettensorff(parent->getId(),
	62	parent->getDaug(0)->getId(),
	63	q2,
	64	parent->getDaug(0)->mass(),
65	&hf,
66	&kf,
67	&bpf,
68	&bmf);
69
70
71	double costhl_flag = 1.0;
72
73	if(parent->getId()==D0\|\|parent->getId()==D0B\|\|
74	parent->getId()==DP\|\|parent->getId()==DM) {
75	costhl_flag = -1.0;
76	}
77	if(parent->getId()==DSP\|\|parent->getId()==DSM) {
78	costhl_flag = -1.0;
79	}
80	hf = hf * costhl_flag;
81
82	EvtVector4R p4b;
83	p4b.set(parent->mass(),0.0,0.0,0.0);
84
85	EvtVector4R p4meson = parent->getDaug(0)->getP4();
86
87	EvtVector4C l1,l2;
88
89	EvtId l_num = parent->getDaug(1)->getId();
90
91	EvtVector4C ep_meson_b[5];
92
93	ep_meson_b[0] = ((parent->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
94	ep_meson_b[1] = ((parent->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
95	ep_meson_b[2] = ((parent->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
96	ep_meson_b[3] = ((parent->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
97	ep_meson_b[4] = ((parent->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
98
99	EvtVector4R pp,pm;
100
101	pp=p4b+p4meson;
102	pm=p4b-p4meson;
103
104	//lange - October 31,2002 - try to lessen the mass dependence of probmax
105	double q2max = p4b.mass2() + p4meson.mass2() - 2.0p4b.mass()p4meson.mass();
106	double q2maxin=1.0/q2max;
107
108	EvtComplex ep_meson_bb[5];
109
110	ep_meson_bb[0]=ep_meson_b[0]*(p4b);
111	ep_meson_bb[1]=ep_meson_b[1]*(p4b);
112	ep_meson_bb[2]=ep_meson_b[2]*(p4b);
113	ep_meson_bb[3]=ep_meson_b[3]*(p4b);
114	ep_meson_bb[4]=ep_meson_b[4]*(p4b);
115
116
117	EvtVector4C tds0,tds1,tds2,tds3,tds4;
118
119	EvtTensor4C tds;
120	if (l_num==EM\|\|l_num==MUM\|\|l_num==TAUM){
0ca57c2f	121	EvtTensor4C tdual=EvtComplex(0.0,hf)*dual(EvtGenFunctions::directProd(pp,pm));
da0e9ce3	122	tds0=tdual.cont2(ep_meson_b[0])
	123	-kf*ep_meson_b[0]
	124	-bpfep_meson_bb[0]pp-bmfep_meson_bb[0]pm;
	125	tds0*=q2maxin;
	126
	127	tds1=tdual.cont2(ep_meson_b[1])
	128	-kf*ep_meson_b[1]
	129	-bpfep_meson_bb[1]pp-bmfep_meson_bb[1]pm;
	130	tds1*=q2maxin;
	131
	132	tds2=tdual.cont2(ep_meson_b[2])
	133	-kf*ep_meson_b[2]
	134	-bpfep_meson_bb[2]pp-bmfep_meson_bb[2]pm;
	135	tds2*=q2maxin;
	136
	137	tds3=tdual.cont2(ep_meson_b[3])
	138	-kf*ep_meson_b[3]
	139	-bpfep_meson_bb[3]pp-bmfep_meson_bb[3]pm;
	140	tds3*=q2maxin;
	141
	142	tds4=tdual.cont2(ep_meson_b[4])
	143	-kf*ep_meson_b[4]
	144	-bpfep_meson_bb[4]pp-bmfep_meson_bb[4]pm;
	145	tds4*=q2maxin;
	146
	147
	148	l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0),
	149	parent->getDaug(2)->spParentNeutrino());
	150	l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1),
	151	parent->getDaug(2)->spParentNeutrino());
	152	}
	153	else{
	154	if (l_num==EP\|\|l_num==MUP\|\|l_num==TAUP){
0ca57c2f	155	EvtTensor4C tdual=EvtComplex(0.0,-hf)*dual(EvtGenFunctions::directProd(pp,pm));
da0e9ce3	156	tds0=tdual.cont2(ep_meson_b[0])
	157	-kf*ep_meson_b[0]
	158	-bpfep_meson_bb[0]pp-bmfep_meson_bb[0]pm;
	159	tds0*=q2maxin;
	160
	161	tds1=tdual.cont2(ep_meson_b[1])
	162	-kf*ep_meson_b[1]
	163	-bpfep_meson_bb[1]pp-bmfep_meson_bb[1]pm;
	164	tds1*=q2maxin;
	165
	166	tds2=tdual.cont2(ep_meson_b[2])
	167	-kf*ep_meson_b[2]
	168	-bpfep_meson_bb[2]pp-bmfep_meson_bb[2]pm;
	169	tds2*=q2maxin;
	170
	171	tds3=tdual.cont2(ep_meson_b[3])
	172	-kf*ep_meson_b[3]
	173	-bpfep_meson_bb[3]pp-bmfep_meson_bb[3]pm;
	174	tds3*=q2maxin;
	175
	176	tds4=tdual.cont2(ep_meson_b[4])
	177	-kf*ep_meson_b[4]
	178	-bpfep_meson_bb[4]pp-bmfep_meson_bb[4]pm;
	179	tds4*=q2maxin;
	180
	181	l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
	182	parent->getDaug(1)->spParent(0));
	183	l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(),
	184	parent->getDaug(1)->spParent(1));
	185	}
	186	else{
	187	report(ERROR,"EvtGen") << "dfnb89agngri wrong lepton number\n";
	188	}
	189	}
	190
	191	amp.vertex(0,0,l1*tds0);
	192	amp.vertex(0,1,l2*tds0);
	193
	194	amp.vertex(1,0,l1*tds1);
	195	amp.vertex(1,1,l2*tds1);
	196
	197	amp.vertex(2,0,l1*tds2);
	198	amp.vertex(2,1,l2*tds2);
	199
	200	amp.vertex(3,0,l1*tds3);
	201	amp.vertex(3,1,l2*tds3);
	202
	203	amp.vertex(4,0,l1*tds4);
	204	amp.vertex(4,1,l2*tds4);
	205
	206	return;
	207
	208	}
	209