[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtBHadronic.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: EvtBHadronic.cc
//
// Description: Model for hadronic B decays. Uses naive factorization. 
//
// Modification history:
//
//    RYD     June 14, 1997         Module created
//
//------------------------------------------------------------------------
//
#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenModels/EvtISGW2FF.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenModels/EvtISGW2FF.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenModels/EvtBHadronic.hh"
#include "EvtGenBase/EvtReport.hh"
#include <string>
using std::endl;

EvtBHadronic::~EvtBHadronic() {}

std::string EvtBHadronic::getName(){

  return "BHADRONIC";    

}


EvtDecayBase* EvtBHadronic::clone(){

  return new EvtBHadronic;

}


void EvtBHadronic::init(){

  // check that there are 2 argument
  checkNArg(2);


}

void EvtBHadronic::decay( EvtParticle *p){

  //added by Lange Jan4,2000
  static EvtId B0=EvtPDL::getId("B0");
  static EvtId D0=EvtPDL::getId("D0");
  static EvtId DST0=EvtPDL::getId("D*0");
  static EvtId D3P10=EvtPDL::getId("D'_10");
  static EvtId D3P20=EvtPDL::getId("D_2*0");
  static EvtId D3P00=EvtPDL::getId("D_0*0");
  static EvtId D1P10=EvtPDL::getId("D_10");

  static EvtISGW2FF ffmodel;

  p->initializePhaseSpace(getNDaug(),getDaugs());

  EvtVector4R p4[MAX_DAUG];
  p->mass();


  int i,j;

  for ( i=0; i<getNDaug(); i++) {
    p4[i]=p->getDaug(i)->getP4();
  }

  int bcurrent,wcurrent;
  int nbcurrent=0;
  int nwcurrent=0;

  bcurrent=(int)getArg(0);
  wcurrent=(int)getArg(1);

  EvtVector4C jb[5],jw[5];
  EvtTensor4C g,tds;

  EvtVector4R p4b;
  p4b.set(p->mass(),0.0,0.0,0.0);

  EvtVector4R q;
  double q2;

  EvtComplex ep_meson_bb[5];
  double f,gf,ap,am;
  double fp,fm;
  double hf,kf,bp,bm;
  EvtVector4R pp,pm;
  EvtVector4C ep_meson_b[5];

  switch (bcurrent) {
  
  // D0
  case 1:
    q=p4b-p4[0];
    q2=q*q;
    nbcurrent=1;
    ffmodel.getscalarff(B0,D0,EvtPDL::getMeanMass(D0),
			EvtPDL::getMeanMass(getDaugs()[1]),&fp,&fm);
    jb[0]=EvtVector4C(fp*(p4b+p4[0])-fm*q);
    break;
    // D*
  case 2:
    q=p4b-p4[0];
    q2=q*q;
    nbcurrent=3;
    ffmodel.getvectorff(B0,DST0,EvtPDL::getMeanMass(DST0),q2,&f,&gf,&ap,&am);

    g.setdiag(1.0,-1.0,-1.0,-1.0);
    tds = -f*g 
      -ap*(EvtGenFunctions::directProd(p4b,p4b)+EvtGenFunctions::directProd(p4b,p4[0]))
      -gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
      -am*((EvtGenFunctions::directProd(p4b,p4b)-
            EvtGenFunctions::directProd(p4b,p4[0])));
    jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
    jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
    jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
    break;
    // D1
  case 3:
    q=p4b-p4[0];
    q2=q*q;
    nbcurrent=3;
    ffmodel.getvectorff(B0,D3P10,EvtPDL::getMeanMass(D3P10),q2,&f,&gf,&ap,&am);

    g.setdiag(1.0,-1.0,-1.0,-1.0);
    tds = -f*g 
      -ap*(EvtGenFunctions::directProd(p4b,p4b)+
           EvtGenFunctions::directProd(p4b,p4[0]))
      -gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
      -am*((EvtGenFunctions::directProd(p4b,p4b)-
            EvtGenFunctions::directProd(p4b,p4[0])));
    jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
    jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
    jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());

    break;
    // D2*
  case 4:
    q=p4b-p4[0];
    q2=q*q;
    nbcurrent=5;
    ffmodel.gettensorff(B0,D3P20,EvtPDL::getMeanMass(D3P20),q2,&hf,&kf,&bp,&bm);
    g.setdiag(1.0,-1.0,-1.0,-1.0);
    
    ep_meson_b[0] = ((p->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
    ep_meson_b[1] = ((p->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
    ep_meson_b[2] = ((p->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
    ep_meson_b[3] = ((p->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
    ep_meson_b[4] = ((p->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
  
    pp=p4b+p4[0];
    pm=p4b-p4[0];
    
    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);
    
    jb[0]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[0])
       -kf*ep_meson_b[0]
       -bp*ep_meson_bb[0]*pp-bm*ep_meson_bb[0]*pm;
     
    jb[1]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[1])
       -kf*ep_meson_b[1]
       -bp*ep_meson_bb[1]*pp-bm*ep_meson_bb[1]*pm;
    
    jb[2]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[2])
       -kf*ep_meson_b[2]
       -bp*ep_meson_bb[2]*pp-bm*ep_meson_bb[2]*pm;
    
    jb[3]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[3])
       -kf*ep_meson_b[3]
       -bp*ep_meson_bb[3]*pp-bm*ep_meson_bb[3]*pm;
    
    jb[4]=EvtComplex(0.0,(1.0)*hf)*dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[4])
       -kf*ep_meson_b[4]
       -bp*ep_meson_bb[4]*pp-bm*ep_meson_bb[4]*pm;
    break;
    // D_0*
  case 5:
    q=p4b-p4[0];
    q2=q*q;
    double f,gf,ap,am;
    nbcurrent=3;
    ffmodel.getvectorff(B0,D1P10,EvtPDL::getMeanMass(D1P10),q2,&f,&gf,&ap,&am);
    g.setdiag(1.0,-1.0,-1.0,-1.0);
    tds = -f*g 
      -ap*(EvtGenFunctions::directProd(p4b,p4b)+
           EvtGenFunctions::directProd(p4b,p4[0]))
      +gf*EvtComplex(0.0,1.0)*dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
      -am*((EvtGenFunctions::directProd(p4b,p4b)-
            EvtGenFunctions::directProd(p4b,p4[0])));
    jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
    jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
    jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
    break;
    // D_1
  case 6:
    q=p4b-p4[0];
    q2=q*q;
    nbcurrent=1;
    ffmodel.getscalarff(B0,D3P00,EvtPDL::getMeanMass(D3P00),q2,&fp,&fm);
    jb[0]=fp*(p4b+p4[0])+fm*q;
    break;
  default:
    report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown hadronic current."<<endl;
    
  }  

  double norm;

  switch (wcurrent) {


  case 1: case 3: case 4:
    nwcurrent=1;
    jw[0]=p4[getNDaug()-1];
    break;

  case 2: case 5: case 6:
    nwcurrent=3;
    norm=1.0/sqrt(p4[1].get(0)*p4[1].get(0)/p4[1].mass2()-1.0);
    jw[0]=norm*p->getDaug(getNDaug()-1)->epsParent(0);
    jw[1]=norm*p->getDaug(getNDaug()-1)->epsParent(1);
    jw[2]=norm*p->getDaug(getNDaug()-1)->epsParent(2);
    break;


  default:
    report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown W current."<<endl;

  }  

  if (nbcurrent==1&&nwcurrent==1){
    vertex(jb[0]*jw[0]);
    return;
  }

  if (nbcurrent==1){
    for(j=0;j<nwcurrent;j++){
      vertex(j,jb[0]*jw[j]);
    }
    return;
  }

  if (nwcurrent==1){
    for(j=0;j<nbcurrent;j++){
      vertex(j,jb[j]*jw[0]);
    }
    return;
  }

  for(j=0;j<nbcurrent;j++){
    for(i=0;i<nwcurrent;i++){
      vertex(j,i,jb[j]*jw[i]);
    }
  }
  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: EvtBHadronic.cc
	12	//
	13	// Description: Model for hadronic B decays. Uses naive factorization.
	14	//
	15	// Modification history:
	16	//
	17	// RYD June 14, 1997 Module created
	18	//
	19	//------------------------------------------------------------------------
	20	//
	21	#include "EvtGenBase/EvtPatches.hh"
	22	#include <stdlib.h>
	23	#include "EvtGenBase/EvtParticle.hh"
	24	#include "EvtGenModels/EvtISGW2FF.hh"
	25	#include "EvtGenBase/EvtGenKine.hh"
	26	#include "EvtGenBase/EvtPDL.hh"
	27	#include "EvtGenModels/EvtISGW2FF.hh"
	28	#include "EvtGenBase/EvtTensor4C.hh"
	29	#include "EvtGenBase/EvtVector4C.hh"
	30	#include "EvtGenBase/EvtVector4R.hh"
	31	#include "EvtGenModels/EvtBHadronic.hh"
	32	#include "EvtGenBase/EvtReport.hh"
	33	#include <string>
	34	using std::endl;
	35
	36	EvtBHadronic::~EvtBHadronic() {}
	37
	38	std::string EvtBHadronic::getName(){
	39
	40	return "BHADRONIC";
	41
	42	}
	43
	44
	45	EvtDecayBase* EvtBHadronic::clone(){
	46
	47	return new EvtBHadronic;
	48
	49	}
	50
	51
	52	void EvtBHadronic::init(){
	53
	54	// check that there are 2 argument
	55	checkNArg(2);
	56
	57
	58	}
	59
	60	void EvtBHadronic::decay( EvtParticle *p){
	61
	62	//added by Lange Jan4,2000
	63	static EvtId B0=EvtPDL::getId("B0");
	64	static EvtId D0=EvtPDL::getId("D0");
65	static EvtId DST0=EvtPDL::getId("D*0");
66	static EvtId D3P10=EvtPDL::getId("D'_10");
67	static EvtId D3P20=EvtPDL::getId("D_2*0");
68	static EvtId D3P00=EvtPDL::getId("D_0*0");
69	static EvtId D1P10=EvtPDL::getId("D_10");
70
71	static EvtISGW2FF ffmodel;
72
73	p->initializePhaseSpace(getNDaug(),getDaugs());
74
75	EvtVector4R p4[MAX_DAUG];
0ca57c2f	76	p->mass();
da0e9ce3	77
	78
	79	int i,j;
	80
	81	for ( i=0; i<getNDaug(); i++) {
	82	p4[i]=p->getDaug(i)->getP4();
	83	}
	84
	85	int bcurrent,wcurrent;
	86	int nbcurrent=0;
	87	int nwcurrent=0;
	88
	89	bcurrent=(int)getArg(0);
	90	wcurrent=(int)getArg(1);
	91
	92	EvtVector4C jb[5],jw[5];
	93	EvtTensor4C g,tds;
	94
	95	EvtVector4R p4b;
	96	p4b.set(p->mass(),0.0,0.0,0.0);
	97
	98	EvtVector4R q;
	99	double q2;
	100
	101	EvtComplex ep_meson_bb[5];
	102	double f,gf,ap,am;
	103	double fp,fm;
	104	double hf,kf,bp,bm;
	105	EvtVector4R pp,pm;
	106	EvtVector4C ep_meson_b[5];
	107
	108	switch (bcurrent) {
	109
	110	// D0
	111	case 1:
	112	q=p4b-p4[0];
	113	q2=q*q;
	114	nbcurrent=1;
	115	ffmodel.getscalarff(B0,D0,EvtPDL::getMeanMass(D0),
	116	EvtPDL::getMeanMass(getDaugs()[1]),&fp,&fm);
	117	jb[0]=EvtVector4C(fp(p4b+p4[0])-fmq);
	118	break;
	119	// D*
	120	case 2:
	121	q=p4b-p4[0];
	122	q2=q*q;
	123	nbcurrent=3;
	124	ffmodel.getvectorff(B0,DST0,EvtPDL::getMeanMass(DST0),q2,&f,&gf,&ap,&am);
	125
	126	g.setdiag(1.0,-1.0,-1.0,-1.0);
	127	tds = -f*g
0ca57c2f	128	-ap*(EvtGenFunctions::directProd(p4b,p4b)+EvtGenFunctions::directProd(p4b,p4[0]))
	129	-gfEvtComplex(0.0,1.0)dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
	130	-am*((EvtGenFunctions::directProd(p4b,p4b)-
	131	EvtGenFunctions::directProd(p4b,p4[0])));
da0e9ce3	132	jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
	133	jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
	134	jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
	135	break;
	136	// D1
	137	case 3:
	138	q=p4b-p4[0];
	139	q2=q*q;
	140	nbcurrent=3;
	141	ffmodel.getvectorff(B0,D3P10,EvtPDL::getMeanMass(D3P10),q2,&f,&gf,&ap,&am);
	142
	143	g.setdiag(1.0,-1.0,-1.0,-1.0);
	144	tds = -f*g
0ca57c2f	145	-ap*(EvtGenFunctions::directProd(p4b,p4b)+
	146	EvtGenFunctions::directProd(p4b,p4[0]))
	147	-gfEvtComplex(0.0,1.0)dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
	148	-am*((EvtGenFunctions::directProd(p4b,p4b)-
	149	EvtGenFunctions::directProd(p4b,p4[0])));
da0e9ce3	150	jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
	151	jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
	152	jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
	153
	154	break;
	155	// D2*
	156	case 4:
	157	q=p4b-p4[0];
	158	q2=q*q;
	159	nbcurrent=5;
	160	ffmodel.gettensorff(B0,D3P20,EvtPDL::getMeanMass(D3P20),q2,&hf,&kf,&bp,&bm);
	161	g.setdiag(1.0,-1.0,-1.0,-1.0);
	162
	163	ep_meson_b[0] = ((p->getDaug(0)->epsTensorParent(0)).cont2(p4b)).conj();
	164	ep_meson_b[1] = ((p->getDaug(0)->epsTensorParent(1)).cont2(p4b)).conj();
	165	ep_meson_b[2] = ((p->getDaug(0)->epsTensorParent(2)).cont2(p4b)).conj();
	166	ep_meson_b[3] = ((p->getDaug(0)->epsTensorParent(3)).cont2(p4b)).conj();
	167	ep_meson_b[4] = ((p->getDaug(0)->epsTensorParent(4)).cont2(p4b)).conj();
	168
	169	pp=p4b+p4[0];
	170	pm=p4b-p4[0];
	171
	172	ep_meson_bb[0]=ep_meson_b[0]*(p4b);
	173	ep_meson_bb[1]=ep_meson_b[1]*(p4b);
	174	ep_meson_bb[2]=ep_meson_b[2]*(p4b);
	175	ep_meson_bb[3]=ep_meson_b[3]*(p4b);
	176	ep_meson_bb[4]=ep_meson_b[4]*(p4b);
	177
0ca57c2f	178	jb[0]=EvtComplex(0.0,(1.0)hf)dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[0])
da0e9ce3	179	-kf*ep_meson_b[0]
	180	-bpep_meson_bb[0]pp-bmep_meson_bb[0]pm;
	181
0ca57c2f	182	jb[1]=EvtComplex(0.0,(1.0)hf)dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[1])
da0e9ce3	183	-kf*ep_meson_b[1]
	184	-bpep_meson_bb[1]pp-bmep_meson_bb[1]pm;
	185
0ca57c2f	186	jb[2]=EvtComplex(0.0,(1.0)hf)dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[2])
da0e9ce3	187	-kf*ep_meson_b[2]
	188	-bpep_meson_bb[2]pp-bmep_meson_bb[2]pm;
	189
0ca57c2f	190	jb[3]=EvtComplex(0.0,(1.0)hf)dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[3])
da0e9ce3	191	-kf*ep_meson_b[3]
	192	-bpep_meson_bb[3]pp-bmep_meson_bb[3]pm;
	193
0ca57c2f	194	jb[4]=EvtComplex(0.0,(1.0)hf)dual(EvtGenFunctions::directProd(pp,pm)).cont2(ep_meson_b[4])
da0e9ce3	195	-kf*ep_meson_b[4]
	196	-bpep_meson_bb[4]pp-bmep_meson_bb[4]pm;
	197	break;
	198	// D_0*
	199	case 5:
	200	q=p4b-p4[0];
	201	q2=q*q;
	202	double f,gf,ap,am;
	203	nbcurrent=3;
	204	ffmodel.getvectorff(B0,D1P10,EvtPDL::getMeanMass(D1P10),q2,&f,&gf,&ap,&am);
	205	g.setdiag(1.0,-1.0,-1.0,-1.0);
	206	tds = -f*g
0ca57c2f	207	-ap*(EvtGenFunctions::directProd(p4b,p4b)+
	208	EvtGenFunctions::directProd(p4b,p4[0]))
	209	+gfEvtComplex(0.0,1.0)dual(EvtGenFunctions::directProd(p4[0]+p4b,p4b-p4[0]))
	210	-am*((EvtGenFunctions::directProd(p4b,p4b)-
	211	EvtGenFunctions::directProd(p4b,p4[0])));
da0e9ce3	212	jb[0]=tds.cont1(p->getDaug(0)->epsParent(0).conj());
	213	jb[1]=tds.cont1(p->getDaug(0)->epsParent(1).conj());
	214	jb[2]=tds.cont1(p->getDaug(0)->epsParent(2).conj());
	215	break;
	216	// D_1
	217	case 6:
	218	q=p4b-p4[0];
	219	q2=q*q;
	220	nbcurrent=1;
	221	ffmodel.getscalarff(B0,D3P00,EvtPDL::getMeanMass(D3P00),q2,&fp,&fm);
	222	jb[0]=fp(p4b+p4[0])+fmq;
	223	break;
	224	default:
	225	report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown hadronic current."<<endl;
	226
	227	}
	228
	229	double norm;
	230
	231	switch (wcurrent) {
	232
	233
	234	case 1: case 3: case 4:
	235	nwcurrent=1;
	236	jw[0]=p4[getNDaug()-1];
	237	break;
	238
	239	case 2: case 5: case 6:
	240	nwcurrent=3;
	241	norm=1.0/sqrt(p4[1].get(0)*p4[1].get(0)/p4[1].mass2()-1.0);
	242	jw[0]=norm*p->getDaug(getNDaug()-1)->epsParent(0);
	243	jw[1]=norm*p->getDaug(getNDaug()-1)->epsParent(1);
	244	jw[2]=norm*p->getDaug(getNDaug()-1)->epsParent(2);
	245	break;
	246
	247
	248	default:
	249	report(ERROR,"EvtGen")<<"In EvtBHadronic, unknown W current."<<endl;
	250
	251	}
	252
	253	if (nbcurrent==1&&nwcurrent==1){
	254	vertex(jb[0]*jw[0]);
	255	return;
	256	}
	257
	258	if (nbcurrent==1){
	259	for(j=0;j<nwcurrent;j++){
	260	vertex(j,jb[0]*jw[j]);
	261	}
	262	return;
	263	}
	264
	265	if (nwcurrent==1){
	266	for(j=0;j<nbcurrent;j++){
	267	vertex(j,jb[j]*jw[0]);
	268	}
	269	return;
	270	}
	271
	272	for(j=0;j<nbcurrent;j++){
	273	for(i=0;i<nwcurrent;i++){
	274	vertex(j,i,jb[j]*jw[i]);
	275	}
276	}
277	return;
278
279	}
280
281
282
283
284
285
286