[u/mrichter/AliRoot.git] / TEvtGen / EvtGenModels / EvtYmSToYnSpipiCLEO.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: EvtGen/EvtYmSToYnSpipiCLEO.hh
//
// Description: This model is based on matrix element method used by
//              CLEO in Phys.Rev.D76:072001,2007 to model the dipion mass
//              and helicity angle distribution in the decays Y(mS) -> pi pi Y(nS),
//              where m,n are integers and m>n and m<4.
//              This model has two parameters, Re(B/A) and Im(B/A), which
//              are coefficients of the matrix element's terms determined by
//              the CLEO fits.
//
// Example:
//
// Decay  Upsilon(3S)
//  1.0000    Upsilon      pi+  pi-     YMSTOYNSPIPICLEO -2.523 1.189;
// Enddecay
// Decay  Upsilon(3S)
//  1.0000    Upsilon(2S)  pi+  pi-     YMSTOYNSPIPICLEO -0.395 0.001;
// Enddecay
// Decay  Upsilon(2S)
//  1.0000    Upsilon      pi+  pi-     YMSTOYNSPIPICLEO -0.753 0.000;
// Enddecay
//
//   --> the order of parameters is: Re(B/A) Im(B/A)
//
// Modification history:
//
//    SEKULA  Jan. 28, 2008         Module created
//
//------------------------------------------------------------------------


#include "EvtGenBase/EvtPatches.hh"
#include <stdlib.h>
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtGenKine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtVector4C.hh"
#include "EvtGenBase/EvtTensor4C.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenModels/EvtYmSToYnSpipiCLEO.hh"
#include <string>
using std::endl;

EvtYmSToYnSpipiCLEO::~EvtYmSToYnSpipiCLEO() {}

std::string EvtYmSToYnSpipiCLEO::getName(){

  return "YMSTOYNSPIPICLEO";     

}


EvtDecayBase* EvtYmSToYnSpipiCLEO::clone(){

  return new EvtYmSToYnSpipiCLEO;

}

void EvtYmSToYnSpipiCLEO::init(){

  static EvtId PIP=EvtPDL::getId("pi+");
  static EvtId PIM=EvtPDL::getId("pi-");
  static EvtId PI0=EvtPDL::getId("pi0");

  // check that there are 2 arguments
  checkNArg(2);
  checkNDaug(3);

  checkSpinParent(EvtSpinType::VECTOR);
  checkSpinDaughter(0,EvtSpinType::VECTOR);


  if ((!(getDaug(1)==PIP&&getDaug(2)==PIM))&&
      (!(getDaug(1)==PI0&&getDaug(2)==PI0))) {
    report(ERROR,"EvtGen") << "EvtYmSToYnSpipiCLEO generator expected "
                           << " pi+ and pi- (or pi0 and pi0) "
			   << "as 2nd and 3rd daughter. "<<endl;
    report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
    ::abort();
  }

}

void EvtYmSToYnSpipiCLEO::initProbMax() {
  setProbMax(2.0);
}      

