00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 #include "EvtGenBase/EvtPatches.hh"
00025 #include <stdlib.h>
00026 #include "EvtGenBase/EvtParticle.hh"
00027 #include "EvtGenBase/EvtGenKine.hh"
00028 #include "EvtGenBase/EvtPDL.hh"
00029 #include "EvtGenBase/EvtVector4C.hh"
00030 #include "EvtGenBase/EvtVector4R.hh"
00031 #include "EvtGenBase/EvtTensor4C.hh"
00032 #include "EvtGenBase/EvtDiracParticle.hh"
00033 #include "EvtGenBase/EvtScalarParticle.hh"
00034 #include "EvtGenBase/EvtVectorParticle.hh"
00035 #include "EvtGenBase/EvtTensorParticle.hh"
00036 #include "EvtGenBase/EvtPhotonParticle.hh"
00037 #include "EvtGenBase/EvtNeutrinoParticle.hh"
00038 #include "EvtGenBase/EvtStringParticle.hh"
00039 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
00040 #include "EvtGenBase/EvtHighSpinParticle.hh"
00041 #include "EvtGenBase/EvtReport.hh"
00042 #include "EvtGenBase/EvtHelSys.hh"
00043 #include "EvtGenModels/EvtAngH2.hh"
00044 #include "EvtGenBase/EvtRandom.hh"
00045 #include <string>
00046
00047 #include "TH1.h"
00048 #include "TAxis.h"
00049 #include "TH2.h"
00050 #include "TFile.h"
00051 #include "TApplication.h"
00052 #include "TROOT.h"
00053
00054
00055
00056
00057 using std::endl;
00058
00059 EvtAngH2::~EvtAngH2() {}
00060
00061 void EvtAngH2::getName(std::string& model_name){
00062
00063 model_name="AngH2";
00064
00065 }
00066
00067 EvtDecayBase* EvtAngH2::clone(){
00068
00069 return new EvtAngH2;
00070
00071 }
00072
00073
00074 void EvtAngH2::init(){
00075
00076
00077 checkNArg(0);
00078 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00079 }
00080 void EvtAngH2::initProbMax(){
00081
00082 noProbMax();
00083
00084 }
00085
00086 void EvtAngH2::decay( EvtParticle *p ){
00087
00088 const char* fl=setFileName();
00089 const char* hp=setHpoint();
00090 int* dp;
00091
00092 TFile f(fl);
00093 TH2F *hid = (TH2F*)f.Get(hp);
00094 TAxis* xaxis = hid->GetXaxis();
00095 TAxis* yaxis = hid->GetYaxis();
00096
00097 int BINSx =xaxis->GetLast();
00098 int BINSy =yaxis->GetLast();
00099 int BINS =BINSx*BINSy;
00100 double yvalue,ymax=0.0;
00101 int i,j,binxy;
00102
00103 for(i=1;i<BINSx+1;i++){
00104 for(j=1;j<BINSy+1;j++){
00105 binxy=hid->GetBin(i,j);
00106 yvalue=hid->GetBinContent(binxy);
00107
00108 if(yvalue>ymax) ymax=yvalue;
00109 }
00110 }
00111
00112 loop:
00113 p->initializePhaseSpace(getNDaug(),getDaugs());
00114
00115 EvtParticle *id1,*id2,*id3,*id4,*s1;
00116 EvtVector4R pd1,pd2,pd3,pd4,ps;
00117 double xcostheta,ycostheta;
00118
00119 id1 =p->getDaug(0);
00120 id2 =p->getDaug(1);
00121 id3 =p->getDaug(2);
00122
00123
00124 pd1 =id1->getP4Lab();
00125 pd2 =id2->getP4Lab();
00126 pd3 =id3->getP4Lab();
00127
00128 xcostheta=pd1.get(3)/pd1.d3mag();
00129 ycostheta=pd2.get(3)/pd2.d3mag();
00130 if( EvtPDL::name(p->getDaug(0)->getId()) == EvtPDL::name(p->getDaug(1)->getId()) ){
00131 if(pd1.get(0)>pd2.get(0)){
00132 xcostheta=pd2.get(3)/pd2.d3mag();
00133 ycostheta=pd1.get(3)/pd1.d3mag();
00134 }
00135 }
00136 int xbin = hid->GetXaxis()->FindBin(xcostheta);
00137 int ybin = hid->GetYaxis()->FindBin(ycostheta);
00138 int xybin= hid->GetBin(xbin,ybin);
00139 double zvalue=hid->GetBinContent(xybin);
00140 double xratio=zvalue/ymax;
00141
00142 double rd1=EvtRandom::Flat(0.0, 1.0);
00143 if(rd1>xratio) goto loop;
00144 return ;
00145 }
00146
00147