00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015 #include "EvtGenBase/EvtHis2F.hh"
00016 #include <stdlib.h>
00017 #include <string>
00018 #include <fstream>
00019 #include <iostream>
00020
00021
00022 using namespace std;
00023 using std::endl;
00024
00025 EvtHis2F::~EvtHis2F(){}
00026
00027 const char* EvtHis2F::getFile(){
00028 return datafile;
00029 }
00030
00031 const char* EvtHis2F::getHTitle(){
00032 return datatitle;
00033 }
00034
00035 void EvtHis2F::setFile(const char* dtfile){
00036 datafile=dtfile;
00037 return;
00038 }
00039
00040 void EvtHis2F::setHTitle(const char* htitle){
00041 datatitle=htitle;
00042 return;
00043 }
00044
00045 void EvtHis2F::setHmc(TH2F *H2){
00046 HMC=(TH2F*) H2->Clone("HMC");
00047 return;
00048 }
00049
00050 void EvtHis2F::setHdt(TH2F *H2){
00051 HDATA=(TH2F*) H2->Clone("HDATA");
00052 return;
00053 }
00054
00055 TH2F* EvtHis2F::getHmc(){
00056 return HMC;
00057 }
00058
00059 TH2F* EvtHis2F::getHdt(){
00060 return HDATA;
00061 }
00062
00063 TH2F* EvtHis2F::getHwt(){
00064 return HWT;
00065 }
00066
00067 void EvtHis2F::show(TH2F *h2){
00068 TAxis* Xaxis=h2->GetXaxis();
00069 TAxis* Yaxis=h2->GetYaxis();
00070
00071 int bx =Xaxis->GetLast();
00072 int by =Yaxis->GetLast();
00073
00074 for(int i=1;i<bx+1;i++){
00075 for(int j=1;j<by+1;j++){
00076 int ij=h2->GetBin(i,j);
00077 double xij=h2->GetBinContent(ij);
00078 std::cout<<"i,j,cell,value "<<i<<" "<<j<<" "<<ij<<" "<<xij<<std::endl;
00079 }
00080 }
00081 }
00082
00083 void EvtHis2F::showFrame(TH2F *h2){
00084 TAxis* Xaxis=h2->GetXaxis();
00085 TAxis* Yaxis=h2->GetYaxis();
00086
00087 int bx =Xaxis->GetLast();
00088 int by =Yaxis->GetLast();
00089
00090 double x1 =Xaxis->GetBinLowEdge(1);
00091 double y1 =Yaxis->GetBinLowEdge(1);
00092 double x2 =Xaxis->GetBinUpEdge(bx);
00093 double y2 =Yaxis->GetBinUpEdge(by);
00094 std::cout<<"N_x, x_low, x_up; N_y, y_low, y_up "
00095 <<bx<<" "<<x1<<" "<<x2<<" "<<by<<" "<<y1<<" "<<y2<<std::endl;
00096 }
00097
00098 void EvtHis2F::init(){
00099 const char* dtfile =getFile();
00100 const char* htitle =getHTitle();
00101
00102 TFile dataf(dtfile);
00103 HDATA = (TH2F*)dataf.Get(htitle);
00104
00105 double ndata=HDATA->Integral();
00106 std::cout<<"Number events in HDATA= "<<ndata<<std::endl;
00107
00108 TAxis* Xaxis=HDATA->GetXaxis();
00109 TAxis* Yaxis=HDATA->GetYaxis();
00110
00111 BINSx =Xaxis->GetLast();
00112 BINSy =Yaxis->GetLast();
00113
00114 xlow=Xaxis->GetBinLowEdge(1);
00115 ylow=Yaxis->GetBinLowEdge(1);
00116 xup =Xaxis->GetBinUpEdge(BINSx);
00117 yup =Yaxis->GetBinUpEdge(BINSy);
00118
00119 std::cout<<"BINSx,xlow,xup,BINSy,ylow,yup: "<<BINSx<<" ,"<<xlow<<", "<<xup<<", "<<BINSy<<", "<<ylow<<", "<<yup<<std::endl;
00120 HMC = new TH2F("myHMC","",BINSx,xlow,xup,BINSy,ylow,yup);
00121 HWT = new TH2F("myHWT","",BINSx,xlow,xup,BINSy,ylow,yup);
00122 HMC ->SetDirectory(0);
00123 HWT ->SetDirectory(0);
00124 HDATA->SetDirectory(0);
00125 icount=0;
00126 return;
00127 }
00128
00129 void EvtHis2F::HFill(double m12s,double m13s){
00130 HMC->Fill(m12s,m13s,1.0);
00131 icount++;
00132
00133 if(icount==100000) showFrame(HMC);
00134 return;
00135 }
00136
00137 void EvtHis2F::HReweight(){
00138 std::cout<<"Now I reweight the MC to the data Dalitz Plot"<<std::endl;
00139 double ndata=HDATA->Integral();
00140 std::cout<<"Number events in Dalitz plot = "<<ndata<<std::endl;
00141 double nmc=HMC->Integral();
00142 std::cout<<"Number of events in HMC reweight= "<<nmc<<std::endl;
00143 HWT->Divide(HDATA,HMC);
00144 ndata=HWT->Integral();
00145 std::cout<<"Number of events after reweight= "<<ndata<<std::endl;
00146 double norm=nmc/ndata;
00147 HWT->Scale(norm);
00148 nmc=HMC->Integral();
00149 std::cout<<"Number of events in HMC after scale= "<<nmc<<std::endl;
00150 }
00151
00152 void EvtHis2F::setBins(TH2F *h2){
00153 TAxis* Xaxis=h2->GetXaxis();
00154 TAxis* Yaxis=h2->GetYaxis();
00155 BINSx =Xaxis->GetLast();
00156 BINSy =Yaxis->GetLast();
00157 xlow=Xaxis->GetBinLowEdge(1);
00158 ylow=Yaxis->GetBinLowEdge(1);
00159 xup =Xaxis->GetBinUpEdge(BINSx);
00160 yup =Yaxis->GetBinUpEdge(BINSy);
00161 }
00162 void EvtHis2F::setZmax(TH2F *hwt){
00163 double ymax=0;
00164 setBins(hwt);
00165
00166 for(int dx=1;dx <= BINSx;dx++){
00167 for(int dy=1;dy <= BINSy;dy++){
00168 int nxy=hwt->GetBin(dx,dy);
00169 double ncell=hwt->GetBinContent(nxy);
00170 if(ncell<=0) continue;
00171 if(ncell>ymax){ymax=ncell;}
00172 }
00173 }
00174 zmax=ymax;
00175 }
00176
00177 void EvtHis2F::setZmax(){
00178 double ymax=0;
00179 setBins(HWT);
00180
00181 for(int dx=1;dx <= BINSx;dx++){
00182 for(int dy=1;dy <= BINSy;dy++){
00183 int nxy=HWT->GetBin(dx,dy);
00184 double ncell=HWT->GetBinContent(nxy);
00185 if(ncell<=0) continue;
00186 if(ncell>ymax){ymax=ncell;}
00187 }
00188 }
00189 zmax=ymax;
00190 }
00191
00192 double EvtHis2F::getZmax(){
00193 return zmax;
00194 }
00195
00196 double EvtHis2F::getZvalue(double xmass2,double ymass2){
00197 int xbin = HWT->GetXaxis()->FindBin(xmass2);
00198 int ybin = HWT->GetYaxis()->FindBin(ymass2);
00199 int xybin= HWT->GetBin(xbin,ybin);
00200 double zvalue=HWT->GetBinContent(xybin);
00201
00202 return zvalue;
00203 }
00204
00205 bool EvtHis2F::AR(double xmass2,double ymass2){
00206
00207 bool accept;
00208 accept=false;
00209 if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) {return accept;}
00210 double zvalue=getZvalue(xmass2,ymass2);
00211 double ratio=zvalue/zmax;
00212 double rndm=EvtRandom::Flat(0.0, 1.0);
00213
00214 if(ratio > rndm){accept=true;} else {accept=false;}
00215 return accept;
00216 }
00217
00218 bool EvtHis2F::AR(double xmass2,double ymass2, double zmax, TH2F *h2){
00219
00220 bool accept;
00221 accept=false;
00222 setBins(h2);
00223 if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) {return accept;}
00224 int xbin = h2->GetXaxis()->FindBin(xmass2);
00225 int ybin = h2->GetYaxis()->FindBin(ymass2);
00226 int xybin= h2->GetBin(xbin,ybin);
00227 double zvalue=h2->GetBinContent(xybin);
00228
00229 double ratio=zvalue/zmax;
00230 double rndm=EvtRandom::Flat(0.0, 1.0);
00231
00232 if(ratio > rndm){accept=true;} else {accept=false;}
00233 return accept;
00234 }
00235
00236 void EvtHis2F::init(TH2F *hmc, TH2F *hdt, TH2F *hwt){
00237 setHmc(hmc);
00238 setHdt(hdt);
00239 setHwt(hwt);
00240 HReweight();
00241 setZmax();
00242 }
00243