BOSS_7.0.2: EvtDalitzReso Class Reference

Enumerator:

Undefined
PicPic
PizPiz
PiPi
KcKc
KzKz
KK
EtaPic
EtaPiz
PicPicKK
WA76

Definition at line 51 of file EvtDalitzReso.hh.

00051 {Undefined=0,PicPic=1,PizPiz,PiPi,KcKc,KzKz,KK,EtaPic,EtaPiz,PicPicKK,WA76};

enum EvtDalitzReso::NumType

Enumerator:

NBW
RBW_ZEMACH
RBW_KUEHN
RBW_CLEO
RBW_ZEMACH2
GS_CLEO
K_MATRIX
RBW_CLEO_ZEMACH
GS_CLEO_ZEMACH
LASS
K_MATRIX_I
K_MATRIX_II
GAUSS_CLEO
GAUSS_CLEO_ZEMACH

Definition at line 38 of file EvtDalitzReso.hh.

00038                { NBW            = 0 , RBW_ZEMACH        = 1 , RBW_KUEHN  = 2 , RBW_CLEO        = 3 ,
00039                  RBW_ZEMACH2    = 4 , GS_CLEO           = 5 , K_MATRIX   = 6 , RBW_CLEO_ZEMACH = 7 ,
00040                  GS_CLEO_ZEMACH = 8 , LASS              = 9 , K_MATRIX_I = 10, K_MATRIX_II     = 11,
00041                  GAUSS_CLEO     = 12, GAUSS_CLEO_ZEMACH = 13 };

EvtDalitzReso::EvtDalitzReso ( ) [inline]

Definition at line 53 of file EvtDalitzReso.hh.

Referenced by clone().

00053 {};

EvtDalitzReso::EvtDalitzReso	(	const EvtDalitzPlot &	dp,
		EvtCyclic3::Pair	pairAng,
		EvtCyclic3::Pair	pairRes,
		EvtSpinType::spintype	spin,
		double	m0,
		double	g0,
		NumType	typeN
	)

EvtDalitzReso::EvtDalitzReso	(	const EvtDalitzPlot &	dp,
		EvtCyclic3::Pair	pairAng,
		EvtCyclic3::Pair	pairRes,
		EvtSpinType::spintype	spin,
		double	m0,
		double	g0,
		NumType	typeN,
		double	m0_mix,
		double	g0_mix,
		double	delta_mix,
		EvtComplex	amp_mix
	)

EvtDalitzReso::EvtDalitzReso	(	const EvtDalitzPlot &	dp,
		EvtCyclic3::Pair	pairAng,
		EvtCyclic3::Pair	pairRes,
		EvtSpinType::spintype	spin,
		double	m0,
		NumType	typeN,
		double	g1,
		double	g2,
		CouplingType	coupling2
	)

EvtDalitzReso::EvtDalitzReso	(	const EvtDalitzPlot &	dp,
		EvtCyclic3::Pair	pairRes,
		std::string	nameIndex,
		NumType	typeN,
		EvtComplex	fr12prod,
		EvtComplex	fr13prod,
		EvtComplex	fr14prod,
		EvtComplex	fr15prod,
		double	s0prod
	)

EvtDalitzReso::EvtDalitzReso	(	const EvtDalitzPlot &	dp,
		EvtCyclic3::Pair	pairRes,
		double	m0,
		double	g0,
		double	a,
		double	r,
		double	B,
		double	phiB,
		double	R,
		double	phiR
	)

EvtDalitzReso::EvtDalitzReso ( const EvtDalitzReso & other )

Definition at line 155 of file EvtDalitzReso.cc.

00156   : _dp(other._dp),
00157     _pairAng(other._pairAng),
00158     _pairRes(other._pairRes),
00159     _spin(other._spin),
00160     _typeN(other._typeN),
00161     _m0(other._m0),_g0(other._g0),
00162     _vb(other._vb),_vd(other._vd),
00163     _massFirst(other._massFirst),_massSecond(other._massSecond),
00164     _m0_mix(other._m0_mix),_g0_mix(other._g0_mix),_delta_mix(other._delta_mix),_amp_mix(other._amp_mix),
00165     _g1(other._g1),_g2(other._g2),_coupling2(other._coupling2),
00166     _kmatrix_index(other._kmatrix_index),
00167     _fr12prod(other._fr12prod),_fr13prod(other._fr13prod),_fr14prod(other._fr14prod),_fr15prod(other._fr15prod),
00168     _s0prod(other._s0prod),
00169     _a(other._a),_r(other._r),_Blass(other._Blass),_phiB(other._phiB),_R(other._R),_phiR(other._phiR)
00170 {}

EvtDalitzReso::~EvtDalitzReso ( )

Definition at line 173 of file EvtDalitzReso.cc.

00174 {}

double EvtDalitzReso::angDep ( const EvtDalitzPoint & p ) [private]

Definition at line 468 of file EvtDalitzReso.cc.

References _pairAng, _spin, EvtdFunction::d(), EvtSpinType::getSpin2(), Bes_Common::INFO, report(), and x.

Referenced by numerator().

00469 { 
00470   // Angular dependece for factorizable amplitudes  
00471   // unphysical cosines indicate we are in big trouble
00472   double cosTh = x.cosTh(_pairAng,_pairRes);  // angle between common(reso,ang) and other(reso)
00473   if(fabs(cosTh) > 1.) {
00474     report(INFO,"EvtGen") << "cosTh " << cosTh << std::endl; 
00475     assert(0);
00476   }
00477   
00478   // in units of half-spin
00479   return EvtdFunction::d(EvtSpinType::getSpin2(_spin),2*0,2*0,acos(cosTh));
00480 }

EvtDalitzReso* EvtDalitzReso::clone ( ) const [inline]

Definition at line 77 of file EvtDalitzReso.hh.

References EvtDalitzReso().

00077 { return new EvtDalitzReso(*this); }

EvtComplex EvtDalitzReso::evaluate ( const EvtDalitzPoint & p )

Definition at line 177 of file EvtDalitzReso.cc.

