BOSS 7.0.3
BESIII Offline Software System
Loading...
Searching...
No Matches
EvtHis2F.cc
Go to the documentation of this file.
1//--------------------------------------------------------------------------
2//
3// Environment:
4//
5// Module: EvtGen/EvtHis2F.hh
6//
7// Description: Class to handle H2F histogram
8//
9// Modification history:
10//
11// PING RG September 8, 2010 Module created
12//
13//------------------------------------------------------------------------
14
16#include <stdlib.h>
17#include <string>
18#include <fstream>
19#include <iostream>
20
21
22 using namespace std;
23 using std::endl;
24
26
27const char* EvtHis2F::getFile(){
28 return datafile;
29}
30
31const char* EvtHis2F::getHTitle(){
32 return datatitle;
33}
34
35void EvtHis2F::setFile(const char* dtfile){
36 datafile=dtfile;
37 return;
38}
39
40void EvtHis2F::setHTitle(const char* htitle){
41 datatitle=htitle;
42 return;
43}
44
45void EvtHis2F::setHmc(TH2F *H2){
46 HMC=(TH2F*) H2->Clone("HMC");
47 return;
48}
49
50void EvtHis2F::setHdt(TH2F *H2){
51 HDATA=(TH2F*) H2->Clone("HDATA");
52 return;
53}
54
56 return HMC;
57}
58
60 return HDATA;
61}
62
64 return HWT;
65}
66
67void EvtHis2F::show(TH2F *h2){
68 TAxis* Xaxis=h2->GetXaxis();
69 TAxis* Yaxis=h2->GetYaxis();
70
71 int bx =Xaxis->GetLast();
72 int by =Yaxis->GetLast();
73
74 for(int i=1;i<bx+1;i++){
75 for(int j=1;j<by+1;j++){
76 int ij=h2->GetBin(i,j);
77 double xij=h2->GetBinContent(ij);
78 std::cout<<"i,j,cell,value "<<i<<" "<<j<<" "<<ij<<" "<<xij<<std::endl;
79 }
80 }
81}
82
83void EvtHis2F::showFrame(TH2F *h2){
84 TAxis* Xaxis=h2->GetXaxis();
85 TAxis* Yaxis=h2->GetYaxis();
86
87 int bx =Xaxis->GetLast();
88 int by =Yaxis->GetLast();
89
90 double x1 =Xaxis->GetBinLowEdge(1);
91 double y1 =Yaxis->GetBinLowEdge(1);
92 double x2 =Xaxis->GetBinUpEdge(bx);
93 double y2 =Yaxis->GetBinUpEdge(by);
94 std::cout<<"N_x, x_low, x_up; N_y, y_low, y_up "
95 <<bx<<" "<<x1<<" "<<x2<<" "<<by<<" "<<y1<<" "<<y2<<std::endl;
96}
97
99 const char* dtfile =getFile(); // get data root file as input
100 const char* htitle =getHTitle(); // get Dalitz plot title
101
102 TFile dataf(dtfile);
103 HDATA = (TH2F*)dataf.Get(htitle);
104
105 double ndata=HDATA->Integral();
106 std::cout<<"Number events in HDATA= "<<ndata<<std::endl;
107
108 TAxis* Xaxis=HDATA->GetXaxis();
109 TAxis* Yaxis=HDATA->GetYaxis();
110
111 BINSx =Xaxis->GetLast();
112 BINSy =Yaxis->GetLast();
113
114 xlow=Xaxis->GetBinLowEdge(1);
115 ylow=Yaxis->GetBinLowEdge(1);
116 xup =Xaxis->GetBinUpEdge(BINSx);
117 yup =Yaxis->GetBinUpEdge(BINSy);
118
119 std::cout<<"BINSx,xlow,xup,BINSy,ylow,yup: "<<BINSx<<" ,"<<xlow<<", "<<xup<<", "<<BINSy<<", "<<ylow<<", "<<yup<<std::endl;
120 HMC = new TH2F("myHMC","",BINSx,xlow,xup,BINSy,ylow,yup);
121 HWT = new TH2F("myHWT","",BINSx,xlow,xup,BINSy,ylow,yup);
122 HMC ->SetDirectory(0);
123 HWT ->SetDirectory(0);
124 HDATA->SetDirectory(0);
125 icount=0;
126 return;
127}
128
129void EvtHis2F::HFill(double m12s,double m13s){
130 HMC->Fill(m12s,m13s,1.0);
131 icount++;
132 // if(icount==100000)show(HMC);
133 if(icount==100000) showFrame(HMC);
134 return;
135}
136
137void EvtHis2F::HReweight(){ //reweight Dalitz Plot after get HDATA and HMC
138 std::cout<<"Now I reweight the MC to the data Dalitz Plot"<<std::endl;
139 double ndata=HDATA->Integral();
140 std::cout<<"Number events in Dalitz plot = "<<ndata<<std::endl;
141 double nmc=HMC->Integral();
142 std::cout<<"Number of events in HMC reweight= "<<nmc<<std::endl;
143 HWT->Divide(HDATA,HMC);
144 ndata=HWT->Integral();
145 std::cout<<"Number of events after reweight= "<<ndata<<std::endl;
146 double norm=nmc/ndata;
147 HWT->Scale(norm);
148 nmc=HMC->Integral();
149 std::cout<<"Number of events in HMC after scale= "<<nmc<<std::endl;
150}
151
152void EvtHis2F::setBins(TH2F *h2){
153 TAxis* Xaxis=h2->GetXaxis();
154 TAxis* Yaxis=h2->GetYaxis();
155 BINSx =Xaxis->GetLast();
156 BINSy =Yaxis->GetLast();
157 xlow=Xaxis->GetBinLowEdge(1);
158 ylow=Yaxis->GetBinLowEdge(1);
159 xup =Xaxis->GetBinUpEdge(BINSx);
160 yup =Yaxis->GetBinUpEdge(BINSy);
161}
162void EvtHis2F::setZmax(TH2F *hwt){
163 double ymax=0;
164 setBins(hwt);
165
166 for(int dx=1;dx <= BINSx;dx++){
167 for(int dy=1;dy <= BINSy;dy++){
168 int nxy=hwt->GetBin(dx,dy);
169 double ncell=hwt->GetBinContent(nxy);
170 if(ncell<=0) continue;
171 if(ncell>ymax){ymax=ncell;}
172 }
173 }
174 zmax=ymax;
175}
176
178 double ymax=0;
179 setBins(HWT);
180
181 for(int dx=1;dx <= BINSx;dx++){
182 for(int dy=1;dy <= BINSy;dy++){
183 int nxy=HWT->GetBin(dx,dy);
184 double ncell=HWT->GetBinContent(nxy);
185 if(ncell<=0) continue;
186 if(ncell>ymax){ymax=ncell;}
187 }
188 }
189 zmax=ymax;
190}
191
193 return zmax;
194}
195
196double EvtHis2F::getZvalue(double xmass2,double ymass2){
197 int xbin = HWT->GetXaxis()->FindBin(xmass2);
198 int ybin = HWT->GetYaxis()->FindBin(ymass2);
199 int xybin= HWT->GetBin(xbin,ybin);
200 double zvalue=HWT->GetBinContent(xybin);
201// std::cout<<"xmass, ymass, zvalue "<<xmass2<<" "<<ymass2<<" "<<zvalue<<std::endl;
202 return zvalue;
203}
204
205bool EvtHis2F::AR(double xmass2,double ymass2){//accept: true. reject: false
206 //call this function after initiation and filling MC histogram HMC
207 bool accept;
208 accept=false;
209 if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) {return accept;}
210 double zvalue=getZvalue(xmass2,ymass2);
211 double ratio=zvalue/zmax;
212 double rndm=EvtRandom::Flat(0.0, 1.0);
213 // std::cout<<"zvalue, zmax , rndm "<<zvalue<<" "<<zmax<<" "<<rndm<<std::endl;
214 if(ratio > rndm){accept=true;} else {accept=false;}
215 return accept;
216}
217
218bool EvtHis2F::AR(double xmass2,double ymass2, double zmax, TH2F *h2){//accept: true. reject: false
219 //call this function after initiation and filling MC histogram HMC
220 bool accept;
221 accept=false;
222 setBins(h2);
223 if(xmass2<xlow || xmass2>xup || ymass2<ylow || ymass2>yup) {return accept;}
224 int xbin = h2->GetXaxis()->FindBin(xmass2);
225 int ybin = h2->GetYaxis()->FindBin(ymass2);
226 int xybin= h2->GetBin(xbin,ybin);
227 double zvalue=h2->GetBinContent(xybin);
228
229 double ratio=zvalue/zmax;
230 double rndm=EvtRandom::Flat(0.0, 1.0);
231 // std::cout<<"zvalue, zmax , rndm "<<zvalue<<" "<<zmax<<" "<<rndm<<std::endl;
232 if(ratio > rndm){accept=true;} else {accept=false;}
233 return accept;
234}
235
236void EvtHis2F::init(TH2F *hmc, TH2F *hdt, TH2F *hwt){
237 setHmc(hmc);
238 setHdt(hdt);
239 setHwt(hwt);
240 HReweight();
241 setZmax();
242}
243
void show(TH2F *h2)
Definition: EvtHis2F.cc:67
const char * getHTitle()
Definition: EvtHis2F.cc:31
void HReweight()
Definition: EvtHis2F.cc:137
void setBins(TH2F *h2)
Definition: EvtHis2F.cc:152
void setHwt(TH2F *H2)
Definition: EvtHis2F.hh:50
TH2F * getHmc()
Definition: EvtHis2F.cc:55
void setHdt(TH2F *H2)
Definition: EvtHis2F.cc:50
void init()
Definition: EvtHis2F.cc:98
virtual ~EvtHis2F()
Definition: EvtHis2F.cc:25
double getZmax()
Definition: EvtHis2F.cc:192
const char * getFile()
Definition: EvtHis2F.cc:27
void setHmc(TH2F *H2)
Definition: EvtHis2F.cc:45
bool AR(double xmass2, double ymass2)
Definition: EvtHis2F.cc:205
TH2F * getHwt()
Definition: EvtHis2F.cc:63
void HFill(double xmass2, double ymass2)
Definition: EvtHis2F.cc:129
void setHTitle(const char *htitle)
Definition: EvtHis2F.cc:40
void setZmax()
Definition: EvtHis2F.cc:177
double getZvalue(double m12_square, double m13_square)
Definition: EvtHis2F.cc:196
TH2F * getHdt()
Definition: EvtHis2F.cc:59
void setFile(const char *dtfile)
Definition: EvtHis2F.cc:35
void showFrame(TH2F *h2)
Definition: EvtHis2F.cc:83
static double Flat()
Definition: EvtRandom.cc:73