TEvtGen/EvtGenModels/EvtYmSToYnSpipiCLEO.cxx

   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 + 2*Re(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    = (mX*mX - 2.0 * m_pi * m_pi);
 209
 210       double deltaEmax =
 211         - 2.0 *
 212         sqrt( P_nS.get(0)*P_nS.get(0) - M_nS*M_nS ) *
 213         sqrt( 0.25 - pow(m_pi/mX,2.0));
 214
 215       double sumE = (M_mS*M_mS - M_nS*M_nS + mX*mX)/(2.0 * M_mS);
 216
 217       double E1E2 = 0.25 * ( pow(sumE, 2.0) - pow( deltaEmax * costheta, 2.0) );
 218
 219       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;
 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_mS*M_mS - pow(M_nS + mX,2.0)) * (M_mS*M_mS - pow(M_nS - mX,2.0)) ) * // p(*)_X
 231         sqrt(mX*mX - 4*m_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