References _coupling2, _dp, _g0, _g0_mix, _g1, _g2, _kmatrix_index, _m0, _m0_mix, _massFirst, _massSecond, _spin, _typeN, _vb, _vd, EvtDalitzPlot::bigM(), EtaPic, EtaPiz, EvtTwoBodyVertex::formFactor(), Fvector(), GAUSS_CLEO, GAUSS_CLEO_ZEMACH, EvtPDL::getId(), EvtPDL::getMass(), GS_CLEO, GS_CLEO_ZEMACH, K_MATRIX, K_MATRIX_I, K_MATRIX_II, KcKc, KK, KzKz, lass(), LASS, EvtDalitzPlot::m(), mixFactor(), NBW, numerator(), EvtCyclic3::other(), EvtTwoBodyVertex::pD(), PicPic, PicPicKK, PiPi, PizPiz, propBreitWigner(), propBreitWignerRel(), propBreitWignerRelCoupled(), propGauss(), propGounarisSakurai(), psFactor(), Undefined, WA76, EvtTwoBodyVertex::widthFactor(), and x.

Referenced by EvtD0mixDalitz::dalitzKsKK(), and EvtD0mixDalitz::dalitzKsPiPi().

00178 {
00179   double m = sqrt(x.q(_pairRes));
00180 
00181   // do use always hash table (speed up fitting)
00182   if (_typeN==K_MATRIX || _typeN==K_MATRIX_I || _typeN==K_MATRIX_II)
00183     return Fvector( m*m, _kmatrix_index );
00184 
00185   if (_typeN==LASS)
00186     return lass(m*m);
00187 
00188   EvtComplex amp(1.0,0.0);
00189 
00190   if (_dp.bigM() != x.bigM()) _vb = EvtTwoBodyVertex(_m0,_dp.m(EvtCyclic3::other(_pairRes)),x.bigM(),_spin); 
00191   EvtTwoBodyKine vb(m,x.m(EvtCyclic3::other(_pairRes)),x.bigM());
00192   EvtTwoBodyKine vd(_massFirst,_massSecond,m);   
00193 
00194   EvtComplex prop(0,0);
00195   if (_typeN==NBW) {
00196     prop = propBreitWigner(_m0,_g0,m);
00197   } else if (_typeN==GAUSS_CLEO || _typeN==GAUSS_CLEO_ZEMACH) {
00198     prop = propGauss(_m0,_g0,m);
00199   } else {
00200     if (_coupling2==Undefined) {  
00201       // single BW
00202       double g = (_g0<=0. || _vd.pD()<=0.)? -_g0 : _g0*_vd.widthFactor(vd);  // running width
00203       if (_typeN==GS_CLEO || _typeN==GS_CLEO_ZEMACH) {
00204         // Gounaris-Sakurai (GS)
00205         assert(_massFirst==_massSecond);
00206         prop = propGounarisSakurai(_m0,fabs(_g0),_vd.pD(),m,g,vd.p());
00207       } else {
00208         // standard relativistic BW
00209         prop = propBreitWignerRel(_m0,g,m);
00210       }
00211     } else {    
00212       // coupled width BW
00213       EvtComplex G1,G2;
00214       switch (_coupling2) { 
00215       case PicPic: {
00216         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00217         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00218         G2 = _g2*_g2*psFactor(mPic,mPic,m);
00219         break;
00220       }
00221       case PizPiz: {
00222         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00223         static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
00224         G2 = _g2*_g2*psFactor(mPiz,mPiz,m);
00225         break;
00226       }
00227       case PiPi: {
00228         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00229         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00230         static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
00231         G2 = _g2*_g2*psFactor(mPic,mPic,mPiz,mPiz,m);
00232         break;
00233       }
00234       case KcKc: {
00235         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00236         static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
00237         G2 = _g2*_g2*psFactor(mKc,mKc,m);
00238         break;
00239       }
00240       case KzKz: {
00241         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00242         static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
00243         G2 = _g2*_g2*psFactor(mKz,mKz,m);
00244         break;
00245       }
00246       case KK: {
00247         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00248         static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
00249         static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
00250         G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
00251         break;
00252       }
00253       case EtaPic: {
00254         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00255         static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
00256         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00257         G2 = _g2*_g2*psFactor(mEta,mPic,m);
00258         break;
00259       }
00260       case EtaPiz: {
00261         G1 = _g1*_g1*psFactor(_massFirst,_massSecond,m);
00262         static double mEta = EvtPDL::getMass( EvtPDL::getId( "eta" ) );
00263         static double mPiz = EvtPDL::getMass( EvtPDL::getId( "pi0" ) );
00264         G2 = _g2*_g2*psFactor(mEta,mPiz,m);
00265         break;
00266       }
00267       case PicPicKK: {
00268         static double mPic = EvtPDL::getMass( EvtPDL::getId( "pi+" ) );
00269         //G1 = _g1*_g1*psFactor(mPic,mPic,m);
00270         G1 = _g1*psFactor(mPic,mPic,m);
00271         static double mKc = EvtPDL::getMass( EvtPDL::getId( "K+" ) );
00272         static double mKz = EvtPDL::getMass( EvtPDL::getId( "K0" ) );
00273         //G2 = _g2*_g2*psFactor(mKc,mKc,mKz,mKz,m);
00274         G2 = _g2*psFactor(mKc,mKc,mKz,mKz,m);
00275         break;
00276       }
00277       default:
00278         std::cout << "EvtDalitzReso:evaluate(): PANIC, wrong coupling2 state." << std::endl;
00279         assert(0);
00280         break;
00281       }
00282       // calculate standard couple BW propagator
00283       if (_coupling2 != WA76)
00284         prop = _g1*propBreitWignerRelCoupled(_m0,G1,G2,m);
00285     } 
00286   }
00287   amp *= prop;
00288 
00289   // Compute form-factors (Blatt-Weisskopf penetration factor)
00290   amp *= _vb.formFactor(vb);  
00291   amp *= _vd.formFactor(vd);  
00292 
00293   // Compute numerator (angular distribution)
00294   amp *= numerator(x,vb,vd);  
00295 
00296   // Compute electromagnetic mass mixing factor
00297   if (_m0_mix>0.) {
00298     EvtComplex prop_mix;
00299     if (_typeN==NBW) {
00300       prop_mix = propBreitWigner(_m0_mix,_g0_mix,m);
00301     } else {
00302       assert(_g1<0.); // running width only
00303       double g_mix = _g0_mix*_vd.widthFactor(vd);
00304       prop_mix = propBreitWignerRel(_m0_mix,g_mix,m);
00305     }
00306     amp *= mixFactor(prop,prop_mix);
00307   }
00308 
00309   return amp;
00310 }

