TEvtGen/EvtGenModels/EvtSVVNONCPEIGEN.cpp

   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) 2001      Royal Holloway, University of London
  10 //
  11 // Module: EvtSVVNONCPEIGEN.cc
  12 //
  13 // Description: Routine to decay scalar -> vector vector
  14 //              and has CP violation.
  15 //
  16 //              This model does all the ckm-suppressed decays and mixing for you. It randomly 'overwrites'
  17 //              any reco or tagging state as set in the Y(4S) decay model (VSS_(B)MIX) with its own generated states.
  18 //
  19 //              As such, the corresponding dec file requires only one decay-mode description, for example:
  20 //              Decay MyB0
  21 //              1.000    rho+ MyD*-       SVV_NONCPEIGEN dm beta gamma 0.322 0.31 0.941 0 0.107 1.42 0.02 0 0.02 0 0.02 0 ;
  22 //              EndDecay
  23 //              and furthermore Y(4S) only needs to decay to B0's (or B0bar's).
  24 //              The decay above should be a CKM-favored mode (eg. B0->D*-rho+ or B0bar->D*+rho-).
  25 //              All ckm-suppressed decays and the mixing are derived from this line in the ::Decay function.
  26 //
  27 //              There are 15 or 27 arguments. The first three are dm, phase1
  28 //              and phase2. dm is the B0-B0bar mass difference. Phases 1
  29 //              and 2 are the CKM weak phases relevant for the particular mode,
  30 //              eg for B-->DstRho phase1 is beta and phase2 is gamma.
  31 //
  32 //              The next arguments are the 2 amplitudes (= 12 input parameters)
  33 //              in the order: A_f, Abar_f. In the example above, the 'A_f' amplitude now
  34 //              stands for the ckm-favored decay 'B0->D*-rho+', and 'Abar_f' stands for 'B0bar->D*-rho+'
  35 //
  36 //              Each amplitude has its 3 helicity states in the order +, 0, -, which are each
  37 //              specified by a magnitude and a strong phase.
  38 //
  39 //              The last 2 arguments A_fbar and Abar_fbar (=12 input parameters) are not necessary,
  40 //              but can included if one wants to set them differently from A_f, Abar_f.
  41 //
  42 //              Mind you that Hbar_+- = H_-+ (ignoring the weak phase, which flips sign).
  43 //              It is custumary to select one set of helicity states (eg H_+-) and to adopt these for
  44 //              the CP-conjugate decays as well (ie. depict Hbar_-+ with H_+-), which is the interpretation
  45 //              we use for the input-parameters above.
  46 //              However, the angular decay in EvtGen is just a formula in which helicity amplitudes are 'plugged' in,
  47 //              making no difference between B0 or B0bar decays. In the model below we (thus) account for the +-
  48 //              flipping between B0 and B0bar.
  49 //
  50 //
  51 // Modification history:
  52 //    Ajit Kurup 9 March 2001        Module created (from EvtSVSNONCPEIGEN)
  53 //    Max Baak 01/16/2004            Fix of Helicity amplitude ordering.
  54 //                                   Decay also works for B0bar decays.
  55 //------------------------------------------------------------------------
  56 //
  57 #include "EvtGenBase/EvtPatches.hh"
  58 #include <stdlib.h>
  59 #include "EvtGenBase/EvtParticle.hh"
  60 #include "EvtGenBase/EvtRandom.hh"
  61 #include "EvtGenBase/EvtGenKine.hh"
  62 #include "EvtGenBase/EvtCPUtil.hh"
  63 #include "EvtGenBase/EvtPDL.hh"
  64 #include "EvtGenBase/EvtReport.hh"
  65 #include "EvtGenBase/EvtVector4C.hh"
  66 #include "EvtGenModels/EvtSVVNONCPEIGEN.hh"
  67 #include <string>
  68 #include "EvtGenModels/EvtSVVHelAmp.hh"
  69 #include "EvtGenBase/EvtConst.hh"
  70
  71 EvtSVVNONCPEIGEN::~EvtSVVNONCPEIGEN() {}
  72
  73 std::string EvtSVVNONCPEIGEN::getName(){
  74
  75   return "SVV_NONCPEIGEN";
  76
  77 }
  78
  79
  80 EvtDecayBase* EvtSVVNONCPEIGEN::clone(){
  81
  82   return new EvtSVVNONCPEIGEN;
  83
  84 }
  85
  86 void EvtSVVNONCPEIGEN::init(){
  87
  88   // check that there are 27 arguments
  89   checkNArg(27,15);
  90   checkNDaug(2);
  91
  92   checkSpinDaughter(0,EvtSpinType::VECTOR);
  93   checkSpinDaughter(1,EvtSpinType::VECTOR);
  94
  95   //  The ordering of A_f is :
  96   //  A_f[0-2] = A_f
  97   //  A_f[3-5] = Abar_f
  98   //  A_f[6-8] = A_fbar
  99   //  A_f[9-11] = Abar_fbar
 100   //
 101   //  Each of the 4 amplitudes include the 3 different helicity states in
 102   //  the order +, 0, -. See more about helicity amplitude ordering in ::decay
 103
 104   int i=0;
 105   int j=(getNArg()-3)/2;
 106
 107   for(i=0; i<j; ++i){
 108     _A_f[i] = getArg((2*i)+3) * EvtComplex( cos(getArg((2*i)+4)),sin(getArg((2*i)+4)) );
 109   }
 110
 111   //  If only 6 amplitudes are specified, calculate the last 6 from the first 6:
 112   if(6 == j){
 113     for(i = 0; i < 3; ++i){
 114       _A_f[6+i] = _A_f[3+i];
 115       _A_f[9+i] = _A_f[i];
 116     }
 117   }
 118 }
 119
 120 void EvtSVVNONCPEIGEN::initProbMax() {
 121   double probMax = 0;
 122   for (int i = 0; i < 12; ++i){
 123     double amp = abs(_A_f[i]);
 124     probMax += amp * amp;
 125   }
 126
 127   setProbMax(probMax);
 128 }
 129
 130 void EvtSVVNONCPEIGEN::decay( EvtParticle *p){
 131
 132   //added by Lange Jan4,2000
 133   static EvtId B0=EvtPDL::getId("B0");
 134   static EvtId B0B=EvtPDL::getId("anti-B0");
 135
 136   double t;
 137   EvtId other_b;
 138   EvtId daugs[2];
 139
 140
 141   // MB: flip selects the final of the decay
 142   int flip = ((p->getId() == B0) ? 0 : 1);
 143   daugs[0]=getDaug(0);
 144   daugs[1]=getDaug(1);
 145   p->initializePhaseSpace(2,daugs);
 146
 147   EvtCPUtil::getInstance()->OtherB(p,t,other_b,0.5);
 148
 149   EvtComplex amp[3];
 150
 151   double dmt2 = getArg(0) * t / (2 * EvtConst::c);
 152   double phiCKM = (2.0 * getArg(1) + getArg(2));   // 2b+g
 153   EvtComplex ePlusIPhi(cos(phiCKM), sin(phiCKM));
 154   EvtComplex eMinusIPhi(cos(-phiCKM), sin(-phiCKM));
 155
 156   // flip == 0 : D*-rho+
 157   // flip == 1 : D*+rho-
 158
 159   if (!flip) {
 160     if (other_b==B0B){
 161       // At t=0 we have a B0
 162       for (int i=0; i<3; ++i) {
 163         amp[i] = _A_f[i]*cos(dmt2) + eMinusIPhi*EvtComplex(0.0,sin(dmt2))*_A_f[i+3];
 164       }
 165     }
 166     if (other_b==B0){
 167       // At t=0 we have a B0bar
 168       for(int i=0; i<3; ++i) {
 169         amp[i] = _A_f[i]*ePlusIPhi*EvtComplex(0.0,sin(dmt2)) + _A_f[i+3]*cos(dmt2);
 170       }
 171     }
 172   } else{
 173     if (other_b==B0B){
 174       // At t=0 we have a B0
 175
 176       // M.Baak 01/16/2004
 177       // Note: \bar{H}+- = H-+
 178       // If one wants to use the correct helicities for B0 and B0bar decays but the same formula-notation (as done in EvtSVV_HelAmp),
 179       // count the B0bar helicities backwards. (Equivalently, one could flip the chi angle.)
 180
 181       for(int i=0; i<3; ++i) {
 182         amp[i] = _A_f[8-i]*cos(dmt2) + eMinusIPhi*EvtComplex(0.0,sin(dmt2))*_A_f[11-i];
 183       }
 184     }
 185     if (other_b==B0){
 186       // At t=0 we have a B0bar
 187       for(int i=0; i<3; ++i) {
 188         amp[i] = _A_f[8-i] * ePlusIPhi * EvtComplex(0.0,sin(dmt2)) + _A_f[11-i]*cos(dmt2);
 189       }
 190     }
 191   }
 192
 193   EvtSVVHelAmp::SVVHel(p,_amp2,daugs[0],daugs[1],amp[0],amp[1],amp[2]);
 194
 195   return ;
 196 }
 197
 198 std::string EvtSVVNONCPEIGEN::getParamName(int i) {
 199   switch(i) {
 200   case 0:
 201     return "deltaM";
 202   case 1:
 203     return "weakPhase1";
 204   case 2:
 205     return "weakPhase2";
 206   case 3:
 207     return "AfPlusHelAmp";
 208   case 4:
 209     return "AfPlusHelAmpPhase";
 210   case 5:
 211     return "AfZeroHelAmp";
 212   case 6:
 213     return "AfZeroHelAmpPhase";
 214   case 7:
 215     return "AfMinusHelAmp";
 216   case 8:
 217     return "AfMinusHelAmpPhase";
 218   case 9:
 219     return "AbarfPlusHelAmp";
 220   case 10:
 221     return "AbarfPlusHelAmpPhase";
 222   case 11:
 223     return "AbarfZeroHelAmp";
 224   case 12:
 225     return "AbarfZeroHelAmpPhase";
 226   case 13:
 227     return "AbarfMinusHelAmp";
 228   case 14:
 229     return "AbarfMinusHelAmpPhase";
 230   case 15:
 231     return "AfbarPlusHelAmp";
 232   case 16:
 233     return "AfbarPlusHelAmpPhase";
 234   case 17:
 235     return "AfbarZeroHelAmp";
 236   case 18:
 237     return "AfbarZeroHelAmpPhase";
 238   case 19:
 239     return "AfbarMinusHelAmp";
 240   case 20:
 241     return "AfbarMinusHelAmpPhase";
 242   case 21:
 243     return "AbarfbarPlusHelAmp";
 244   case 22:
 245     return "AbarfbarPlusHelAmpPhase";
 246   case 23:
 247     return "AbarfbarZeroHelAmp";
 248   case 24:
 249     return "AbarfbarZeroHelAmpPhase";
 250   case 25:
 251     return "AbarfbarMinusHelAmp";
 252   case 26:
 253     return "AbarfbarMinusHelAmpPhase";
 254   default:
 255     return "";
 256   }
 257 }
 258
 259 std::string EvtSVVNONCPEIGEN::getParamDefault(int i) {
 260   switch(i) {
 261   case 3:
 262     return "1.0";
 263   case 4:
 264     return "0.0";
 265   case 5:
 266     return "1.0";
 267   case 6:
 268     return "0.0";
 269   case 7:
 270     return "1.0";
 271   case 8:
 272     return "0.0";
 273   case 9:
 274     return "1.0";
 275   case 10:
 276     return "0.0";
 277   case 11:
 278     return "1.0";
 279   case 12:
 280     return "0.0";
 281   case 13:
 282     return "1.0";
 283   case 14:
 284     return "0.0";
 285   default:
 286     return "";
 287   }
 288 }