[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtBcVNpi.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: EvtBcVNpi.cc
//
// Description: Module to implement Bc -> psi + (n pi) decays
//
// Modification history:
//
//    AVL     July 6, 2012        Module created
//
//------------------------------------------------------------------------
// 
#include "EvtGenBase/EvtPatches.hh"
#include <iostream>
#include <iomanip>
#include <fstream>
#include <ctype.h>
#include <stdlib.h>
#include <string.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenModels/EvtTauHadnu.hh"
#include "EvtGenBase/EvtDiracSpinor.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtIdSet.hh"
#include "EvtGenBase/EvtParser.hh"

#include "EvtGenModels/EvtBcVNpi.hh"
#include "EvtGenModels/EvtWnPi.hh"


EvtBcVNpi::~EvtBcVNpi() {
//   cout<<"BcVNpi::destructor : nCall = "<<nCall<<" getProbMax(-1) = "<<getProbMax(-1)<<endl;
  
}

std::string EvtBcVNpi::getName(){ return "BC_VNPI";}

EvtDecayBase* EvtBcVNpi::clone() {  return new EvtBcVNpi;}

//======================================================
void EvtBcVNpi::init(){
    //cout<<"BcVNpi::init()"<<endl;
    
    checkNArg(1);
    checkSpinParent(EvtSpinType::SCALAR);
    checkSpinDaughter(0,EvtSpinType::VECTOR);
    for (int i=1; i<=(getNDaug()-1);i++) {
      checkSpinDaughter(i,EvtSpinType::SCALAR);
    };

    if(getNDaug()<2 || getNDaug()>6) {
      report(ERROR,"EvtGen") << "Have not yet implemented this final state in BcVNpi model" << endl;
      report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
      for ( int id=0; id<(getNDaug()-1); id++ ) 
	report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
      return;
    }

  
//     for(int i=0; i<getNDaug(); i++)
//       cout<<"BcVNpi::init \t\t daughter "<<i<<" : "<<getDaug(i).getId()<<"   "<<EvtPDL::name(getDaug(i)).c_str()<<endl;

   idVector = getDaug(0).getId();
    whichfit = int(getArg(0)+0.1);
//     cout<<"BcVNpi: whichfit ="<<whichfit<<"  idVector="<<idVector<<endl;
    ffmodel = new EvtBCVFF(idVector,whichfit);
    
    wcurr = new EvtWnPi();
    
    nCall = 0;
}

//======================================================
void EvtBcVNpi::initProbMax() {
//     cout<<"BcVNpi::initProbMax()"<<endl;
    if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==6) setProbMax(720000.);
    else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==6) setProbMax(471817.);
    else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==4) setProbMax(42000.);
    else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==4) setProbMax(16000.);
    
    else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==4) setProbMax(1200.);
    else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==4) setProbMax(2600.);
    else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==6) setProbMax(40000.);
    else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==6) setProbMax(30000.);
}

//======================================================
void EvtBcVNpi::decay( EvtParticle *root_particle ) {
   ++nCall;
//     cout<<"BcVNpi::decay()"<<endl;
    root_particle->initializePhaseSpace(getNDaug(),getDaugs());

    EvtVector4R
	    p4b(root_particle->mass(), 0., 0., 0.),                  // Bc momentum
	    p4meson=root_particle->getDaug(0)->getP4(),     		   // J/psi momenta
	    Q=p4b-p4meson;
    double Q2=Q.mass2();


// check pi-mesons and calculate hadronic current
    EvtVector4C hardCur;
//     bool foundHadCurr=false;
    if( getNDaug() == 2) {
      hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() );
//       foundHadCurr=true;
    }
    else if( getNDaug() == 3) {
      hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() , 
				 root_particle->getDaug(2)->getP4() 
			       );
//       foundHadCurr=true;   
    }
    else if( getNDaug() == 4) {
      hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() , 
				 root_particle->getDaug(2)->getP4(), 
				 root_particle->getDaug(3)->getP4() 
			       );
//       foundHadCurr=true;         
    }
    else if( getNDaug() == 6) // Bc -> psi pi+ pi+ pi- pi- pi+ from [Kuhn, Was, hep-ph/0602162
    {

		hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(),
					  root_particle->getDaug(2)->getP4(),
					  root_particle->getDaug(3)->getP4(),
					  root_particle->getDaug(4)->getP4(),
					  root_particle->getDaug(5)->getP4()
				 );
// 		foundHadCurr=true;
    }	
    else {
	    report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
	    report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
	    int id;
	    for ( id=0; id<(getNDaug()-1); id++ ) 
	    report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
	    ::abort();
    };  