EvtComplex EvtDalitzReso::Fvector	(	double	s,
		int	index
	)			`[private]`

Definition at line 491 of file EvtDalitzReso.cc.

References _fr12prod, _fr13prod, _fr14prod, _fr15prod, _s0prod, _typeN, EvtCyclic3::A, EvtCyclic3::C, F, genRecEmupikp::i, EvtMatrix< T >::inverse(), ganga-rec::j, K_MATRIX, K_MATRIX_I, K_MATRIX_II, meta, mpi, PRECISION, and EvtMatrix< T >::setRange().

Referenced by evaluate().

00492 {
00493   assert(index>=1 && index<=6);
00494 
00495   //Define the complex coupling constant
00496   //The convection is as follow
00497   //i=0 --> pi+ pi-
00498   //i=1 --> KK
00499   //i=2 --> 4pi
00500   //i=3 --> eta eta
00501   //i=4 --> eta eta'
00502   //The first index is the resonace-pole index
00503       
00504   double g[5][5]; // Coupling constants. The first index is the pole index. The second index is the decay channel
00505   double ma[5];   // Pole masses. The unit is in GeV
00506 
00507   int solution = (_typeN==K_MATRIX)? 3 : (    (_typeN==K_MATRIX_I)? 1 : ( (_typeN==K_MATRIX_II)? 2 : 0 )    ) ;
00508   if (solution==0) { std::cout << "EvtDalitzReso::Fvector() error. Kmatrix solution incorrectly chosen ! " << std::endl; abort(); } 
00509 
00510   if (solution == 3 ) {
00511 
00512     // coupling constants
00513     //pi+pi- channel
00514     g[0][0]=0.22889;
00515     g[1][0]=0.94128;
00516     g[2][0]=0.36856;
00517     g[3][0]=0.33650;
00518     g[4][0]=0.18171;
00519     //K+K- channel
00520     g[0][1]=-0.55377;
00521     g[1][1]=0.55095;
00522     g[2][1]=0.23888;
00523     g[3][1]=0.40907;
00524     g[4][1]=-0.17558;
00525     //4pi channel
00526     g[0][2]=0;
00527     g[1][2]=0;
00528     g[2][2]=0.55639;
00529     g[3][2]=0.85679;
00530     g[4][2]=-0.79658;
00531     //eta eta channel
00532     g[0][3]=-0.39899;
00533     g[1][3]=0.39065;
00534     g[2][3]=0.18340;
00535     g[3][3]=0.19906;
00536     g[4][3]=-0.00355;
00537     //eta eta' channel
00538     g[0][4]=-0.34639;
00539     g[1][4]=0.31503;
00540     g[2][4]=0.18681;
00541     g[3][4]=-0.00984;
00542     g[4][4]=0.22358;
00543 
00544     // Pole masses
00545     ma[0]=0.651;      
00546     ma[1]=1.20360;
00547     ma[2]=1.55817;
00548     ma[3]=1.21000;
00549     ma[4]=1.82206;
00550 
00551   } else if (solution == 1) { // solnI.txt 
00552     
00553     // coupling constants
00554     //pi+pi- channel
00555     g[0][0]=0.31896;
00556     g[1][0]=0.85963;
00557     g[2][0]=0.47993;
00558     g[3][0]=0.45121;
00559     g[4][0]=0.39391;
00560     //K+K- channel
00561     g[0][1]=-0.49998;
00562     g[1][1]=0.52402;
00563     g[2][1]=0.40254;
00564     g[3][1]=0.42769;
00565     g[4][1]=-0.30860;
00566     //4pi channel
00567     g[0][2]=0;
00568     g[1][2]=0;
00569     g[2][2]=1.0;
00570     g[3][2]=1.15088;
00571     g[4][2]=0.33999;
00572     //eta eta channel
00573     g[0][3]=-0.21554;
00574     g[1][3]=0.38093;
00575     g[2][3]=0.21811;
00576     g[3][3]=0.22925;
00577     g[4][3]=0.06919;
00578     //eta eta' channel
00579     g[0][4]=-0.18294;
00580     g[1][4]=0.23788;
00581     g[2][4]=0.05454;
00582     g[3][4]=0.06444;
00583     g[4][4]=0.32620;
00584 
00585     // Pole masses
00586     ma[0]=0.7369;
00587     ma[1]=1.24347;
00588     ma[2]=1.62681;
00589     ma[3]=1.21900;
00590     ma[4]=1.74932;
00591 
00592   } else if (solution == 2) { // solnIIa.txt 
00593     
00594     // coupling constants
00595     //pi+pi- channel
00596     g[0][0]=0.26014;
00597     g[1][0]=0.95289;
00598     g[2][0]=0.46244;
00599     g[3][0]=0.41848;
00600     g[4][0]=0.01804;
00601     //K+K- channel
00602     g[0][1]=-0.57849;
00603     g[1][1]=0.55887;
00604     g[2][1]=0.31712;
00605     g[3][1]=0.49910;
00606     g[4][1]=-0.28430;
00607     //4pi channel
00608     g[0][2]=0;
00609     g[1][2]=0;
00610     g[2][2]=0.70340;
00611     g[3][2]=0.96819;
00612     g[4][2]=-0.90100;
00613     //eta eta channel
00614     g[0][3]=-0.32936;
00615     g[1][3]=0.39910;
00616     g[2][3]=0.22963;
00617     g[3][3]=0.24415;
00618     g[4][3]=-0.07252;
00619     //eta eta' channel
00620     g[0][4]=-0.30906;
00621     g[1][4]=0.31143;
00622     g[2][4]=0.19802;
00623     g[3][4]=-0.00522;
00624     g[4][4]=0.17097;
00625 
00626     // Pole masses
00627     ma[0]=0.67460;
00628     ma[1]=1.21094;
00629     ma[2]=1.57896;
00630     ma[3]=1.21900;
00631     ma[4]=1.86602;
00632   } 
00633 
00634   //Now define the K-matrix pole
00635   double  rho1sq,rho2sq,rho4sq,rho5sq;
00636   EvtComplex rho[5];
00637   double f[5][5];
00638 
00639   //Initalize the mass of the resonance
00640   double mpi=0.13957;
00641   double mK=0.493677;     //using charged K value
00642   double meta=0.54775;    //using PDG value
00643   double metap=0.95778;   //using PDG value
00644     
00645   //Initialize the matrix to value zero
00646   EvtComplex K[5][5];
00647   for(int i=0;i<5;i++) { 
00648     for(int j=0;j<5;j++) {
00649       K[i][j]=EvtComplex(0,0);
00650       f[i][j]=0;
00651     }
00652   }
00653 
00654   //Input the _f[i][j] scattering data
00655   double s_scatt=0.0 ; 
00656   if (solution == 3) 
00657     s_scatt=-3.92637; 
00658   else if (solution == 1) 
00659     s_scatt= -5.0 ;
00660   else if (solution == 2) 
00661     s_scatt= -5.0 ; 
00662   double sa=1.0;
00663   double sa_0=-0.15;
00664   if (solution == 3) {
00665     f[0][0]=0.23399;  // f^scatt
00666     f[0][1]=0.15044;
00667     f[0][2]=-0.20545;
00668     f[0][3]=0.32825;
00669     f[0][4]=0.35412;
00670   }else if (solution == 1) {
00671     f[0][0]=0.04214;  // f^scatt
00672     f[0][1]=0.19865;
00673     f[0][2]=-0.63764;
00674     f[0][3]=0.44063;
00675     f[0][4]=0.36717;
00676   }else if (solution == 2) {
00677     f[0][0]=0.26447;  // f^scatt
00678     f[0][1]=0.10400;
00679     f[0][2]=-0.35445;
00680     f[0][3]=0.31596;
00681     f[0][4]=0.42483;
00682   }   
00683   f[1][0]=f[0][1];
00684   f[2][0]=f[0][2];
00685   f[3][0]=f[0][3];
00686   f[4][0]=f[0][4];
00687 
00688   //Now construct the phase-space factor
00689   //For eta-eta' there is no difference term
00690   rho1sq = 1. - pow( mpi + mpi, 2 ) / s;   //pi+ pi- phase factor
00691   if( rho1sq >= 0 )
00692     rho[ 0 ] = EvtComplex( sqrt( rho1sq ), 0 );
00693   else
00694     rho[ 0 ] = EvtComplex( 0, sqrt( -rho1sq ) );  
00695 
00696   rho2sq = 1. - pow( mK + mK, 2 ) / s;
00697   if( rho2sq >= 0 )
00698     rho[ 1 ] = EvtComplex( sqrt( rho2sq ), 0 );
00699   else
00700     rho[ 1 ] = EvtComplex( 0, sqrt( -rho2sq ) );
00701 
00702   //using the A&S 4pi phase space Factor:
00703   //Shit, not continue
00704   if( s <= 1 )
00705     {
00706       double real   = 1.2274 + .00370909 / ( s * s ) - .111203 / s - 6.39017 * s + 16.8358*s*s - 21.8845*s*s*s + 11.3153*s*s*s*s;
00707       double cont32 = sqrt(1.0-(16.0*mpi*mpi));
00708       rho[ 2 ] = EvtComplex( cont32 * real, 0 );
00709     }
00710   else
00711     rho[ 2 ] = EvtComplex( sqrt( 1. - 16. * mpi * mpi / s ), 0 );
00712 
00713   rho4sq = 1. - pow( meta + meta, 2 ) / s;
00714   if( rho4sq >= 0 )
00715     rho[ 3 ] = EvtComplex( sqrt( rho4sq ), 0 );
00716   else
00717     rho[ 3 ] = EvtComplex( 0, sqrt( -rho4sq ) );
00718 
00719   rho5sq = 1. - pow( meta + metap, 2 ) / s;
00720   if( rho5sq >= 0 )
00721     rho[ 4 ] = EvtComplex( sqrt( rho5sq ), 0 );
00722   else
00723     rho[ 4 ] = EvtComplex( 0, sqrt( -rho5sq ) );
00724 
00725   double smallTerm = 1; // Factor to prevent divergences.
00726 
00727   // Check if some pole may arise problems.
00728   for ( int pole = 0; pole < 5; pole++ )
00729     if ( fabs( pow( ma[ pole ], 2 ) - s ) < PRECISION )
00730       smallTerm = pow( ma[ pole ], 2 ) - s;
00731 
00732   //now sum all the pole
00733   //equation (3) in the E791 K-matrix paper
00734   for(int i=0;i<5;i++) { 
00735     for(int j=0;j<5;j++) {  
00736       for (int pole_index=0;pole_index<5;pole_index++) {
00737         double A=g[pole_index][i]*g[pole_index][j];
00738         double B=ma[pole_index]*ma[pole_index]-s;
00739 
00740         if ( fabs( B ) < PRECISION )
00741           K[ i ][ j ] += EvtComplex( A    , 0 );
00742         else
00743           K[ i ][ j ] += EvtComplex( A / B, 0 ) * smallTerm;
00744       }
00745     }
00746   }
00747 
00748   //now add the SVT part
00749   for(int i=0;i<5;i++) { 
00750     for(int j=0;j<5;j++) {
00751       double C=f[i][j]*(1.0-s_scatt);
00752       double D=(s-s_scatt);
00753       K[ i ][ j ] += EvtComplex( C / D, 0 ) * smallTerm;
00754     }
00755   }
00756 
00757   //Fix the bug in the FOCUS paper
00758   //Include the Alder zero term:
00759   for(int i=0;i<5;i++) { 
00760     for(int j=0;j<5;j++) {
00761       double E=(s-(sa*mpi*mpi*0.5))*(1.0-sa_0);
00762       double F=(s-sa_0);    
00763       K[ i ][ j ] *= EvtComplex(E/F,0);
00764     }
00765   }
00766 
00767   //This is not correct!
00768   //(1-ipK) != (1-iKp)
00769   static EvtMatrix< EvtComplex > mat;
00770   mat.setRange( 5 ); // Try to do in only the first time. DEFINE ALLOCATION IN CONSTRUCTOR.
00771 
00772   for ( int row = 0; row < 5; row++ )
00773     for ( int col = 0; col < 5; col++ )
00774       mat( row, col ) = ( row == col ) * smallTerm - EvtComplex( 0., 1. ) * K[ row ][ col ] * rho[ col ];
00775 
00776 
00777   EvtMatrix< EvtComplex >* matInverse = mat.inverse();  //The 1st row of the inverse matrix. This matrix is {(I-iKp)^-1}_0j
00778   vector< EvtComplex > U1j;
00779   for ( int j = 0; j < 5; j++ )
00780     U1j.push_back( (*matInverse)[ 0 ][ j ] );
00781 
00782   delete matInverse;
00783 
00784   //this calculates final F0 factor
00785   EvtComplex value( 0, 0 );
00786   if (index<=5) {
00787     //this calculates the beta_idx Factors
00788     for(int j=0;j<5;j++) {        // sum for 5 channel
00789       EvtComplex top    = U1j[j]*g[index-1][j];
00790       double     bottom = ma[index-1]*ma[index-1]-s;
00791 
00792       if ( fabs( bottom ) < PRECISION )
00793         value += top;
00794       else
00795         value += top / bottom * smallTerm;
00796     }
00797   } else {
00798     //this calculates fprod Factors
00799     value += U1j[0];
00800     value += U1j[1]*_fr12prod;
00801     value += U1j[2]*_fr13prod;
00802     value += U1j[3]*_fr14prod;
00803     value += U1j[4]*_fr15prod;
00804 
00805     value *= (1-_s0prod)/(s-_s0prod) * smallTerm;
00806   }
00807 
00808   return value;
00809 }

