Converting TEvtGen to native cmake

[u/mrichter/AliRoot.git] / TEvtGen / EvtGen / EvtGenModels / EvtWnPi.cpp

diff --git a/TEvtGen/EvtGen/EvtGenModels/EvtWnPi.cpp b/TEvtGen/EvtGen/EvtGenModels/EvtWnPi.cpp
new file mode 100644 (file)
index 0000000..b5a0403
--- /dev/null
+++ b/TEvtGen/EvtGen/EvtGenModels/EvtWnPi.cpp
@@ -0,0 +1,137 @@
+//--------------------------------------------------------------------------
+//
+// 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: EvtTVP.cc
+//
+// Description: Routine to calculate W -> (n pi) current
+//                     according to [Kuhn, Was, Acta.Phys.Polon B39 (2008) 147]
+//
+// Modification history:
+//     AVL     6 July, 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/EvtWnPi.hh"
+
+using namespace std;
+
+// W+ -> pi_ current
+EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1) { 
+  return q1;
+}
+
+// W+ -> pi+ pi0 current
+EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2) { 
+  return BWr(q1+q2)*(q1-q2);
+}
+
+// W+ -> pi+ pi+ pi- current
+EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3) {
+ EvtVector4R Q=q1+q2+q3;
+ double Q2=Q.mass2();
+ return BWa(Q)*( (q1-q3) - (Q*(Q*(q1-q3))/Q2)*BWr(q2+q3) +
+             (q2-q3) - (Q*(Q*(q2-q3))/Q2)*BWr(q1+q3) ); 
+}
+
+// W+ -> pi+ pi+ pi- pi- pi+ current with symmetrization
+EvtVector4C EvtWnPi::WCurrent(EvtVector4R q1, EvtVector4R q2, EvtVector4R q3, EvtVector4R q4, EvtVector4R q5) {
+//  double Q2 = Qtot*Qtot;
+//  return q1-Qtot*(q1*Qtot)/Q2;
+  EvtVector4C V = JB(q1, q2, q3, q4, q5) + JB(q5, q2, q3, q4, q1) + JB(q1, q5, q3, q4, q2) +
+                 JB(q1,q2,q4,q3,q5)+JB(q5,q2,q4,q3,q1)+JB(q1,q5,q4,q3,q2); 
+//  cout<<"BC2: Qtot="<<Qtot<<", V="<<V<<endl;
+  return V;
+}
+
+
+// a1 -> pi+ pi+ pi- BW
+EvtComplex EvtWnPi::BWa(EvtVector4R q) {
+  double const _mA1=1.26, _GA1=0.4;
+  EvtComplex I(0,1);
+  double Q2 = q.mass2();
+  double GA1=_GA1*pi3G(Q2)/pi3G(_mA1*_mA1);
+  EvtComplex denBA1(_mA1*_mA1 - Q2,-1.*_mA1*GA1);
+  return _mA1*_mA1 / denBA1; 
+}
+
+
+EvtComplex EvtWnPi::BWf(EvtVector4R q) {
+  double const mf=0.8, Gf=0.6;
+  EvtComplex I(0,1);
+  double Q2 = q.mass2();
+  return mf*mf/(mf*mf-Q2-I*mf*Gf); 
+}
+
+EvtComplex EvtWnPi::BWr(EvtVector4R q) {
+       double _mRho = 0.775, _gammaRho=0.149, _mRhopr=1.364,  _gammaRhopr=0.400, _beta=-0.108;
+       double m1=EvtPDL::getMeanMass(EvtPDL::getId("pi+")), m2=EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
+       double mQ2=q.mass2();
+       
+       // momenta in the rho->pipi decay
+       double dRho= _mRho*_mRho - m1*m1 - m2*m2;
+       double pPiRho = (1.0/_mRho)*sqrt((dRho*dRho)/4.0 - m1*m1*m2*m2);
+       
+       double dRhopr= _mRhopr*_mRhopr - m1*m1 - m2*m2;
+       double pPiRhopr = (1.0/_mRhopr)*sqrt((dRhopr*dRhopr)/4.0 - m1*m1*m2*m2);
+       
+       double dQ= mQ2 - m1*m1 - m2*m2;
+       double pPiQ = (1.0/sqrt(mQ2))*sqrt((dQ*dQ)/4.0 - m1*m1*m2*m2);
+       
+       
+       double gammaRho = _gammaRho*_mRho/sqrt(mQ2)*pow((pPiQ/pPiRho),3);
+       EvtComplex BRhoDem(_mRho*_mRho - mQ2,-1.0*_mRho*gammaRho);
+       EvtComplex BRho= _mRho*_mRho / BRhoDem;
+       
+       double gammaRhopr = _gammaRhopr*_mRhopr/sqrt(mQ2)*pow((pPiQ/pPiRhopr),3);
+       EvtComplex BRhoprDem(_mRhopr*_mRhopr - mQ2,-1.0*_mRho*gammaRhopr);
+       EvtComplex BRhopr= _mRhopr*_mRhopr / BRhoprDem;
+       
+       return (BRho + _beta*BRhopr)/(1+_beta);
+}
+
+double EvtWnPi::pi3G(double m2) {
+ double mPi = EvtPDL::getMeanMass(EvtPDL::getId("pi+"));
+ double _mRho = 0.775;
+  if ( m2 > (_mRho+mPi) ) {
+    return m2*(1.623 + 10.38/m2 - 9.32/(m2*m2) + 0.65/(m2*m2*m2));
+  }
+  else {
+    double t1=m2-9.0*mPi*mPi;
+    return 4.1*pow(t1,3.0)*(1.0 - 3.3*t1+5.8*t1*t1);
+  };
+}
+
+EvtVector4C EvtWnPi::JB( EvtVector4R p1, EvtVector4R p2, EvtVector4R p3, EvtVector4R p4, EvtVector4R p5) {
+  EvtVector4R Qtot = p1+p2+p3+p4+p5, Qa=p1+p2+p3;
+  EvtTensor4C T= (1/Qtot.mass2())*EvtGenFunctions::directProd(Qtot,Qtot) - EvtTensor4C::g();
+  EvtVector4R V13 = Qa*( p2*(p1-p3) )/Qa.mass2() - (p1-p3);
+  EvtVector4R V23 = Qa*( p1*(p2-p3) )/Qa.mass2() - (p2-p3);
+  return BWa(Qtot)*BWa(Qa)*BWf(p4+p5)*(
+    T.cont1(V13)*BWr(p1+p3) + T.cont1(V23)*BWr(p2+p3)
+  );
+}