// calculate Bc -> V W form-factors
	double a1f, a2f, vf, a0f;
	double m_meson = root_particle->getDaug(0)->mass();
	double m_b = root_particle->mass();
	ffmodel->getvectorff(root_particle->getId(),
				root_particle->getDaug(0)->getId(),
				Q2,
				m_meson,
				&a1f, 
				&a2f, 
				&vf, 
				&a0f);
	double a3f = ((m_b+m_meson)/(2.0*m_meson))*a1f -
	      ((m_b-m_meson)/(2.0*m_meson))*a2f;

// calculate Bc -> V W current
	EvtTensor4C H;
	H = a1f*(m_b+m_meson)*EvtTensor4C::g();
	H.addDirProd((-a2f/(m_b+m_meson))*p4b,p4b+p4meson);
	H+=EvtComplex(0.0,vf/(m_b+m_meson))*dual(EvtGenFunctions::directProd(p4meson+p4b,p4b-p4meson));
	H.addDirProd((a0f-a3f)*2.0*(m_meson/Q2)*p4b,p4b-p4meson);
	EvtVector4C  Heps=H.cont2(hardCur);
	
	for(int i=0; i<4; i++) {
		EvtVector4C  eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
		EvtComplex amp=eps*Heps;
		vertex(i,amp);
	};

}
Commit	Line	Data
0ca57c2f	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: EvtBcVNpi.cc
	12	//
	13	// Description: Module to implement Bc -> psi + (n pi) decays
	14	//
	15	// Modification history:
	16	//
	17	// AVL July 6, 2012 Module created
	18	//
	19	//------------------------------------------------------------------------
	20	//
	21	#include "EvtGenBase/EvtPatches.hh"
	22	#include <iostream>
	23	#include <iomanip>
	24	#include <fstream>
	25	#include <ctype.h>
	26	#include <stdlib.h>
	27	#include <string.h>
	28	#include "EvtGenBase/EvtParticle.hh"
	29	#include "EvtGenBase/EvtPDL.hh"
	30	#include "EvtGenBase/EvtGenKine.hh"
	31	#include "EvtGenModels/EvtTauHadnu.hh"
	32	#include "EvtGenBase/EvtDiracSpinor.hh"
	33	#include "EvtGenBase/EvtReport.hh"
	34	#include "EvtGenBase/EvtVector4C.hh"
	35	#include "EvtGenBase/EvtTensor4C.hh"
	36	#include "EvtGenBase/EvtIdSet.hh"
	37	#include "EvtGenBase/EvtParser.hh"
	38
	39	#include "EvtGenModels/EvtBcVNpi.hh"
	40	#include "EvtGenModels/EvtWnPi.hh"
	41
	42
	43
	44	EvtBcVNpi::~EvtBcVNpi() {
	45	// cout<<"BcVNpi::destructor : nCall = "<<nCall<<" getProbMax(-1) = "<<getProbMax(-1)<<endl;
	46
	47	}
	48
	49	std::string EvtBcVNpi::getName(){ return "BC_VNPI";}
	50
	51	EvtDecayBase* EvtBcVNpi::clone() { return new EvtBcVNpi;}
	52
	53	//======================================================
	54	void EvtBcVNpi::init(){
	55	//cout<<"BcVNpi::init()"<<endl;
	56
	57	checkNArg(1);
	58	checkSpinParent(EvtSpinType::SCALAR);
	59	checkSpinDaughter(0,EvtSpinType::VECTOR);
	60	for (int i=1; i<=(getNDaug()-1);i++) {
	61	checkSpinDaughter(i,EvtSpinType::SCALAR);
	62	};
	63
	64	if(getNDaug()<2 \|\| getNDaug()>6) {
65	report(ERROR,"EvtGen") << "Have not yet implemented this final state in BcVNpi model" << endl;
66	report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
67	for ( int id=0; id<(getNDaug()-1); id++ )
68	report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
69	return;
70	}
71
72
73	// for(int i=0; i<getNDaug(); i++)
74	// cout<<"BcVNpi::init \t\t daughter "<<i<<" : "<<getDaug(i).getId()<<" "<<EvtPDL::name(getDaug(i)).c_str()<<endl;
75
76	idVector = getDaug(0).getId();
77	whichfit = int(getArg(0)+0.1);
78	// cout<<"BcVNpi: whichfit ="<<whichfit<<" idVector="<<idVector<<endl;
79	ffmodel = new EvtBCVFF(idVector,whichfit);
80
81	wcurr = new EvtWnPi();
82
83	nCall = 0;
84	}
85
86	//======================================================
87	void EvtBcVNpi::initProbMax() {
88	// cout<<"BcVNpi::initProbMax()"<<endl;
89	if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==6) setProbMax(720000.);
90	else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==6) setProbMax(471817.);
91	else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 1 && getNDaug()==4) setProbMax(42000.);
92	else if(idVector == EvtPDL::getId("J/psi").getId() && whichfit == 2 && getNDaug()==4) setProbMax(16000.);
93
94	else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==4) setProbMax(1200.);
95	else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==4) setProbMax(2600.);
96	else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 1 && getNDaug()==6) setProbMax(40000.);
97	else if(idVector == EvtPDL::getId("psi(2S)").getId() && whichfit == 2 && getNDaug()==6) setProbMax(30000.);
98	}
99
100	//======================================================
101	void EvtBcVNpi::decay( EvtParticle *root_particle ) {
102	++nCall;
103	// cout<<"BcVNpi::decay()"<<endl;
104	root_particle->initializePhaseSpace(getNDaug(),getDaugs());
105
106	EvtVector4R
107	p4b(root_particle->mass(), 0., 0., 0.), // Bc momentum
108	p4meson=root_particle->getDaug(0)->getP4(), // J/psi momenta
109	Q=p4b-p4meson;
110	double Q2=Q.mass2();
111
112
113	// check pi-mesons and calculate hadronic current
114	EvtVector4C hardCur;
115	// bool foundHadCurr=false;
116	if( getNDaug() == 2) {
117	hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() );
118	// foundHadCurr=true;
119	}
120	else if( getNDaug() == 3) {
121	hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
122	root_particle->getDaug(2)->getP4()
123	);
124	// foundHadCurr=true;
125	}
126	else if( getNDaug() == 4) {
127	hardCur = wcurr->WCurrent( root_particle->getDaug(1)->getP4() ,
128	root_particle->getDaug(2)->getP4(),
129	root_particle->getDaug(3)->getP4()
130	);
131	// foundHadCurr=true;
132	}
133	else if( getNDaug() == 6) // Bc -> psi pi+ pi+ pi- pi- pi+ from [Kuhn, Was, hep-ph/0602162
134	{
135
136	hardCur = wcurr->WCurrent(root_particle->getDaug(1)->getP4(),
137	root_particle->getDaug(2)->getP4(),
138	root_particle->getDaug(3)->getP4(),
139	root_particle->getDaug(4)->getP4(),
140	root_particle->getDaug(5)->getP4()
141	);
142	// foundHadCurr=true;
143	}
144	else {
145	report(ERROR,"EvtGen") << "Have not yet implemented this final state in BCNPI model" << endl;
146	report(ERROR,"EvtGen") << "Ndaug="<<getNDaug() << endl;
147	int id;
148	for ( id=0; id<(getNDaug()-1); id++ )
149	report(ERROR,"EvtGen") << "Daug " << id << " "<<EvtPDL::name(getDaug(id)).c_str() << endl;
150	::abort();
151	};
152
153	// calculate Bc -> V W form-factors
154	double a1f, a2f, vf, a0f;
155	double m_meson = root_particle->getDaug(0)->mass();
156	double m_b = root_particle->mass();
157	ffmodel->getvectorff(root_particle->getId(),
158	root_particle->getDaug(0)->getId(),
159	Q2,
160	m_meson,
161	&a1f,
162	&a2f,
163	&vf,
164	&a0f);
165	double a3f = ((m_b+m_meson)/(2.0m_meson))a1f -
166	((m_b-m_meson)/(2.0m_meson))a2f;
167
168	// calculate Bc -> V W current
169	EvtTensor4C H;
170	H = a1f(m_b+m_meson)EvtTensor4C::g();
171	H.addDirProd((-a2f/(m_b+m_meson))*p4b,p4b+p4meson);
172	H+=EvtComplex(0.0,vf/(m_b+m_meson))*dual(EvtGenFunctions::directProd(p4meson+p4b,p4b-p4meson));
173	H.addDirProd((a0f-a3f)2.0(m_meson/Q2)*p4b,p4b-p4meson);
174	EvtVector4C Heps=H.cont2(hardCur);
175
176	for(int i=0; i<4; i++) {
177	EvtVector4C eps=root_particle->getDaug(0)->epsParent(i).conj(); // psi-meson polarization vector
178	EvtComplex amp=eps*Heps;
179	vertex(i,amp);
180	};
181
182	}
183