double EvtDalitzReso::GS_d	(	const double &	m0,
		const double &	k0
	)			`[inline, private]`

Definition at line 393 of file EvtDalitzReso.cc.

References _massFirst, and EvtConst::pi.

Referenced by propGounarisSakurai().

00394 {return 3./EvtConst::pi*_massFirst*_massFirst/(k0*k0)*log((m0+2.*k0)/(2.*_massFirst)) + 
00395    m0/(2.*EvtConst::pi*k0) - _massFirst*_massFirst*m0/(EvtConst::pi*k0*k0*k0) ;}

double EvtDalitzReso::GS_dhods	(	const double &	m0,
		const double &	k0
	)			`[inline, private]`

Definition at line 390 of file EvtDalitzReso.cc.

References GS_h(), and EvtConst::pi.

Referenced by GS_f().

00391 {return GS_h(m0,k0)*( 0.125/(k0*k0) - 0.5/(m0*m0) ) + 0.5/(EvtConst::pi*m0*m0) ;}

double EvtDalitzReso::GS_f	(	const double &	m0,
		const double &	g0,
		const double &	k0,
		const double &	m,
		const double &	k
	)			`[inline, private]`

Definition at line 378 of file EvtDalitzReso.cc.

References GS_dhods(), and GS_h().

Referenced by propGounarisSakurai().

