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/EvtHAngSam3.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 EvtHAngSam3::~EvtHAngSam3() {}
00058
00059 void EvtHAngSam3::getName(std::string& model_name){
00060
00061 model_name="HAngSam3";
00062
00063 }
00064
00065 EvtDecayBase* EvtHAngSam3::clone(){
00066
00067 return new EvtHAngSam3;
00068
00069 }
00070
00071
00072 void EvtHAngSam3::init(){
00073
00074
00075 checkNArg(0);
00076 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId());
00077 }
00078 void EvtHAngSam3::initProbMax(){
00079
00080 noProbMax();
00081
00082 }
00083
00084 void EvtHAngSam3::decay( EvtParticle *p ){
00085
00086 const char* fl=setFileName();
00087 const char* hp=setHpoint();
00088 int* dp;
00089
00090 TFile f(fl);
00091 TH1F *hid = (TH1F*)f.Get(hp);
00092 TAxis* xaxis = hid->GetXaxis();
00093
00094 double BLE =xaxis->GetBinLowEdge(1);
00095 int BINS =xaxis->GetLast();
00096 double yvalue[20000];
00097 static double ymax=0.0;
00098 int i;
00099 double Ntotal=0,yc[20000];
00100 static int icount=0;
00101
00102 if(icount==0){
00103 for(i=1;i<BINS+1;i++){
00104 yvalue[i]=hid->GetBinContent(i);
00105 if(yvalue[i]>ymax) ymax=yvalue[i];
00106 }
00107 icount++;
00108 ymax=ymax*1.2;
00109 }
00110
00111 loop:
00112 p->initializePhaseSpace(getNDaug(),getDaugs());
00113
00114 EvtParticle *d1,*d2,*d3;
00115 EvtVector4R pd1,pd2,pd3,pp,nor;
00116 double costheta,cosphi;
00117
00118 d1 =p->getDaug(0);
00119 d2 =p->getDaug(1);
00120 pd1=d1->getP4();
00121 pd2=d2->getP4();
00122 pp =p->getP4();
00123 double ppmag=pp.d3mag();
00124 nor=pd1.cross(pd2);
00125
00126 if(ppmag<0.0001){costheta=nor.get(3)/nor.d3mag();} else
00127 {costheta=nor.dot(pp)/nor.d3mag()/pp.d3mag();}
00128
00129 int xbin =hid->FindBin(costheta);
00130 double xratio=(hid->GetBinContent(xbin))/ymax;
00131
00132 double rd1=EvtRandom::Flat(0.0, 1.0);
00133 if(rd1>xratio) goto loop;
00134 return ;
00135 }
00136
00137