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 "EvtGenBase/EvtParticle.hh"
00025 #include "EvtGenBase/EvtGenKine.hh"
00026 #include "EvtGenBase/EvtPDL.hh"
00027 #include "EvtGenBase/EvtVector4C.hh"
00028 #include "EvtGenBase/EvtVector4R.hh"
00029 #include "EvtGenBase/EvtTensor4C.hh"
00030 #include "EvtGenBase/EvtDiracParticle.hh"
00031 #include "EvtGenBase/EvtScalarParticle.hh"
00032 #include "EvtGenBase/EvtVectorParticle.hh"
00033 #include "EvtGenBase/EvtTensorParticle.hh"
00034 #include "EvtGenBase/EvtPhotonParticle.hh"
00035 #include "EvtGenBase/EvtNeutrinoParticle.hh"
00036 #include "EvtGenBase/EvtStringParticle.hh"
00037 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
00038 #include "EvtGenBase/EvtHighSpinParticle.hh"
00039 #include "EvtGenBase/EvtReport.hh"
00040 #include "EvtGenBase/EvtHelSys.hh"
00041 #include "EvtGenModels/EvtBody3.hh"
00042 #include "EvtGenBase/EvtRandom.hh"
00043 #include <string>
00044
00045 #include "TH1.h"
00046 #include "TAxis.h"
00047 #include "TH2.h"
00048 #include "TFile.h"
00049 #include "TApplication.h"
00050 #include "TROOT.h"
00051
00052
00053
00054
00055 using std::endl;
00056
00057 EvtBody3::~EvtBody3() {}
00058
00059 void EvtBody3::getName(std::string& model_name){
00060
00061 model_name="Body3";
00062
00063 }
00064
00065 EvtDecayBase* EvtBody3::clone(){
00066
00067 return new EvtBody3;
00068
00069 }
00070
00071
00072 void EvtBody3::init(){
00073
00074
00075 checkNArg(0);
00076 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00077 }
00078 void EvtBody3::initProbMax(){
00079
00080 noProbMax();
00081
00082 }
00083
00084 void EvtBody3::decay( EvtParticle *p ){
00085
00086 const char* fl=setFileName();
00087 const char* hp=setHpoint();
00088 int* dp;
00089 int nd1,nd2,ii,sn;
00090
00091
00092 sn=setDaugAngNo();
00093
00094 if(sn==2){
00095 nd1=0;
00096 nd2=1;
00097 }
00098
00099 if(sn==0){
00100 nd1=1;
00101 nd2=2;
00102 }
00103
00104 if(sn==1){
00105 nd1=0;
00106 nd2=2;
00107 }
00108 const char* DA1=setDaugAng(nd1);
00109 const char* DA2=setDaugAng(nd2);
00110
00111
00112 dp=setDaugPair();
00113 int d1=dp[0];
00114 int d2=dp[1];
00115 int d3=dp[2];
00116 int d4=dp[3];
00117
00118
00119 TFile f(fl);
00120 TH1F *did1= (TH1F*)f.Get(DA1);
00121 TH1F *did2= (TH1F*)f.Get(DA2);
00122 TH2F *hid = (TH2F*)f.Get(hp);
00123
00124 TAxis* d1x=did1->GetXaxis();
00125 TAxis* d2x=did2->GetXaxis();
00126 TAxis* xaxis = hid->GetXaxis();
00127 TAxis* yaxis = hid->GetYaxis();
00128
00129 int BINSd1 =did1->GetXaxis()->GetLast();
00130 int BINSd2 =did2->GetXaxis()->GetLast();
00131 int BINSx =xaxis->GetLast();
00132 int BINSy =yaxis->GetLast();
00133 int BINS =BINSx*BINSy;
00134
00135 double av1,av2,avm1,avm2;
00136 avm1=0.;
00137 avm2=0.;
00138 double yvalue,ymax=0.0;
00139 int i,j,binxy;
00140
00141
00142 for(i=1;i<BINSd1+1;i++){
00143 av1=did1->GetBinContent(i);
00144 if(av1>avm1) avm1=av1;
00145 }
00146
00147
00148 for(i=1;i<BINSd2+1;i++){
00149 av2=did2->GetBinContent(i);
00150 if(av2>avm2) avm2=av2;
00151 }
00152
00153
00154
00155 for(i=1;i<BINSx+1;i++){
00156 for(j=1;j<BINSy+1;j++){
00157 binxy=hid->GetBin(i,j);
00158 yvalue=hid->GetBinContent(binxy);
00159
00160 if(yvalue>ymax) ymax=yvalue;
00161 }
00162 }
00163
00164 loop:
00165 p->initializePhaseSpace(getNDaug(),getDaugs());
00166
00167 EvtParticle *id1,*id2,*id3,*id4,*s1;
00168 EvtVector4R pd1,pd2,pd3,pd4,ps;
00169 EvtVector4R dp1,dp2;
00170 double xmass2,ymass2;
00171
00172 id1 =p->getDaug(d1);
00173 id2 =p->getDaug(d2);
00174 id3 =p->getDaug(d3);
00175 id4 =p->getDaug(d4);
00176
00177 pd1 =id1->getP4Lab();
00178 pd2 =id2->getP4Lab();
00179 pd3 =id3->getP4Lab();
00180 pd4 =id4->getP4Lab();
00181
00182 dp1 =p->getDaug(nd1)->getP4Lab();
00183 dp2 =p->getDaug(nd2)->getP4Lab();
00184
00185 xmass2=(pd1+pd2).mass2();
00186 ymass2=(pd3+pd4).mass2();
00187
00188 int xbin = hid->GetXaxis()->FindBin(xmass2);
00189 int ybin = hid->GetYaxis()->FindBin(ymass2);
00190 int xybin= hid->GetBin(xbin,ybin);
00191 double zvalue=hid->GetBinContent(xybin);
00192 double xratio=zvalue/ymax;
00193 if(xratio ==0 ) goto loop;
00194 double rd1=EvtRandom::Flat(0.0, 1.0);
00195 if(rd1>xratio) goto loop;
00196
00197 double da1=dp1.get(3)/dp1.d3mag();
00198 double da2=dp2.get(3)/dp2.d3mag();
00199
00200 int dbin1=did1->FindBin(da1);
00201 int dbin2=did2->FindBin(da2);
00202
00203 double dr1=(did1->GetBinContent(dbin1))/avm1;
00204 double dr2=(did2->GetBinContent(dbin2))/avm2;
00205 if(dr1==0 || dr2==0) goto loop;
00206 rd1=EvtRandom::Flat(0.0, 1.0);
00207 if(rd1>dr1) goto loop;
00208
00209 rd1=EvtRandom::Flat(0.0, 1.0);
00210 if(rd1>dr2) goto loop;
00211
00212 return ;
00213 }
00214
00215