00379 {
00380   // m: sqrt(s)
00381   // m0: nominal resonance mass
00382   // k: momentum of pion in resonance rest frame (at m)
00383   // k0: momentum of pion in resonance rest frame (at nominal resonance mass)
00384   return g0*m0*m0/(k0*k0*k0)*( k*k*(GS_h(m,k)-GS_h(m0,k0)) + (m0*m0-m*m)*k0*k0*GS_dhods(m0,k0) );
00385 }

double EvtDalitzReso::GS_h	(	const double &	m,
		const double &	k
	)			`[inline, private]`

Definition at line 387 of file EvtDalitzReso.cc.

References _massFirst, and EvtConst::pi.

Referenced by GS_dhods(), and GS_f().

00388 {return 2./EvtConst::pi*k/m*log((m+2.*k)/(2.*_massFirst)) ;}

EvtComplex EvtDalitzReso::lass ( double s ) [private]

Definition at line 813 of file EvtDalitzReso.cc.

References _a, _Blass, _g0, _m0, _massFirst, _massSecond, _phiB, _phiR, _R, _r, _vd, cos(), EvtTwoBodyKine::p(), q, sin(), T, and EvtTwoBodyVertex::widthFactor().

Referenced by evaluate().

00814 {
00815   EvtTwoBodyKine vd(_massFirst,_massSecond, sqrt(s));
00816   double q = vd.p();
00817   double GammaM = _g0*_vd.widthFactor(vd);  // running width;
00818 
00819   //calculate the background phase motion
00820   double cot_deltaB = 1.0/(_a*q) + 0.5*_r*q;
00821   double deltaB = atan( 1.0/cot_deltaB);
00822   double totalB = deltaB + _phiB ;
00823 
00824   //calculate the resonant phase motion
00825   double deltaR = atan((_m0*GammaM/(_m0*_m0 - s)));
00826   double totalR = deltaR + _phiR ;
00827 
00828   //sum them up
00829   EvtComplex  bkgB,resT;
00830   bkgB = EvtComplex(_Blass*sin(totalB),0)*EvtComplex(cos(totalB),sin(totalB));
00831   resT = EvtComplex(_R*sin(deltaR),0)*EvtComplex(cos(totalR),sin(totalR))*EvtComplex(cos(2*totalB),sin(2*totalB));
00832   EvtComplex T = bkgB + resT;      
00833 
00834   return T;
00835 }

EvtComplex EvtDalitzReso::mixFactor	(	EvtComplex	prop,
		EvtComplex	prop_mix
	)			`[private]`

Definition at line 483 of file EvtDalitzReso.cc.

References _amp_mix, _delta_mix, _m0, and _m0_mix.

Referenced by evaluate().

00484 {
00485   double Delta = _delta_mix*(_m0+_m0_mix);
00486   return 1/(1-Delta*Delta*prop*prop_mix)*(1+_amp_mix*Delta*prop_mix);
00487 }

EvtComplex EvtDalitzReso::numerator	(	const EvtDalitzPoint &	p,
		const EvtTwoBodyKine &	vb,
		const EvtTwoBodyKine &	vd
	)			`[private]`

