TEvtGen/EvtGenBase/EvtSemiLeptonicAmp.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) 1998      Caltech, UCSB
  10 //
  11 // Module: EvtSemiLeptonicAmp.cc
  12 //
  13 // Description: Base class for semileptonic decays
  14 //
  15 // Modification history:
  16 //
  17 //    DJL       April 17,1998       Module created
  18 //
  19 //------------------------------------------------------------------------
  20 //
  21 #include "EvtGenBase/EvtPatches.hh"
  22 #include "EvtGenBase/EvtPatches.hh"
  23 #include "EvtGenBase/EvtParticle.hh"
  24 #include "EvtGenBase/EvtGenKine.hh"
  25 #include "EvtGenBase/EvtPDL.hh"
  26 #include "EvtGenBase/EvtReport.hh"
  27 #include "EvtGenBase/EvtVector4C.hh"
  28 #include "EvtGenBase/EvtTensor4C.hh"
  29 #include "EvtGenBase/EvtDiracSpinor.hh"
  30 #include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh"
  31 #include "EvtGenBase/EvtId.hh"
  32 #include "EvtGenBase/EvtAmp.hh"
  33 #include "EvtGenBase/EvtScalarParticle.hh"
  34 #include "EvtGenBase/EvtVectorParticle.hh"
  35 #include "EvtGenBase/EvtTensorParticle.hh"
  36
  37 double EvtSemiLeptonicAmp::CalcMaxProb( EvtId parent, EvtId meson,
  38                                         EvtId lepton, EvtId nudaug,
  39                      EvtSemiLeptonicFF *FormFactors ) {
  40
  41   //This routine takes the arguements parent, meson, and lepton
  42   //number, and a form factor model, and returns a maximum
  43   //probability for this semileptonic form factor model.  A
  44   //brute force method is used.  The 2D cos theta lepton and
  45   //q2 phase space is probed.
  46
  47   //Start by declaring a particle at rest.
  48
  49   //It only makes sense to have a scalar parent.  For now.
  50   //This should be generalized later.
  51
  52   EvtScalarParticle *scalar_part;
  53   EvtParticle *root_part;
  54
  55   scalar_part=new EvtScalarParticle;
  56
  57   //cludge to avoid generating random numbers!
  58   scalar_part->noLifeTime();
  59
  60   EvtVector4R p_init;
  61
  62   p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0);
  63   scalar_part->init(parent,p_init);
  64   root_part=(EvtParticle *)scalar_part;
  65   root_part->setDiagonalSpinDensity();
  66
  67   EvtParticle *daughter, *lep, *trino;
  68
  69   EvtAmp amp;
  70
  71   EvtId listdaug[3];
  72   listdaug[0] = meson;
  73   listdaug[1] = lepton;
  74   listdaug[2] = nudaug;
  75
  76   amp.init(parent,3,listdaug);
  77
  78   root_part->makeDaughters(3,listdaug);
  79   daughter=root_part->getDaug(0);
  80   lep=root_part->getDaug(1);
  81   trino=root_part->getDaug(2);
  82
  83   //cludge to avoid generating random numbers!
  84   daughter->noLifeTime();
  85   lep->noLifeTime();
  86   trino->noLifeTime();
  87
  88
  89   //Initial particle is unpolarized, well it is a scalar so it is
  90   //trivial
  91   EvtSpinDensity rho;
  92   rho.setDiag(root_part->getSpinStates());
  93
  94   double mass[3];
  95
  96   double m = root_part->mass();
  97
  98   EvtVector4R p4meson, p4lepton, p4nu, p4w;
  99   double q2max;
 100
 101   double q2, elepton, plepton;
 102   int i,j;
 103   double erho,prho,costl;
 104
 105   double maxfoundprob = 0.0;
 106   double prob = -10.0;
 107   int massiter;
 108
 109   for (massiter=0;massiter<3;massiter++){
 110
 111     mass[0] = EvtPDL::getMeanMass(meson);
 112     mass[1] = EvtPDL::getMeanMass(lepton);
 113     mass[2] = EvtPDL::getMeanMass(nudaug);
 114     if ( massiter==1 ) {
 115       mass[0] = EvtPDL::getMinMass(meson);
 116     }
 117     if ( massiter==2 ) {
 118       mass[0] = EvtPDL::getMaxMass(meson);
 119       if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001;
 120     }
 121
 122     q2max = (m-mass[0])*(m-mass[0]);
 123
 124     //loop over q2
 125
 126     for (i=0;i<25;i++) {
 127       q2 = ((i+0.5)*q2max)/25.0;
 128
 129       erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m);
 130
 131       prho = sqrt(erho*erho-mass[0]*mass[0]);
 132
 133       p4meson.set(erho,0.0,0.0,-1.0*prho);
 134       p4w.set(m-erho,0.0,0.0,prho);
 135
 136       //This is in the W rest frame
 137       elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2));
 138       plepton = sqrt(elepton*elepton-mass[1]*mass[1]);
 139
 140       double probctl[3];
 141
 142       for (j=0;j<3;j++) {
 143
 144         costl = 0.99*(j - 1.0);
 145
 146         //These are in the W rest frame. Need to boost out into
 147         //the B frame.
 148         p4lepton.set(elepton,0.0,
 149                   plepton*sqrt(1.0-costl*costl),plepton*costl);
 150         p4nu.set(plepton,0.0,
 151                  -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl);
 152
 153         EvtVector4R boost((m-erho),0.0,0.0,1.0*prho);
 154         p4lepton=boostTo(p4lepton,boost);
 155         p4nu=boostTo(p4nu,boost);
 156
 157         //Now initialize the daughters...
 158
 159         daughter->init(meson,p4meson);
 160         lep->init(lepton,p4lepton);
 161         trino->init(nudaug,p4nu);
 162
 163         CalcAmp(root_part,amp,FormFactors);
 164
 165         //Now find the probability at this q2 and cos theta lepton point
 166         //and compare to maxfoundprob.
 167
 168         //Do a little magic to get the probability!!
 169         prob = rho.normalizedProb(amp.getSpinDensity());
 170
 171         probctl[j]=prob;
 172       }
 173
 174       //probclt contains prob at ctl=-1,0,1.
 175       //prob=a+b*ctl+c*ctl^2
 176
 177       double a=probctl[1];
 178       double b=0.5*(probctl[2]-probctl[0]);
 179       double c=0.5*(probctl[2]+probctl[0])-probctl[1];
 180
 181       prob=probctl[0];
 182       if (probctl[1]>prob) prob=probctl[1];
 183       if (probctl[2]>prob) prob=probctl[2];
 184
 185       if (fabs(c)>1e-20){
 186         double ctlx=-0.5*b/c;
 187         if (fabs(ctlx)<1.0){
 188           double probtmp=a+b*ctlx+c*ctlx*ctlx;
 189           if (probtmp>prob) prob=probtmp;
 190         }
 191
 192       }
 193
 194       if ( prob > maxfoundprob ) {
 195         maxfoundprob = prob;
 196       }
 197
 198     }
 199     if ( EvtPDL::getWidth(meson) <= 0.0 ) {
 200       //if the particle is narrow dont bother with changing the mass.
 201       massiter = 4;
 202     }
 203
 204   }
 205   root_part->deleteTree();
 206
 207   maxfoundprob *=1.1;
 208   return maxfoundprob;
 209
 210 }
 211