void EvtYmSToYnSpipiCLEO::decay( EvtParticle *p){


  // We want to simulate the following process:
  //
  // Y(mS) -> Y(nS) X, X -> pi+ pi- (pi0 pi0)
  //
  // The CLEO analysis assumed such an intermediate process
  // were occurring, and wrote down the matrix element
  // and its components according to this assumption. 
  //
  //


  double ReB_over_A = getArg(0);
  double ImB_over_A = getArg(1);

  p->makeDaughters(getNDaug(),getDaugs());
  EvtParticle *v,*s1,*s2;
  v=p->getDaug(0);
  s1=p->getDaug(1);
  s2=p->getDaug(2);

  double m_pi = s1->getP4().mass();
  double M_mS = p->getP4().mass();
  double M_nS = v->getP4().mass();

// //   report(INFO,"EvtYmSToYnSpipiCLEO")  << "M_nS = " << v->getP4().mass() << endl;

  EvtVector4R P_nS;
  EvtVector4R P_pi1;
  EvtVector4R P_pi2;

  // Perform a simple accept/reject until we get a configuration of the
  // dipion system that passes
  bool acceptX = false;

  while( false == acceptX ) 
    {

      // Begin by generating a random X mass between the kinematic
      // boundaries, 2*m_pi and M(mS) - M(nS)

      double mX = EvtRandom::Flat(2.0 * m_pi, M_mS-M_nS);

      //   report(INFO,"EvtYmSToYnSpipiCLEO")  << "m_X = " << mX << endl;

      // Now create a two-body decay from the Y(mS) in its rest frame
      // of Y(mS) -> Y(nS) + X

      double masses[2];
      masses[0] = M_nS;
      masses[1] = mX;

      EvtVector4R p4[2];

      EvtGenKine::PhaseSpace( 2, masses, p4, M_mS );

      P_nS = p4[0];
      EvtVector4R P_X  = p4[1];

      // Now create the four-vectors for the two pions in the X
      // rest frame, X -> pi pi
  
      masses[0] = s1->mass();
      masses[1] = s2->mass();

      EvtGenKine::PhaseSpace( 2, masses, p4, P_X.mass() );

      // compute cos(theta), the polar helicity angle between a pi+ and
      // the direction opposite the Y(mS) in the X rest frame. If the pions are pi0s, then
      // choose the one where cos(theta) = [0:1].

      EvtVector4R P_YmS_X = boostTo(p->getP4(), P_X);
      double costheta = - p4[0].dot(P_YmS_X)/(p4[0].d3mag()*P_YmS_X.d3mag());
      if (EvtPDL::name(s1->getId()) == "pi0") {
	if (costheta < 0) {
	  costheta = - p4[1].dot(P_YmS_X)/(p4[1].d3mag()*P_YmS_X.d3mag());
	}
      }
      if (EvtPDL::name(s1->getId()) == "pi-") {
	costheta = - p4[1].dot(P_YmS_X)/(p4[1].d3mag()*P_YmS_X.d3mag());
      }
  
      // //   report(INFO,"EvtYmSToYnSpipiCLEO")  << "cos(theta) = " << costheta << endl;
    
    
      // Now boost the pion four vectors into the Y(mS) rest frame
      P_pi1 = boostTo(p4[0],P_YmS_X);
      P_pi2 = boostTo(p4[1],P_YmS_X);

      // Use a simple accept-reject to test this dipion system

      // Now compute the components of the matrix-element squared
      //
      // M(x,y)^2 = Q(x,y)^2 + |B/A|^2 * E1E2(x,y)^2 + 2*Re(B/A)*Q(x,y)*E1E2(x,y)
      //
      // x=m_pipi^2 and y = cos(theta), and where 
      //
      //   Q(x,y) = (x^2 + 2*m_pi^2)
      //  
      //   E1E2(x,y) = (1/4) * ( (E1 + E2)^2 - (E2 - E1)^2_max * cos(theta)^2 )
      //
      // and E1 + E2 = M_mS - M_nS and (E2 - E1)_max is the maximal difference
      // in the energy of the two pions allowed for a given mX value.
      //

      double Q    = (mX*mX - 2.0 * m_pi * m_pi);

      double deltaEmax = 
	- 2.0 * 
	sqrt( P_nS.get(0)*P_nS.get(0) - M_nS*M_nS ) *
	sqrt( 0.25 - pow(m_pi/mX,2.0));

      double sumE = (M_mS*M_mS - M_nS*M_nS + mX*mX)/(2.0 * M_mS);

      double E1E2 = 0.25 * ( pow(sumE, 2.0) - pow( deltaEmax * costheta, 2.0) );

      double M2 = Q*Q + (pow(ReB_over_A,2.0) + pow(ImB_over_A,2.0)) * E1E2*E1E2 + 2.0 * ReB_over_A * Q * E1E2;

      // phase space factor
      //
      // this is given as d(PS) = C * p(*)_X * p(X)_{pi+} * d(cosTheta) * d(m_X)
      // 
      // where C is a normalization constant, p(*)_X is the X momentum magnitude in the
      // Y(mS) rest frame, and p(X)_{pi+} is the pi+/pi0 momentum in the X rest frame
      //
  
      double dPS = 
	sqrt( (M_mS*M_mS - pow(M_nS + mX,2.0)) * (M_mS*M_mS - pow(M_nS - mX,2.0)) ) * // p(*)_X
	sqrt(mX*mX - 4*m_pi*m_pi); // p(X)_{pi}

      // the double-differential decay rate dG/(dcostheta dmX)
      double dG = M2 * dPS;

      // Throw a uniform random number from 0 --> probMax and do accept/reject on this
      
      double rnd = EvtRandom::Flat(0.0,getProbMax(0.0));

      if (rnd < dG)
	acceptX = true;

    }


  // initialize the daughters
  v->init(  getDaugs()[0], P_nS);
  s1->init( getDaugs()[1], P_pi1);
  s2->init( getDaugs()[2], P_pi2); 

//   report(INFO,"EvtYmSToYnSpipiCLEO")  << "M_nS = " << v->getP4().mass() << endl;
//   report(INFO,"EvtYmSToYnSpipiCLEO")  << "m_pi = " << s1->getP4().mass() << endl;
//   report(INFO,"EvtYmSToYnSpipiCLEO")  << "m_pi = " << s2->getP4().mass() << endl;
//   report(INFO,"EvtYmSToYnSpipiCLEO")  << "M2 = "   << M2 << endl;
  
  // Pass the polarization of the parent Upsilon
  EvtVector4C ep0,ep1,ep2;  
  
  ep0=p->eps(0);
  ep1=p->eps(1);
  ep2=p->eps(2);


  vertex(0,0,(ep0*v->epsParent(0).conj()));
  vertex(0,1,(ep0*v->epsParent(1).conj()));
  vertex(0,2,(ep0*v->epsParent(2).conj()));
  
  vertex(1,0,(ep1*v->epsParent(0).conj()));
  vertex(1,1,(ep1*v->epsParent(1).conj()));
  vertex(1,2,(ep1*v->epsParent(2).conj()));
  
  vertex(2,0,(ep2*v->epsParent(0).conj()));
  vertex(2,1,(ep2*v->epsParent(1).conj()));
  vertex(2,2,(ep2*v->epsParent(2).conj()));


  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: EvtGen/EvtYmSToYnSpipiCLEO.hh
	12	//
	13	// Description: This model is based on matrix element method used by
	14	// CLEO in Phys.Rev.D76:072001,2007 to model the dipion mass
	15	// and helicity angle distribution in the decays Y(mS) -> pi pi Y(nS),
	16	// where m,n are integers and m>n and m<4.
	17	// This model has two parameters, Re(B/A) and Im(B/A), which
	18	// are coefficients of the matrix element's terms determined by
	19	// the CLEO fits.
	20	//
	21	// Example:
	22	//
	23	// Decay Upsilon(3S)
	24	// 1.0000 Upsilon pi+ pi- YMSTOYNSPIPICLEO -2.523 1.189;
	25	// Enddecay
	26	// Decay Upsilon(3S)
	27	// 1.0000 Upsilon(2S) pi+ pi- YMSTOYNSPIPICLEO -0.395 0.001;
	28	// Enddecay
	29	// Decay Upsilon(2S)
	30	// 1.0000 Upsilon pi+ pi- YMSTOYNSPIPICLEO -0.753 0.000;
	31	// Enddecay
	32	//
	33	// --> the order of parameters is: Re(B/A) Im(B/A)
	34	//
	35	// Modification history:
	36	//
	37	// SEKULA Jan. 28, 2008 Module created
	38	//
	39	//------------------------------------------------------------------------
	40
	41
	42	#include "EvtGenBase/EvtPatches.hh"
	43	#include <stdlib.h>
	44	#include "EvtGenBase/EvtParticle.hh"
	45	#include "EvtGenBase/EvtGenKine.hh"
	46	#include "EvtGenBase/EvtPDL.hh"
	47	#include "EvtGenBase/EvtVector4C.hh"
	48	#include "EvtGenBase/EvtTensor4C.hh"
	49	#include "EvtGenBase/EvtReport.hh"
	50	#include "EvtGenBase/EvtRandom.hh"
	51	#include "EvtGenModels/EvtYmSToYnSpipiCLEO.hh"
	52	#include <string>
	53	using std::endl;
	54
	55	EvtYmSToYnSpipiCLEO::~EvtYmSToYnSpipiCLEO() {}
	56
	57	std::string EvtYmSToYnSpipiCLEO::getName(){
	58
	59	return "YMSTOYNSPIPICLEO";
	60
	61	}
	62
	63
	64	EvtDecayBase* EvtYmSToYnSpipiCLEO::clone(){
65
66	return new EvtYmSToYnSpipiCLEO;
67
68	}
69
70	void EvtYmSToYnSpipiCLEO::init(){
71
72	static EvtId PIP=EvtPDL::getId("pi+");
73	static EvtId PIM=EvtPDL::getId("pi-");
74	static EvtId PI0=EvtPDL::getId("pi0");
75
76	// check that there are 2 arguments
77	checkNArg(2);
78	checkNDaug(3);
79
80	checkSpinParent(EvtSpinType::VECTOR);
81	checkSpinDaughter(0,EvtSpinType::VECTOR);
82
83
84
85	if ((!(getDaug(1)==PIP&&getDaug(2)==PIM))&&
86	(!(getDaug(1)==PI0&&getDaug(2)==PI0))) {
87	report(ERROR,"EvtGen") << "EvtYmSToYnSpipiCLEO generator expected "
88	<< " pi+ and pi- (or pi0 and pi0) "
89	<< "as 2nd and 3rd daughter. "<<endl;
90	report(ERROR,"EvtGen") << "Will terminate execution!"<<endl;
91	::abort();
92	}
93
94	}
95
96	void EvtYmSToYnSpipiCLEO::initProbMax() {
97	setProbMax(2.0);
98	}
99
100	void EvtYmSToYnSpipiCLEO::decay( EvtParticle *p){
101
102
103	// We want to simulate the following process:
104	//
105	// Y(mS) -> Y(nS) X, X -> pi+ pi- (pi0 pi0)
106	//
107	// The CLEO analysis assumed such an intermediate process
108	// were occurring, and wrote down the matrix element
109	// and its components according to this assumption.
110	//
111	//
112
113
114	double ReB_over_A = getArg(0);
115	double ImB_over_A = getArg(1);
116
117	p->makeDaughters(getNDaug(),getDaugs());
118	EvtParticle v,s1,*s2;
119	v=p->getDaug(0);
120	s1=p->getDaug(1);
121	s2=p->getDaug(2);
122
123	double m_pi = s1->getP4().mass();
124	double M_mS = p->getP4().mass();
125	double M_nS = v->getP4().mass();
126
127	// // report(INFO,"EvtYmSToYnSpipiCLEO") << "M_nS = " << v->getP4().mass() << endl;
128
129	EvtVector4R P_nS;
130	EvtVector4R P_pi1;
131	EvtVector4R P_pi2;
132
133	// Perform a simple accept/reject until we get a configuration of the
134	// dipion system that passes
135	bool acceptX = false;
136
137	while( false == acceptX )
138	{
139
140	// Begin by generating a random X mass between the kinematic
141	// boundaries, 2*m_pi and M(mS) - M(nS)
142
143	double mX = EvtRandom::Flat(2.0 * m_pi, M_mS-M_nS);
144
145	// report(INFO,"EvtYmSToYnSpipiCLEO") << "m_X = " << mX << endl;
146
147	// Now create a two-body decay from the Y(mS) in its rest frame
148	// of Y(mS) -> Y(nS) + X
149
150	double masses[2];
151	masses[0] = M_nS;
152	masses[1] = mX;
153
154	EvtVector4R p4[2];
155
156	EvtGenKine::PhaseSpace( 2, masses, p4, M_mS );
157
158	P_nS = p4[0];
159	EvtVector4R P_X = p4[1];
160
161	// Now create the four-vectors for the two pions in the X
162	// rest frame, X -> pi pi
163
164	masses[0] = s1->mass();
165	masses[1] = s2->mass();
166
167	EvtGenKine::PhaseSpace( 2, masses, p4, P_X.mass() );
168
169	// compute cos(theta), the polar helicity angle between a pi+ and
170	// the direction opposite the Y(mS) in the X rest frame. If the pions are pi0s, then
171	// choose the one where cos(theta) = [0:1].
172
173	EvtVector4R P_YmS_X = boostTo(p->getP4(), P_X);
174	double costheta = - p4[0].dot(P_YmS_X)/(p4[0].d3mag()*P_YmS_X.d3mag());
175	if (EvtPDL::name(s1->getId()) == "pi0") {
176	if (costheta < 0) {
177	costheta = - p4[1].dot(P_YmS_X)/(p4[1].d3mag()*P_YmS_X.d3mag());
178	}
179	}
180	if (EvtPDL::name(s1->getId()) == "pi-") {
181	costheta = - p4[1].dot(P_YmS_X)/(p4[1].d3mag()*P_YmS_X.d3mag());
182	}
183
184	// // report(INFO,"EvtYmSToYnSpipiCLEO") << "cos(theta) = " << costheta << endl;
185
186
187
188	// Now boost the pion four vectors into the Y(mS) rest frame
189	P_pi1 = boostTo(p4[0],P_YmS_X);
190	P_pi2 = boostTo(p4[1],P_YmS_X);
191
192	// Use a simple accept-reject to test this dipion system
193
194	// Now compute the components of the matrix-element squared
195	//
196	// M(x,y)^2 = Q(x,y)^2 + \|B/A\|^2 * E1E2(x,y)^2 + 2Re(B/A)Q(x,y)*E1E2(x,y)
197	//
198	// x=m_pipi^2 and y = cos(theta), and where
199	//
200	// Q(x,y) = (x^2 + 2*m_pi^2)
201	//
202	// E1E2(x,y) = (1/4) * ( (E1 + E2)^2 - (E2 - E1)^2_max * cos(theta)^2 )
203	//
204	// and E1 + E2 = M_mS - M_nS and (E2 - E1)_max is the maximal difference
205	// in the energy of the two pions allowed for a given mX value.
206	//
207
208	double Q = (mXmX - 2.0 m_pi * m_pi);
209
210	double deltaEmax =
211	- 2.0 *
212	sqrt( P_nS.get(0)P_nS.get(0) - M_nSM_nS ) *
213	sqrt( 0.25 - pow(m_pi/mX,2.0));
214
215	double sumE = (M_mSM_mS - M_nSM_nS + mXmX)/(2.0 M_mS);
216
217	double E1E2 = 0.25 * ( pow(sumE, 2.0) - pow( deltaEmax * costheta, 2.0) );
218
219	double M2 = QQ + (pow(ReB_over_A,2.0) + pow(ImB_over_A,2.0)) E1E2E1E2 + 2.0 ReB_over_A * Q * E1E2;
220
221	// phase space factor
222	//
223	// this is given as d(PS) = C * p()_X p(X)_{pi+} * d(cosTheta) * d(m_X)
224	//
225	// where C is a normalization constant, p(*)_X is the X momentum magnitude in the
226	// Y(mS) rest frame, and p(X)_{pi+} is the pi+/pi0 momentum in the X rest frame
227	//
228
229	double dPS =
230	sqrt( (M_mSM_mS - pow(M_nS + mX,2.0)) (M_mSM_mS - pow(M_nS - mX,2.0)) ) // p(*)_X
231	sqrt(mXmX - 4m_pi*m_pi); // p(X)_{pi}
232
233	// the double-differential decay rate dG/(dcostheta dmX)
234	double dG = M2 * dPS;
235
236	// Throw a uniform random number from 0 --> probMax and do accept/reject on this
237
238	double rnd = EvtRandom::Flat(0.0,getProbMax(0.0));
239
240	if (rnd < dG)
241	acceptX = true;
242
243	}
244
245
246	// initialize the daughters
247	v->init( getDaugs()[0], P_nS);
248	s1->init( getDaugs()[1], P_pi1);
249	s2->init( getDaugs()[2], P_pi2);
250
251	// report(INFO,"EvtYmSToYnSpipiCLEO") << "M_nS = " << v->getP4().mass() << endl;
252	// report(INFO,"EvtYmSToYnSpipiCLEO") << "m_pi = " << s1->getP4().mass() << endl;
253	// report(INFO,"EvtYmSToYnSpipiCLEO") << "m_pi = " << s2->getP4().mass() << endl;
254	// report(INFO,"EvtYmSToYnSpipiCLEO") << "M2 = " << M2 << endl;
255
256	// Pass the polarization of the parent Upsilon
257	EvtVector4C ep0,ep1,ep2;
258
259	ep0=p->eps(0);
260	ep1=p->eps(1);
261	ep2=p->eps(2);
262
263
264	vertex(0,0,(ep0*v->epsParent(0).conj()));
265	vertex(0,1,(ep0*v->epsParent(1).conj()));
266	vertex(0,2,(ep0*v->epsParent(2).conj()));
267
268	vertex(1,0,(ep1*v->epsParent(0).conj()));
269	vertex(1,1,(ep1*v->epsParent(1).conj()));
270	vertex(1,2,(ep1*v->epsParent(2).conj()));
271
272	vertex(2,0,(ep2*v->epsParent(0).conj()));
273	vertex(2,1,(ep2*v->epsParent(1).conj()));
274	vertex(2,2,(ep2*v->epsParent(2).conj()));
275
276
277	return ;
278
279	}
280
281
282
283