Definition at line 398 of file EvtDalitzReso.cc.

References _m0, _pairAng, _spin, _typeN, _vb, _vd, EvtTwoBodyKine::AB, angDep(), EvtCyclic3::combine(), EvtCyclic3::common(), GAUSS_CLEO, GAUSS_CLEO_ZEMACH, GS_CLEO, GS_CLEO_ZEMACH, NBW, EvtCyclic3::other(), EvtTwoBodyVertex::phaseSpaceFactor(), RBW_CLEO, RBW_CLEO_ZEMACH, RBW_KUEHN, RBW_ZEMACH, RBW_ZEMACH2, EvtSpinType::SCALAR, EvtSpinType::TENSOR, EvtSpinType::VECTOR, and x.

Referenced by evaluate().

00399 {
00400   EvtComplex ret(0.,0.);
00401 
00402   // Non-relativistic Breit-Wigner
00403   if(NBW == _typeN) {
00404     ret = angDep(x);
00405   }
00406 
00407   // Standard relativistic Zemach propagator
00408   else if(RBW_ZEMACH == _typeN) {
00409     ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*angDep(x);
00410   }
00411 
00412   // Standard relativistic Zemach propagator
00413   else if(RBW_ZEMACH2 == _typeN) {
00414     ret = _vd.phaseSpaceFactor(vd,EvtTwoBodyKine::AB)*_vb.phaseSpaceFactor(vb,EvtTwoBodyKine::AB)*angDep(x);
00415     if(_spin == EvtSpinType::VECTOR) {
00416       ret *= -4.;
00417     } else if(_spin == EvtSpinType::TENSOR) {
00418       ret *= 16./3.;
00419     } else if(_spin != EvtSpinType::SCALAR)
00420       assert(0);
00421   }
00422 
00423   // Kuehn-Santamaria normalization:
00424   else if(RBW_KUEHN == _typeN) {
00425     ret = _m0*_m0 * angDep(x);
00426   }  
00427 
00428   // CLEO amplitude 
00429   else if( ( RBW_CLEO        == _typeN ) || ( GS_CLEO           == _typeN ) ||
00430            ( RBW_CLEO_ZEMACH == _typeN ) || ( GS_CLEO_ZEMACH    == _typeN ) ||
00431            ( GAUSS_CLEO      == _typeN ) || ( GAUSS_CLEO_ZEMACH == _typeN)) {
00432 
00433     Index iA = other(_pairAng);           // A = other(BC)
00434     Index iB = common(_pairRes,_pairAng); // B = common(AB,BC)
00435     Index iC = other(_pairRes);           // C = other(AB)
00436     
00437     double M = x.bigM();
00438     double mA = x.m(iA);
00439     double mB = x.m(iB);
00440     double mC = x.m(iC);
00441     double qAB = x.q(combine(iA,iB));
00442     double qBC = x.q(combine(iB,iC));
00443     double qCA = x.q(combine(iC,iA));
00444 
00445     double M2 = M*M;
00446     double m02 = ((RBW_CLEO_ZEMACH == _typeN)||(GS_CLEO_ZEMACH == _typeN)||(GAUSS_CLEO_ZEMACH == _typeN))?  qAB : _m0*_m0;
00447     double mA2 = mA*mA;
00448     double mB2 = mB*mB;
00449     double mC2 = mC*mC;
00450     
00451     if (_spin == EvtSpinType::SCALAR) ret = EvtComplex(1.,0.);
00452     else if(_spin == EvtSpinType::VECTOR) {
00453       ret = qCA - qBC + (M2 - mC2)*(mB2 - mA2)/m02;
00454     } else if(_spin == EvtSpinType::TENSOR) {
00455       double x1 = qBC - qCA + (M2 - mC2)*(mA2 - mB2)/m02;       
00456       double x2 = M2 - mC2;      
00457       double x3 = qAB - 2*M2 - 2*mC2 + x2*x2/m02;      
00458       double x4 = mA2 - mB2;
00459       double x5 = qAB - 2*mB2 - 2*mA2 + x4*x4/m02;
00460       ret = x1*x1 - x3*x5/3.;
00461     } else assert(0);
00462   }
00463   
00464   return ret;
00465 }

EvtComplex EvtDalitzReso::propBreitWigner	(	const double &	m0,
		const double &	g0,
		const double &	m
	)			`[private]`

Definition at line 341 of file EvtDalitzReso.cc.

References EvtConst::twoPi.

Referenced by evaluate().

00342 {
00343   // non-relativistic BW
00344   return sqrt(g0/EvtConst::twoPi)/(m-m0-EvtComplex(0.0,g0/2.));
00345 }

EvtComplex EvtDalitzReso::propBreitWignerRel	(	const double &	m0,
		const EvtComplex &	g0,
		const double &	m
	)			`[private]`

Definition at line 356 of file EvtDalitzReso.cc.

00357 {
00358   // relativistic BW with complex width
00359   return 1./(m0*m0-m*m-EvtComplex(0.,m0)*g0);
00360 }

EvtComplex EvtDalitzReso::propBreitWignerRel	(	const double &	m0,
		const double &	g0,
		const double &	m
	)			`[private]`

Definition at line 348 of file EvtDalitzReso.cc.

Referenced by evaluate().

00349 {
00350   // relativistic BW with real width
00351   return 1./(m0*m0-m*m-EvtComplex(0.,m0*g0));
00352 }

EvtComplex EvtDalitzReso::propBreitWignerRelCoupled	(	const double &	m0,
		const EvtComplex &	g1,
		const EvtComplex &	g2,
		const double &	m
	)			`[private]`

Definition at line 363 of file EvtDalitzReso.cc.

Referenced by evaluate().

00364 {
00365   // relativistic coupled BW
00366   return 1./(m0*m0-m*m-(g1+g2));
00367 }

EvtComplex EvtDalitzReso::propGauss	(	const double &	m0,
		const double &	s0,
		const double &	m
	)			`[private]`

Definition at line 333 of file EvtDalitzReso.cc.

References exp(), and EvtConst::twoPi.

Referenced by evaluate().

00334 {
00335   // Gaussian
00336   double gauss = 1./sqrt(EvtConst::twoPi)/s0*exp(-(m-m0)*(m-m0)/2./(s0*s0));
00337   return EvtComplex(gauss,0.);
00338 }

