00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022 #include "EvtGenBase/EvtPatches.hh"
00023 #include "EvtGenBase/EvtPatches.hh"
00024 #include "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtGenKine.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenBase/EvtReport.hh"
00028 #include "EvtGenBase/EvtVector4C.hh"
00029 #include "EvtGenBase/EvtTensor4C.hh"
00030 #include "EvtGenBase/EvtDiracSpinor.hh"
00031 #include "EvtGenBase/EvtSemiLeptonicScalarAmp.hh"
00032 #include "EvtGenBase/EvtId.hh"
00033 #include "EvtGenBase/EvtAmp.hh"
00034 #include "EvtGenBase/EvtScalarParticle.hh"
00035 #include "EvtGenBase/EvtVectorParticle.hh"
00036 #include "EvtGenBase/EvtTensorParticle.hh"
00037
00038 double EvtSemiLeptonicAmp::CalcMaxProb( EvtId parent, EvtId meson,
00039 EvtId lepton, EvtId nudaug,
00040 EvtSemiLeptonicFF *FormFactors ) {
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 EvtScalarParticle *scalar_part;
00054 EvtParticle *root_part;
00055
00056 scalar_part=new EvtScalarParticle;
00057
00058
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
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
00086 daughter->noLifeTime();
00087 lep->noLifeTime();
00088 trino->noLifeTime();
00089
00090
00091
00092
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
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
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
00149
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
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
00168
00169
00170
00171 prob = rho.NormalizedProb(amp.getSpinDensity());
00172
00173 probctl[j]=prob;
00174 }
00175
00176
00177
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
00197
00198
00199
00200
00201 if ( prob > maxfoundprob ) {
00202 maxfoundprob = prob;
00203 }
00204
00205 }
00206 if ( EvtPDL::getWidth(meson) <= 0.0 ) {
00207
00208 massiter = 4;
00209 }
00210
00211 }
00212 root_part->deleteTree();
00213
00214 maxfoundprob *=1.1;
00215 return maxfoundprob;
00216
00217 }
00218