72 {
73
76 int* dp;
77 int nd1, nd2, ii, sn;
78
80
81 if ( sn == 2 )
82 {
83 nd1 = 0;
84 nd2 = 1;
85 }
86
87 if ( sn == 0 )
88 {
89 nd1 = 1;
90 nd2 = 2;
91 }
92
93 if ( sn == 1 )
94 {
95 nd1 = 0;
96 nd2 = 2;
97 }
100
102 int d1 = dp[0];
103 int d2 = dp[1];
104 int d3 = dp[2];
105 int d4 = dp[3];
106
108 TH1F* did1 = (TH1F*)
f.Get( DA1 );
109 TH1F* did2 = (TH1F*)
f.Get( DA2 );
110 TH2F* hid = (TH2F*)
f.Get( hp );
111
112 TAxis* d1x = did1->GetXaxis();
113 TAxis* d2x = did2->GetXaxis();
114 TAxis* xaxis = hid->GetXaxis();
115 TAxis* yaxis = hid->GetYaxis();
116
117 int BINSd1 = did1->GetXaxis()->GetLast();
118 int BINSd2 = did2->GetXaxis()->GetLast();
119 int BINSx = xaxis->GetLast();
120 int BINSy = yaxis->GetLast();
121 int BINS = BINSx * BINSy;
122
123 double av1, av2, avm1, avm2;
124 avm1 = 0.;
125 avm2 = 0.;
126 double yvalue, ymax = 0.0;
127 int i, j, binxy;
128
129
130 for ( i = 1; i < BINSd1 + 1; i++ )
131 {
132 av1 = did1->GetBinContent( i );
133 if ( av1 > avm1 ) avm1 = av1;
134 }
135
136
137 for ( i = 1; i < BINSd2 + 1; i++ )
138 {
139 av2 = did2->GetBinContent( i );
140 if ( av2 > avm2 ) avm2 = av2;
141 }
142
143
144
145 for ( i = 1; i < BINSx + 1; i++ )
146 {
147 for ( j = 1; j < BINSy + 1; j++ )
148 {
149 binxy = hid->GetBin( i, j );
150 yvalue = hid->GetBinContent( binxy );
151
152 if ( yvalue > ymax ) ymax = yvalue;
153 }
154 }
155
156loop:
158
159 EvtParticle *id1, *id2, *id3, *id4, *s1;
160 EvtVector4R pd1, pd2, pd3, pd4, ps;
161 EvtVector4R dp1, dp2;
162 double xmass2, ymass2;
163
168
173
176
177 xmass2 = ( pd1 + pd2 ).mass2();
178 ymass2 = ( pd3 + pd4 ).mass2();
179
180 int xbin = hid->GetXaxis()->FindBin( xmass2 );
181 int ybin = hid->GetYaxis()->FindBin( ymass2 );
182 int xybin = hid->GetBin( xbin, ybin );
183 double zvalue = hid->GetBinContent( xybin );
184 double xratio = zvalue / ymax;
185 if ( xratio == 0 ) goto loop;
187 if ( rd1 > xratio ) goto loop;
188
189 double da1 = dp1.
get( 3 ) / dp1.
d3mag();
190 double da2 = dp2.
get( 3 ) / dp2.
d3mag();
191
192 int dbin1 = did1->FindBin( da1 );
193 int dbin2 = did2->FindBin( da2 );
194
195 double dr1 = ( did1->GetBinContent( dbin1 ) ) / avm1;
196 double dr2 = ( did2->GetBinContent( dbin2 ) ) / avm2;
197 if ( dr1 == 0 || dr2 == 0 ) goto loop;
199 if ( rd1 > dr1 ) goto loop;
200
202 if ( rd1 > dr2 ) goto loop;
203
204 return;
205}
TFile f("ana_bhabha660a_dqa_mcPat_zy_old.root")
const char * setFileName()
const char * setDaugAng(int i)
EvtParticle * getDaug(int i)
double initializePhaseSpace(int numdaughter, EvtId *daughters, double poleSize=-1., int whichTwo1=0, int whichTwo2=1)