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 <stdlib.h>
00024 #include <iostream>
00025 #include <string>
00026 #include "EvtGenBase/EvtParticle.hh"
00027 #include "EvtGenBase/EvtPDL.hh"
00028 #include "EvtGenBase/EvtGenKine.hh"
00029 #include "EvtGenModels/EvtLNuGamma.hh"
00030 #include "EvtGenBase/EvtDiracSpinor.hh"
00031 #include "EvtGenBase/EvtReport.hh"
00032 #include "EvtGenBase/EvtComplex.hh"
00033 #include "EvtGenBase/EvtVector4C.hh"
00034 #include "EvtGenBase/EvtVector4R.hh"
00035 #include "EvtGenBase/EvtTensor4C.hh"
00036
00037 EvtLNuGamma::EvtLNuGamma(){
00038 _fafvzero = false;
00039 }
00040
00041 EvtLNuGamma::~EvtLNuGamma() {}
00042
00043 void EvtLNuGamma::getName(std::string& model_name){
00044
00045 model_name="LNUGAMMA";
00046
00047 }
00048
00049
00050 EvtDecayBase* EvtLNuGamma::clone(){
00051
00052 return new EvtLNuGamma;
00053
00054 }
00055
00056 void EvtLNuGamma::init(){
00057
00058
00059 checkNArg(3,4);
00060 checkNDaug(3);
00061
00062 if (getNArg() == 4){
00063
00064
00065 if (getArg(3) > 0){
00066 _fafvzero = true;
00067 }
00068 else{
00069 _fafvzero = false;
00070 }
00071 }
00072 else{
00073 _fafvzero = false;
00074 }
00075
00076 checkSpinParent(EvtSpinType::SCALAR);
00077
00078 checkSpinDaughter(0,EvtSpinType::DIRAC);
00079 checkSpinDaughter(1,EvtSpinType::NEUTRINO);
00080 checkSpinDaughter(2,EvtSpinType::PHOTON);
00081 }
00082
00083 void EvtLNuGamma::initProbMax(){
00084
00085 setProbMax(7000.0);
00086
00087 }
00088
00089 void EvtLNuGamma::decay(EvtParticle *p){
00090
00091 static EvtId BM=EvtPDL::getId("B-");
00092 p->initializePhaseSpace(getNDaug(),getDaugs());
00093
00094 EvtComplex myI(0,1);
00095
00096 EvtParticle *lept, *neut,*phot;
00097 lept=p->getDaug(0);
00098 neut=p->getDaug(1);
00099 phot=p->getDaug(2);
00100
00101 EvtVector4C lept1,lept2,photon1,photon2;
00102
00103 if (p->getId()==BM) {
00104 lept1=EvtLeptonVACurrent(lept->spParent(0),neut->spParentNeutrino());
00105 lept2=EvtLeptonVACurrent(lept->spParent(1),neut->spParentNeutrino());
00106 }
00107 else{
00108 lept1=EvtLeptonVACurrent(neut->spParentNeutrino(),lept->spParent(0));
00109 lept2=EvtLeptonVACurrent(neut->spParentNeutrino(),lept->spParent(1));
00110 }
00111
00112
00113
00114 EvtVector4R photp = phot->getP4();
00115 double photE = photp.get(0);
00116
00117 EvtVector4C photone1 = phot->epsParentPhoton(0).conj();
00118 EvtVector4C photone2 = phot->epsParentPhoton(1).conj();
00119
00120 EvtVector4R parVelocity(1,0,0,0);
00121
00122 double fv,fa;
00123
00124 fv = getFormFactor(photE);
00125 if (_fafvzero){
00126 fa = 0.0;
00127 }
00128 else if (p->getId()==BM) {
00129 fa = - fv;
00130 }
00131 else{
00132 fa = fv;
00133 }
00134
00135 EvtVector4C temp1a = dual(directProd(parVelocity,photp)).cont2(photone1);
00136 EvtVector4C temp2a = dual(directProd(parVelocity,photp)).cont2(photone2);
00137
00138 EvtVector4C temp1b = (photone1)*(parVelocity*photp);
00139 EvtVector4C temp1c = (photp)*(photone1*parVelocity);
00140
00141 EvtVector4C temp2b = (photone2)*(parVelocity*photp);
00142 EvtVector4C temp2c = (photp)*(photone2*parVelocity);
00143
00144 photon1 = (temp1a*fv) + (myI*fa*(temp1b - temp1c));
00145 photon2 = (temp2a*fv) + (myI*fa*(temp2b - temp2c));
00146
00147 vertex(0,0,lept1.cont(photon1));
00148 vertex(0,1,lept1.cont(photon2));
00149 vertex(1,0,lept2.cont(photon1));
00150 vertex(1,1,lept2.cont(photon2));
00151
00152 return;
00153
00154 }
00155
00156 double EvtLNuGamma::getFormFactor(double photonEnergy){
00157
00158
00159
00160
00161
00162
00163 double formFactor = 0;
00164 double qu = 2./3.;
00165 double qb = -1./3.;
00166
00167 if (photonEnergy > getArg(0)){
00168 formFactor = (1/photonEnergy) * ((qu*getArg(1)) - (qb/getArg(2)));
00169 }
00170 return formFactor;
00171 }