00001 //-------------------------------------------------------------------------- 00002 // 00003 // Environment: 00004 // This software is part of models developed at BES collaboration 00005 // based on the EvtGen framework. If you use all or part 00006 // of it, please give an appropriate acknowledgement. 00007 // 00008 // Copyright Information: See EvtGen/BesCopyright 00009 // Copyright (A) 2006 Ping Rong-Gang @IHEP 00010 // 00011 // Module: EvtMassH1.cc 00012 // 00013 // Description: Routine to decay a particle using a scatter 00014 // plot forn n-body decays (n>3). 00015 // 00016 // Modification history: 00017 // 00018 // Ping R.-G. Oct. 2011 Module created 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/EvtmH2.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 //#include "CLHEP/config/CLHEP.h" 00052 //#include "CLHEP/config/TemplateFunctions.h" 00053 00054 00055 using std::endl; 00056 00057 EvtmH2::~EvtmH2() {} 00058 00059 void EvtmH2::getName(std::string& model_name){ 00060 00061 model_name="mH2"; 00062 00063 } 00064 00065 EvtDecayBase* EvtmH2::clone(){ 00066 00067 return new EvtmH2; 00068 00069 } 00070 00071 00072 void EvtmH2::init(){ 00073 00074 // check that there are 4 arguments: Invariant mass part. Index: i,j, histor. file name, Hid 00075 checkNArg(2); 00076 EvtSpinType::spintype parenttype = EvtPDL::getSpinType(getParentId()); 00077 nbx = getArg(0); 00078 nby = getArg(1); 00079 } 00080 void EvtmH2::initProbMax(){ 00081 00082 noProbMax(); 00083 00084 } 00085 00086 void EvtmH2::decay( EvtParticle *p ){ 00087 00088 const char* fl=setFileName(); 00089 const char* hp=setHpoint(); 00090 00091 TFile f(fl); 00092 TH2F *hid = (TH2F*)f.Get("mH2"); 00093 TAxis* xaxis = hid->GetXaxis(); 00094 TAxis* yaxis = hid->GetYaxis(); 00095 00096 int BINSx =xaxis->GetLast(); 00097 int BINSy =yaxis->GetLast(); 00098 int BINS =BINSx*BINSy; 00099 double yvalue,ymax=0.0; 00100 int i,j,binxy; 00101 00102 for(i=1;i<BINSx+1;i++){ 00103 for(j=1;j<BINSy+1;j++){ 00104 binxy=hid->GetBin(i,j); 00105 yvalue=hid->GetBinContent(binxy); 00106 // cout<<"binxy,yvalue = "<<binxy<<"; "<<yvalue<<endl; 00107 if(yvalue>ymax) ymax=yvalue; 00108 } 00109 } 00110 00111 loop: 00112 p->initializePhaseSpace(getNDaug(),getDaugs()); 00113 00114 if(p->getNDaug()!= nbx+nby) {std::cout<<"The number of specified particles is not equal the number of decay daughters "<<endl;::abort();} 00115 00116 EvtVector4R pt,pt2; 00117 double xmass,ymass; 00118 00119 pt=p->getDaug(0)->getP4Lab(); 00120 for (int ii=1;ii<nbx;ii++){ 00121 pt=pt+p->getDaug(ii)->getP4Lab(); 00122 } 00123 xmass=pt.mass(); 00124 00125 pt2=p->getDaug(nbx)->getP4Lab(); 00126 for(int jj=nbx+1;jj<nbx+nby;jj++) pt2=pt2+p->getDaug(jj)->getP4Lab(); 00127 ymass=pt2.mass(); 00128 00129 00130 int xbin = hid->GetXaxis()->FindBin(xmass); 00131 int ybin = hid->GetYaxis()->FindBin(ymass); 00132 int xybin= hid->GetBin(xbin,ybin); 00133 double zvalue=hid->GetBinContent(xybin); 00134 double xratio=zvalue/ymax; 00135 // cout<<"***************zvalue,ymax,xratio= "<<zvalue<<"; "<<ymax<<"; "<<xratio<<endl; 00136 double rd1=EvtRandom::Flat(0.0, 1.0); 00137 if(rd1>xratio) goto loop; 00138 return ; 00139 } 00140 00141