EvtComplex EvtDalitzReso::propGounarisSakurai	(	const double &	m0,
		const double &	g0,
		const double &	k0,
		const double &	m,
		const double &	g,
		const double &	k
	)			`[private]`

Definition at line 369 of file EvtDalitzReso.cc.

References GS_d(), and GS_f().

Referenced by evaluate().

00371 {
00372   // Gounaris-Sakurai parameterization of pi+pi- P wave. PRD, Vol61, 112002. PRL, Vol21, 244.
00373   // Expressions taken from BAD637v4, after fixing the imaginary part of the BW denominator: i M_R Gamma_R(s) --> i sqrt(s) Gamma_R(s) 
00374   return (1.+GS_d(m0,k0)*g0/m0)/(m0*m0-m*m-EvtComplex(0.,m*g)+GS_f(m0,g0,k0,m,k));
00375 }

EvtComplex EvtDalitzReso::psFactor	(	double &	ma1,
		double &	mb1,
		double &	ma2,
		double &	mb2,
		double &	m
	)			`[private]`

Definition at line 327 of file EvtDalitzReso.cc.

References psFactor().

00328 {
00329   return 0.5*(psFactor(ma1,mb1,m)+psFactor(ma2,mb2,m));
00330 }

EvtComplex EvtDalitzReso::psFactor	(	double &	ma,
		double &	mb,
		double &	m
	)			`[private]`

Definition at line 313 of file EvtDalitzReso.cc.

References EvtTwoBodyKine::p(), and s.

Referenced by evaluate(), and psFactor().

00314 {
00315   if (m>(ma+mb)) {
00316     EvtTwoBodyKine vd(ma,mb,m);
00317     return EvtComplex(0,2*vd.p()/m);
00318   } else { 
00319     // analytical continuation
00320     double s = m*m;
00321     double phaseFactor_analyticalCont = -0.5*(sqrt(4*ma*ma/s-1)+sqrt(4*mb*mb/s-1)); 
00322     return EvtComplex(phaseFactor_analyticalCont,0);
00323   }
00324 }

void EvtDalitzReso::set_fb ( double R ) [inline]

Definition at line 82 of file EvtDalitzReso.hh.

References _vb, and EvtTwoBodyVertex::set_f().

00082 { _vb.set_f( R ); }

void EvtDalitzReso::set_fd ( double R ) [inline]

Definition at line 81 of file EvtDalitzReso.hh.

References _vd, and EvtTwoBodyVertex::set_f().

00081 { _vd.set_f( R ); }

double EvtDalitzReso::_a [private]

Definition at line 142 of file EvtDalitzReso.hh.

Referenced by lass().

EvtComplex EvtDalitzReso::_amp_mix [private]

Definition at line 130 of file EvtDalitzReso.hh.

Referenced by mixFactor().

double EvtDalitzReso::_Blass [private]

Definition at line 144 of file EvtDalitzReso.hh.

Referenced by lass().

CouplingType EvtDalitzReso::_coupling2 [private]

Definition at line 134 of file EvtDalitzReso.hh.

Referenced by evaluate().

double EvtDalitzReso::_delta_mix [private]

Definition at line 129 of file EvtDalitzReso.hh.

Referenced by mixFactor().

EvtDalitzPlot EvtDalitzReso::_dp [private]

Definition at line 106 of file EvtDalitzReso.hh.

Referenced by evaluate().

EvtComplex EvtDalitzReso::_fr12prod [private]

Definition at line 138 of file EvtDalitzReso.hh.

Referenced by Fvector().

EvtComplex EvtDalitzReso::_fr13prod [private]

Definition at line 138 of file EvtDalitzReso.hh.

Referenced by Fvector().

EvtComplex EvtDalitzReso::_fr14prod [private]

Definition at line 138 of file EvtDalitzReso.hh.

Referenced by Fvector().

EvtComplex EvtDalitzReso::_fr15prod [private]

Definition at line 138 of file EvtDalitzReso.hh.

Referenced by Fvector().

double EvtDalitzReso::_g0 [private]

Definition at line 119 of file EvtDalitzReso.hh.

Referenced by evaluate(), and lass().

double EvtDalitzReso::_g0_mix [private]

Definition at line 129 of file EvtDalitzReso.hh.

Referenced by evaluate().

double EvtDalitzReso::_g1 [private]

Definition at line 133 of file EvtDalitzReso.hh.

Referenced by evaluate().

double EvtDalitzReso::_g2 [private]

Definition at line 133 of file EvtDalitzReso.hh.

Referenced by evaluate().

int EvtDalitzReso::_kmatrix_index [private]

Definition at line 137 of file EvtDalitzReso.hh.

Referenced by evaluate().

double EvtDalitzReso::_m0 [private]

Definition at line 119 of file EvtDalitzReso.hh.

Referenced by evaluate(), lass(), mixFactor(), and numerator().

double EvtDalitzReso::_m0_mix [private]

Definition at line 129 of file EvtDalitzReso.hh.

Referenced by evaluate(), and mixFactor().

double EvtDalitzReso::_massFirst [private]

Definition at line 126 of file EvtDalitzReso.hh.

Referenced by evaluate(), GS_d(), GS_h(), and lass().

double EvtDalitzReso::_massSecond [private]

Definition at line 126 of file EvtDalitzReso.hh.

Referenced by evaluate(), and lass().

EvtCyclic3::Pair EvtDalitzReso::_pairAng [private]

Definition at line 109 of file EvtDalitzReso.hh.

Referenced by angDep(), and numerator().

EvtCyclic3::Pair EvtDalitzReso::_pairRes [private]

Definition at line 110 of file EvtDalitzReso.hh.

double EvtDalitzReso::_phiB [private]

Definition at line 145 of file EvtDalitzReso.hh.

Referenced by lass().

double EvtDalitzReso::_phiR [private]

Definition at line 147 of file EvtDalitzReso.hh.

Referenced by lass().

double EvtDalitzReso::_R [private]

Definition at line 146 of file EvtDalitzReso.hh.

Referenced by lass().

double EvtDalitzReso::_r [private]

