#include <EvtSemiLeptonicScalarAmp.hh>
Inheritance diagram for EvtSemiLeptonicScalarAmp:
Public Member Functions | |
double | CalcMaxProb (EvtId parent, EvtId meson, EvtId lepton, EvtId nudaug, EvtSemiLeptonicFF *FormFactors) |
Private Member Functions | |
void | CalcAmp (EvtParticle *parent, EvtAmp &, EvtSemiLeptonicFF *FormFactors) |
Definition at line 30 of file EvtSemiLeptonicScalarAmp.hh.
void EvtSemiLeptonicScalarAmp::CalcAmp | ( | EvtParticle * | parent, | |
EvtAmp & | amp, | |||
EvtSemiLeptonicFF * | FormFactors | |||
) | [private, virtual] |
Reimplemented from EvtSemiLeptonicAmp.
Definition at line 36 of file EvtSemiLeptonicScalarAmp.cc.
References calibUtil::ERROR, EvtLeptonVACurrent(), EvtParticle::getDaug(), EvtParticle::getId(), EvtPDL::getId(), EvtParticle::getP4(), EvtSemiLeptonicFF::getscalarff(), EvtParticle::mass(), q, report(), EvtVector4R::set(), EvtParticle::spParent(), EvtParticle::spParentNeutrino(), and EvtAmp::vertex().
00038 { 00039 00040 static EvtId EM=EvtPDL::getId("e-"); 00041 static EvtId MUM=EvtPDL::getId("mu-"); 00042 static EvtId TAUM=EvtPDL::getId("tau-"); 00043 static EvtId EP=EvtPDL::getId("e+"); 00044 static EvtId MUP=EvtPDL::getId("mu+"); 00045 static EvtId TAUP=EvtPDL::getId("tau+"); 00046 00047 //Add the lepton and neutrino 4 momenta to find q2 00048 00049 EvtVector4R q = parent->getDaug(1)->getP4() 00050 + parent->getDaug(2)->getP4(); 00051 double q2 = (q.mass2()); 00052 00053 double fpf,f0f; 00054 double mesonmass = parent->getDaug(0)->mass(); 00055 double parentmass = parent->mass(); 00056 00057 FormFactors->getscalarff(parent->getId(), 00058 parent->getDaug(0)->getId(), 00059 q2, 00060 mesonmass, 00061 &fpf, 00062 &f0f); 00063 00064 00065 EvtVector4R p4b; 00066 p4b.set(parent->mass(),0.0,0.0,0.0); 00067 EvtVector4R p4meson = parent->getDaug(0)->getP4(); 00068 double mdiffoverq2; 00069 mdiffoverq2 = parentmass*parentmass - mesonmass*mesonmass; 00070 mdiffoverq2 = mdiffoverq2 / q2; 00071 00072 EvtVector4C l1,l2; 00073 00074 EvtId l_num = parent->getDaug(1)->getId(); 00075 EvtVector4C tds; 00076 00077 if (l_num==EM||l_num==MUM||l_num==TAUM){ 00078 00079 tds = EvtVector4C(fpf*(p4b+p4meson - (mdiffoverq2*(p4b-p4meson)))+ 00080 + f0f*mdiffoverq2*(p4b-p4meson)); 00081 00082 l1=EvtLeptonVACurrent(parent->getDaug(1)->spParent(0), 00083 parent->getDaug(2)->spParentNeutrino()); 00084 l2=EvtLeptonVACurrent(parent->getDaug(1)->spParent(1), 00085 parent->getDaug(2)->spParentNeutrino()); 00086 } 00087 else{ 00088 if (l_num==EP||l_num==MUP||l_num==TAUP){ 00089 00090 tds = EvtVector4C(fpf*(p4b+p4meson - (mdiffoverq2*(p4b-p4meson)))+ 00091 + f0f*mdiffoverq2*(p4b-p4meson)); 00092 00093 l1=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(), 00094 parent->getDaug(1)->spParent(0)); 00095 l2=EvtLeptonVACurrent(parent->getDaug(2)->spParentNeutrino(), 00096 parent->getDaug(1)->spParent(1)); 00097 } 00098 else{ 00099 report(ERROR,"EvtGen") << "dfnb89agngri wrong lepton number\n"; 00100 } 00101 } 00102 00103 amp.vertex(0,l1*tds); 00104 amp.vertex(1,l2*tds); 00105 00106 }
double EvtSemiLeptonicAmp::CalcMaxProb | ( | EvtId | parent, | |
EvtId | meson, | |||
EvtId | lepton, | |||
EvtId | nudaug, | |||
EvtSemiLeptonicFF * | FormFactors | |||
) | [inherited] |
Definition at line 38 of file EvtSemiLeptonicAmp.cc.
References boostTo(), EvtSemiLeptonicAmp::CalcAmp(), EvtParticle::deleteTree(), EvtParticle::getDaug(), EvtPDL::getMass(), EvtPDL::getMaxMass(), EvtPDL::getMeanMass(), EvtPDL::getMinMass(), EvtAmp::getSpinDensity(), EvtParticle::getSpinStates(), EvtPDL::getWidth(), genRecEmupikp::i, EvtParticle::init(), EvtAmp::init(), EvtScalarParticle::init(), ganga-rec::j, EvtParticle::makeDaughters(), EvtParticle::mass(), mass, EvtParticle::noLifeTime(), EvtSpinDensity::NormalizedProb(), EvtVector4R::set(), EvtSpinDensity::SetDiag(), and EvtParticle::setDiagonalSpinDensity().
Referenced by EvtSLPole::initProbMax(), EvtSLBKPole::initProbMax(), EvtHQET2::initProbMax(), and EvtHQET::initProbMax().
00040 { 00041 00042 //This routine takes the arguements parent, meson, and lepton 00043 //number, and a form factor model, and returns a maximum 00044 //probability for this semileptonic form factor model. A 00045 //brute force method is used. The 2D cos theta lepton and 00046 //q2 phase space is probed. 00047 00048 //Start by declaring a particle at rest. 00049 00050 //It only makes sense to have a scalar parent. For now. 00051 //This should be generalized later. 00052 00053 EvtScalarParticle *scalar_part; 00054 EvtParticle *root_part; 00055 00056 scalar_part=new EvtScalarParticle; 00057 00058 //cludge to avoid generating random numbers! 00059 scalar_part->noLifeTime(); 00060 00061 EvtVector4R p_init; 00062 00063 p_init.set(EvtPDL::getMass(parent),0.0,0.0,0.0); 00064 scalar_part->init(parent,p_init); 00065 root_part=(EvtParticle *)scalar_part; 00066 // root_part->set_type(EvtSpinType::SCALAR); 00067 root_part->setDiagonalSpinDensity(); 00068 00069 EvtParticle *daughter, *lep, *trino; 00070 00071 EvtAmp amp; 00072 00073 EvtId listdaug[3]; 00074 listdaug[0] = meson; 00075 listdaug[1] = lepton; 00076 listdaug[2] = nudaug; 00077 00078 amp.init(parent,3,listdaug); 00079 00080 root_part->makeDaughters(3,listdaug); 00081 daughter=root_part->getDaug(0); 00082 lep=root_part->getDaug(1); 00083 trino=root_part->getDaug(2); 00084 00085 //cludge to avoid generating random numbers! 00086 daughter->noLifeTime(); 00087 lep->noLifeTime(); 00088 trino->noLifeTime(); 00089 00090 00091 //Initial particle is unpolarized, well it is a scalar so it is 00092 //trivial 00093 EvtSpinDensity rho; 00094 rho.SetDiag(root_part->getSpinStates()); 00095 00096 double mass[3]; 00097 00098 double m = root_part->mass(); 00099 00100 EvtVector4R p4meson, p4lepton, p4nu, p4w; 00101 double q2max; 00102 00103 double q2, elepton, plepton; 00104 int i,j; 00105 double erho,prho,costl; 00106 00107 double maxfoundprob = 0.0; 00108 double prob = -10.0; 00109 int massiter; 00110 00111 for (massiter=0;massiter<3;massiter++){ 00112 00113 mass[0] = EvtPDL::getMeanMass(meson); 00114 mass[1] = EvtPDL::getMeanMass(lepton); 00115 mass[2] = EvtPDL::getMeanMass(nudaug); 00116 if ( massiter==1 ) { 00117 mass[0] = EvtPDL::getMinMass(meson); 00118 } 00119 if ( massiter==2 ) { 00120 mass[0] = EvtPDL::getMaxMass(meson); 00121 if ( (mass[0]+mass[1]+mass[2])>m) mass[0]=m-mass[1]-mass[2]-0.00001; 00122 } 00123 00124 q2max = (m-mass[0])*(m-mass[0]); 00125 00126 //loop over q2 00127 00128 for (i=0;i<25;i++) { 00129 q2 = ((i+0.5)*q2max)/25.0; 00130 00131 erho = ( m*m + mass[0]*mass[0] - q2 )/(2.0*m); 00132 00133 prho = sqrt(erho*erho-mass[0]*mass[0]); 00134 00135 p4meson.set(erho,0.0,0.0,-1.0*prho); 00136 p4w.set(m-erho,0.0,0.0,prho); 00137 00138 //This is in the W rest frame 00139 elepton = (q2+mass[1]*mass[1])/(2.0*sqrt(q2)); 00140 plepton = sqrt(elepton*elepton-mass[1]*mass[1]); 00141 00142 double probctl[3]; 00143 00144 for (j=0;j<3;j++) { 00145 00146 costl = 0.99*(j - 1.0); 00147 00148 //These are in the W rest frame. Need to boost out into 00149 //the B frame. 00150 p4lepton.set(elepton,0.0, 00151 plepton*sqrt(1.0-costl*costl),plepton*costl); 00152 p4nu.set(plepton,0.0, 00153 -1.0*plepton*sqrt(1.0-costl*costl),-1.0*plepton*costl); 00154 00155 EvtVector4R boost((m-erho),0.0,0.0,1.0*prho); 00156 p4lepton=boostTo(p4lepton,boost); 00157 p4nu=boostTo(p4nu,boost); 00158 00159 //Now initialize the daughters... 00160 00161 daughter->init(meson,p4meson); 00162 lep->init(lepton,p4lepton); 00163 trino->init(nudaug,p4nu); 00164 00165 CalcAmp(root_part,amp,FormFactors); 00166 00167 //Now find the probability at this q2 and cos theta lepton point 00168 //and compare to maxfoundprob. 00169 00170 //Do a little magic to get the probability!! 00171 prob = rho.NormalizedProb(amp.getSpinDensity()); 00172 00173 probctl[j]=prob; 00174 } 00175 00176 //probclt contains prob at ctl=-1,0,1. 00177 //prob=a+b*ctl+c*ctl^2 00178 00179 double a=probctl[1]; 00180 double b=0.5*(probctl[2]-probctl[0]); 00181 double c=0.5*(probctl[2]+probctl[0])-probctl[1]; 00182 00183 prob=probctl[0]; 00184 if (probctl[1]>prob) prob=probctl[1]; 00185 if (probctl[2]>prob) prob=probctl[2]; 00186 00187 if (fabs(c)>1e-20){ 00188 double ctlx=-0.5*b/c; 00189 if (fabs(ctlx)<1.0){ 00190 double probtmp=a+b*ctlx+c*ctlx*ctlx; 00191 if (probtmp>prob) prob=probtmp; 00192 } 00193 00194 } 00195 00196 //report(DEBUG,"EvtGen") << "prob,probctl:"<<prob<<" " 00197 // << probctl[0]<<" " 00198 // << probctl[1]<<" " 00199 // << probctl[2]<<endl; 00200 00201 if ( prob > maxfoundprob ) { 00202 maxfoundprob = prob; 00203 } 00204 00205 } 00206 if ( EvtPDL::getWidth(meson) <= 0.0 ) { 00207 //if the particle is narrow dont bother with changing the mass. 00208 massiter = 4; 00209 } 00210 00211 } 00212 root_part->deleteTree(); 00213 00214 maxfoundprob *=1.1; 00215 return maxfoundprob; 00216 00217 }