00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021 #include "EvtGenBase/EvtPatches.hh"
00022 #include <stdlib.h>
00023 #include "EvtGenBase/EvtParticle.hh"
00024 #include "EvtGenBase/EvtGenKine.hh"
00025 #include "EvtGenBase/EvtPDL.hh"
00026 #include "EvtGenBase/EvtReport.hh"
00027 #include "EvtGenModels/EvtCalHelAmp.hh"
00028 #include "EvtGenModels/EvtGlobalSet.hh"
00029 #include <string>
00030
00031 double EvtCalHelAmp::_H2err=0;
00032 double EvtCalHelAmp::_H1err=0;
00033 double EvtCalHelAmp::_H0err=0;
00034 double EvtCalHelAmp::_H2=0;
00035 double EvtCalHelAmp::_H1=0;
00036 double EvtCalHelAmp::_H0=0;
00037 int EvtCalHelAmp::nevt=0;
00038 EvtCalHelAmp::~EvtCalHelAmp() {}
00039
00040 void EvtCalHelAmp::getName(std::string& model_name){
00041
00042 model_name="CALHELAMP";
00043
00044 }
00045
00046 EvtDecayBase* EvtCalHelAmp::clone(){
00047
00048 return new EvtCalHelAmp;
00049
00050 }
00051
00052
00053 void EvtCalHelAmp::init(){
00054 _H2=0;_H1=0;_H0=0;_H2err=0;_H1err=0;_H0err=0;nevt=0;
00055
00056 checkNArg(4);
00057 nd=getNArg();
00058 VC.resize(nd);
00059 for(int i=0;i<nd;i++) VC[i].resize(nd);
00060 ifstream coverrf("myerr.dat");
00061 if(coverrf==0) {std::cout<<"File of covariant error (myerr.dat)" <<" is not found"<<endl;abort();}
00062 for(int i=0;i<nd;i++){
00063 for(int j=0;j<nd;j++){
00064 double xr;
00065 coverrf>>xr;
00066 VC[i][j]=xr;
00067 }
00068 }
00069
00070 }
00071
00072 void EvtCalHelAmp::initProbMax(){
00073
00074 noProbMax();
00075
00076 }
00077
00078 void EvtCalHelAmp::decay( EvtParticle *p ){
00079 if(getNArg()!=4) {cout<<"The model have 6 parameters: |g00| phi00 |g22| phi22 |g42| phi42"<<endl; abort();}
00080
00081 double weight = p->initializePhaseSpace(getNDaug(),getDaugs());
00082
00083
00084 std::string p0str=EvtPDL::name(p->getId());
00085 std::string p1str=EvtPDL::name(p->getDaug(0)->getId());
00086 std::string p2str=EvtPDL::name(p->getDaug(1)->getId());
00087 if(p1str=="K_2*+" || p2str=="K_2*+" ||p1str=="K_2*0" || p2str=="K_2*0" ) {flag=2;} else
00088 if(p1str=="K_3*+" || p2str=="K_3*+" ||p1str=="K_3*0" || p2str=="K_3*0" ) {flag=3;} else
00089 if(p1str=="Zc+" || p2str=="pi-" ||p1str=="Zc-" || p2str=="pi^+" ) {flag=4;} else
00090 if(p0str=="Zc+" && (p1str=="J/psi" || p2str=="J/psi") ) {flag=5;}
00091
00092 EvtVector4R p1= p->getDaug(0)->getP4();
00093 double r= 2*p1.d3mag();
00094 double R=3.0;
00095 double pi=3.1415926;
00096 EvtBlattWeisskopf B0(0,R,r/2); double b0=B0.get(r/2);
00097 EvtBlattWeisskopf B1(1,R,r/2); double b1=B1.get(r/2);
00098 EvtBlattWeisskopf B2(2,R,r/2); double b2=B2.get(r/2);
00099 EvtBlattWeisskopf B3(3,R,r/2); double b3=B3.get(r/2);
00100 EvtBlattWeisskopf B4(4,R,r/2); double b4=B4.get(r/2);
00101 double g02 = getArg(0);
00102 double phi02 = getArg(1)*2*pi;
00103 double g22 = getArg(2);
00104 double phi22 = getArg(3)*2*pi;
00105
00106
00107 EvtComplex G02(g02*cos(phi02),g02*sin(phi02));
00108 EvtComplex G22(g22*cos(phi22),g22*sin(phi22));
00109
00110 std::vector<EvtComplex> VG; VG.clear();
00111 VG.push_back(G02); VG.push_back(G22);
00112 std::vector<double> VH2,VH1,VH0;
00113 VH2.resize(nd/2);VH1.resize(nd/2);VH0.resize(nd/2);
00114 if(flag==2){
00115 VH2[0]=sqrt(2./5.)*b1*r; VH2[1]= 1/sqrt(10.)*r*r*r*b3;
00116 VH1[0]=-1/sqrt(10.)*b1*r; VH1[1]= r*r*r*b3*sqrt(2./5.);
00117 }else if(flag==3){
00118 VH2[0]=sqrt(5./14.)*b2*r*r; VH2[1]= 1/sqrt(7.)*r*r*r*r*b4;
00119 VH1[0]=sqrt(1./7.)*b2*r*r; VH1[1]= -sqrt(5./14.)*r*r*r*r*b4;
00120 }else if(flag==4 || flag==5){
00121 VH2[0]=sqrt(1./3.)*b0; VH2[1]= 1/sqrt(6.)*r*r*b2;
00122 VH1[0]=sqrt(1./3.)*b0; VH1[1]=-2/sqrt(6.)*r*r*b2;
00123 }else{std::cout<<"Not allowed mode!"<<std::endl; abort();}
00124 EvtComplex H2 = G02*VH2[0] + G22*VH2[1];
00125 EvtComplex H1 = G02*VH1[0] + G22*VH1[1];
00126
00127 _H2 += abs2(H2); _H1 += abs2(H1);
00128 std::vector<double> DH1,DH2;
00129 DH2=firstder(VH2);
00130 DH1=firstder(VH1);
00131
00132
00133
00134
00135 for(int i=0;i<nd;i++){
00136 for(int j=0;j<nd;j++){
00137 _H2err += DH2[i]*DH2[j]*VC[i][j];
00138 _H1err += DH1[i]*DH1[j]*VC[i][j];
00139 }
00140 }
00141
00142 nevt++;
00143
00144 return ;
00145 }
00146
00147
00148 std::vector<double> EvtCalHelAmp::firstder(std::vector<double> vh){
00149
00150 double pi=3.1415926;
00151 std::vector<double> fd;
00152 double g02 = getArg(0);
00153 double phi02 = getArg(1)*2*pi;
00154 double g22 = getArg(2);
00155 double phi22 = getArg(3)*2*pi;
00156 std::vector<double> vpar;
00157 vpar.push_back(g02);vpar.push_back(phi02);vpar.push_back(g22);vpar.push_back(phi22);
00158 EvtComplex G02(vpar[0]*cos(vpar[1]),vpar[0]*sin(vpar[1]));
00159 EvtComplex G22(vpar[2]*cos(vpar[3]),vpar[2]*sin(vpar[3]));
00160 EvtComplex H20 = G02*vh[0] + G22*vh[1];
00161 double hamps0 = abs2(H20);
00162 for(int i=0;i<vpar.size();i++){
00163 vpar.clear();
00164 vpar.push_back(g02);vpar.push_back(phi02);vpar.push_back(g22);vpar.push_back(phi22);
00165 double dev=0.01;
00166 vpar[i] += dev;
00167 EvtComplex G02(vpar[0]*cos(vpar[1]),vpar[0]*sin(vpar[1]));
00168 EvtComplex G22(vpar[2]*cos(vpar[3]),vpar[2]*sin(vpar[3]));
00169 EvtComplex H20 = G02*vh[0] + G22*vh[1];
00170 double hamps2=abs2(H20);
00171 double xder=(hamps2-hamps0)/dev;
00172 fd.push_back(xder);
00173 }
00174 return fd;
00175 }