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/EvtMassH2.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 EvtMassH2::~EvtMassH2() {}
00058
00059 void EvtMassH2::getName(std::string& model_name){
00060
00061 model_name="MassH2";
00062
00063 }
00064
00065 EvtDecayBase* EvtMassH2::clone(){
00066
00067 return new EvtMassH2;
00068
00069 }
00070
00071
00072 void EvtMassH2::init(){
00073
00074
00075 checkNArg(0);
00076 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00077 }
00078 void EvtMassH2::initProbMax(){
00079
00080 noProbMax();
00081
00082 }
00083
00084 void EvtMassH2::decay( EvtParticle *p ){
00085
00086 const char* fl=setFileName();
00087 const char* hp=setHpoint();
00088 int* dp;
00089 dp=setDaugPair();
00090 int d1=dp[0];
00091 int d2=dp[1];
00092 int d3=dp[2];
00093 int d4=dp[3];
00094
00095 TFile f(fl);
00096 TH2F *hid = (TH2F*)f.Get(hp);
00097 TAxis* xaxis = hid->GetXaxis();
00098 TAxis* yaxis = hid->GetYaxis();
00099
00100 int BINSx =xaxis->GetLast();
00101 int BINSy =yaxis->GetLast();
00102 int BINS =BINSx*BINSy;
00103 double yvalue,ymax=0.0;
00104 int i,j,binxy;
00105
00106 for(i=1;i<BINSx+1;i++){
00107 for(j=1;j<BINSy+1;j++){
00108 binxy=hid->GetBin(i,j);
00109 yvalue=hid->GetBinContent(binxy);
00110
00111 if(yvalue>ymax) ymax=yvalue;
00112 }
00113 }
00114
00115 loop:
00116 p->initializePhaseSpace(getNDaug(),getDaugs());
00117
00118 EvtParticle *id1,*id2,*id3,*id4,*s1;
00119 EvtVector4R pd1,pd2,pd3,pd4,ps;
00120 double xmass2,ymass2;
00121
00122 id1 =p->getDaug(d1);
00123 id2 =p->getDaug(d2);
00124 id3 =p->getDaug(d3);
00125 id4 =p->getDaug(d4);
00126
00127 pd1 =id1->getP4Lab();
00128 pd2 =id2->getP4Lab();
00129 pd3 =id3->getP4Lab();
00130 pd4 =id4->getP4Lab();
00131
00132 xmass2=(pd1+pd2).mass2();
00133 ymass2=(pd3+pd4).mass2();
00134
00135 int xbin = hid->GetXaxis()->FindBin(xmass2);
00136 int ybin = hid->GetYaxis()->FindBin(ymass2);
00137 int xybin= hid->GetBin(xbin,ybin);
00138 double zvalue=hid->GetBinContent(xybin);
00139 double xratio=zvalue/ymax;
00140
00141 double rd1=EvtRandom::Flat(0.0, 1.0);
00142 if(rd1>xratio) goto loop;
00143 return ;
00144 }
00145
00146