Definition at line 143 of file EvtDalitzReso.hh.

Referenced by lass().

double EvtDalitzReso::_s0prod [private]

Definition at line 139 of file EvtDalitzReso.hh.

Referenced by Fvector().

EvtSpinType::spintype EvtDalitzReso::_spin [private]

Definition at line 113 of file EvtDalitzReso.hh.

Referenced by angDep(), evaluate(), and numerator().

NumType EvtDalitzReso::_typeN [private]

Definition at line 116 of file EvtDalitzReso.hh.

Referenced by evaluate(), Fvector(), and numerator().

EvtTwoBodyVertex EvtDalitzReso::_vb [private]

Definition at line 122 of file EvtDalitzReso.hh.

Referenced by evaluate(), numerator(), and set_fb().

EvtTwoBodyVertex EvtDalitzReso::_vd [private]

Definition at line 123 of file EvtDalitzReso.hh.

Referenced by evaluate(), lass(), numerator(), and set_fd().

EvtDalitzReso Class Reference

Public Types

Public Member Functions

Private Member Functions

Private Attributes

Detailed Description

Member Enumeration Documentation

Constructor & Destructor Documentation

Member Function Documentation

Member Data Documentation


Public Types
	NBW = 0
	RBW_ZEMACH = 1
	RBW_KUEHN = 2
	RBW_CLEO = 3
	RBW_ZEMACH2 = 4
	GS_CLEO = 5
	K_MATRIX = 6
	RBW_CLEO_ZEMACH = 7
	GS_CLEO_ZEMACH = 8
	LASS = 9
	K_MATRIX_I = 10
	K_MATRIX_II = 11
	GAUSS_CLEO = 12
	GAUSS_CLEO_ZEMACH = 13
	Undefined = 0
	PicPic = 1
	PizPiz
	PiPi
	KcKc
	KzKz
	KK
	EtaPic
	EtaPiz
	PicPicKK
	WA76
enum	NumType { NBW = 0, RBW_ZEMACH = 1, RBW_KUEHN = 2, RBW_CLEO = 3, RBW_ZEMACH2 = 4, GS_CLEO = 5, K_MATRIX = 6, RBW_CLEO_ZEMACH = 7, GS_CLEO_ZEMACH = 8, LASS = 9, K_MATRIX_I = 10, K_MATRIX_II = 11, GAUSS_CLEO = 12, GAUSS_CLEO_ZEMACH = 13 }
enum	CouplingType { Undefined = 0, PicPic = 1, PizPiz, PiPi, KcKc, KzKz, KK, EtaPic, EtaPiz, PicPicKK, WA76 }
Public Member Functions
	EvtDalitzReso ()
	EvtDalitzReso (const EvtDalitzPlot &dp, EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes, EvtSpinType::spintype spin, double m0, double g0, NumType typeN)
	EvtDalitzReso (const EvtDalitzPlot &dp, EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes, EvtSpinType::spintype spin, double m0, double g0, NumType typeN, double m0_mix, double g0_mix, double delta_mix, EvtComplex amp_mix)
	EvtDalitzReso (const EvtDalitzPlot &dp, EvtCyclic3::Pair pairAng, EvtCyclic3::Pair pairRes, EvtSpinType::spintype spin, double m0, NumType typeN, double g1, double g2, CouplingType coupling2)
	EvtDalitzReso (const EvtDalitzPlot &dp, EvtCyclic3::Pair pairRes, std::string nameIndex, NumType typeN, EvtComplex fr12prod, EvtComplex fr13prod, EvtComplex fr14prod, EvtComplex fr15prod, double s0prod)
	EvtDalitzReso (const EvtDalitzPlot &dp, EvtCyclic3::Pair pairRes, double m0, double g0, double a, double r, double B, double phiB, double R, double phiR)
	EvtDalitzReso (const EvtDalitzReso &other)
	~EvtDalitzReso ()
EvtDalitzReso *	clone () const
EvtComplex	evaluate (const EvtDalitzPoint &p)
void	set_fd (double R)
void	set_fb (double R)
Private Member Functions
EvtComplex	psFactor (double &ma, double &mb, double &m)
EvtComplex	psFactor (double &ma1, double &mb1, double &ma2, double &mb2, double &m)
EvtComplex	propGauss (const double &m0, const double &s0, const double &m)
EvtComplex	propBreitWigner (const double &m0, const double &g0, const double &m)
EvtComplex	propBreitWignerRel (const double &m0, const double &g0, const double &m)
EvtComplex	propBreitWignerRel (const double &m0, const EvtComplex &g0, const double &m)
EvtComplex	propBreitWignerRelCoupled (const double &m0, const EvtComplex &g1, const EvtComplex &g2, const double &m)
EvtComplex	propGounarisSakurai (const double &m0, const double &g0, const double &k0, const double &m, const double &g, const double &k)
double	GS_f (const double &m0, const double &g0, const double &k0, const double &m, const double &k)
double	GS_h (const double &m, const double &k)
double	GS_dhods (const double &m0, const double &k0)
double	GS_d (const double &m0, const double &k0)
EvtComplex	numerator (const EvtDalitzPoint &p, const EvtTwoBodyKine &vb, const EvtTwoBodyKine &vd)
double	angDep (const EvtDalitzPoint &p)
EvtComplex	mixFactor (EvtComplex prop, EvtComplex prop_mix)
EvtComplex	Fvector (double s, int index)
EvtComplex	lass (double s)
Private Attributes
EvtDalitzPlot	_dp
EvtCyclic3::Pair	_pairAng
EvtCyclic3::Pair	_pairRes
EvtSpinType::spintype	_spin
NumType	_typeN
double	_m0
double	_g0
EvtTwoBodyVertex	_vb
EvtTwoBodyVertex	_vd
double	_massFirst
double	_massSecond
double	_m0_mix
double	_g0_mix
double	_delta_mix
EvtComplex	_amp_mix
double	_g1
double	_g2
CouplingType	_coupling2
int	_kmatrix_index
EvtComplex	_fr12prod
EvtComplex	_fr13prod
EvtComplex	_fr14prod
EvtComplex	_fr15prod
double	_s0prod
double	_a
double	_r
double	_Blass
double	_phiB
double	_R
double	_phiR