ATLAS Offline Software
Loading...
Searching...
No Matches
generate Namespace Reference

Classes

class  breit_generator
 breit-wigner generator More...
class  experiment
 given a histogram, get a generator for the distribution in that histogram, then run some pseudo-experiments and collect the mean and rms More...
class  generator_base
 base class for a "generator" class More...
class  hist_generator
 generate a distribution according to an input histogram (after smoothing) More...

Functions

void Zero (TH1D *hin)
double Entries (TH1D *hin)
double get_mean (const std::vector< double > &m)
 mean and rms
double get_rms (const std::vector< double > &m)
TH1D * integrate (TH1D *hin)
 generate an integrated distribution
void ScaleIntegral (TH1D *hin)
 scale a distribution by the integral
TH1D * PDF (TH1D *hin)
 generate a (normalised) pdf from a distribution
double Gradient (TH1D *hin)
int iGradient (TH1D *hin)
TH1D * smooth (TH1D *hin, bool sym)
 smooth a distribution about the maximum NB: maximum is not smoothed
double GetEntries (TH1D *h, int ilow, int ihi)
double GetEntries (TH1D *h)
void getRange (TH1D *s, int imax, double frac, int &lowerbin, int &upperbin, double &lowerfrac, double &upperfrac)
double findMean (TH1D *s, double frac=0.95)
int findMax (TH1D *s)
double rmsFrac (TH1D *s, double frac, double mean)
double rmsFrac (TH1D *s, double frac)
double rms95 (TH1D *s, double mean)
 get the fraction of the rms95 value with respect to the rms95 value of a gaussian
double rms95 (TH1D *s)

Variables

static const bool printPad = false
int Nevent_min = 0
int Nevent_max = 5000
const double frac = 0.95

Function Documentation

◆ Entries()

double generate::Entries ( TH1D * hin)

Definition at line 40 of file generate.cxx.

40 {
41 double N = 0;
42 for ( int i=1 ; i<=hin->GetNbinsX() ; i++ ) N += hin->GetBinContent(i);
43 return N;
44}

◆ findMax()

int generate::findMax ( TH1D * s)

find maximum bin

Definition at line 181 of file rmsFrac.cxx.

181 {
182 int imax = 1;
183 for ( int i=2 ; i<=s->GetNbinsX() ; i++ ) {
185 if ( s->GetBinContent(i)>s->GetBinContent(imax) ) imax = i;
186 }
187 return imax;
188}
int imax(int i, int j)

◆ findMean()

double generate::findMean ( TH1D * s,
double frac = 0.95 )

not enough entries for this fraction, ie we want 0.95 then need 5% to be 1 events, so need at least 20 events

maximum 20 iterations, calculating the mean of 95% until stable

Definition at line 84 of file rmsFrac.cxx.

84 {
85
86 const bool interpolate_flag = true;
87
88 // double entries = GetEntries(s,0,int(s->GetNbinsX())+1);
89 double entries = s->GetEffectiveEntries(); // s,0,int(s->GetNbinsX())+1);
90
91 // std::cout << "\nfindMean() " << s->GetName() << "\tentries: " << entries << std::endl;
92
96 if ( (1-frac)*entries<1 ) return 0;
97
98 s->GetXaxis()->SetRange(1,s->GetNbinsX());
99
100 double mean = s->GetMean();
101 double meane = s->GetMeanError();
102
103 int upperbin = 0;
104 int lowerbin = 0;
105
106 double upperfrac = 0;
107 double lowerfrac = 0;
108
109 int imax = 0;
110
111 int it=0;
112 constexpr int it_max=20;
113 constexpr int it_max_report=20;
115 for ( ; it<it_max ; it++ ) {
116
117
118 // std::cout << it << "\tmean " << mean << " +- " << meane << std::endl;
119
120 imax = s->GetXaxis()->FindBin(mean);
121
122 upperbin = imax;
123 lowerbin = imax;
124
125 upperfrac = 0;
126 lowerfrac = 0;
127
128 getRange( s, imax, frac, lowerbin, upperbin, lowerfrac, upperfrac );
129
130 s->GetXaxis()->SetRange(lowerbin, upperbin);
131
132 double m = s->GetMean();
133 double me = s->GetMeanError();
134
135 if ( it>0 ) {
136 if ( mean==m ||
137 std::fabs(mean-m)<me*1e-5 ||
138 std::fabs(mean-m)<meane*1e-5 ) {
139 mean = m;
140 meane = me;
141 // std::cout << "break after " << it << " iterations" << std::endl;
142 break;
143 }
144 }
145 //cppcheck-suppress oppositeInnerCondition
146 if ( it>=it_max_report ) {
147 std::cerr << s->GetName() << "\tMax iterations " << it << " reached" << std::endl;
148 }
149
150 mean = m;
151 meane = me;
152
153 }
154
155 // std::cout << s->GetName() << "\tmean " << mean << " +- " << meane << std::endl;
156
157 if ( interpolate_flag ) {
158
159 s->GetXaxis()->SetRange(lowerbin-1, upperbin+1);
160
161 double m = s->GetMean();
162 // double me = s->GetMeanError();
163
164 s->GetXaxis()->SetRange(1,s->GetNbinsX());
165
166 if ( upperfrac==lowerfrac ) return 0;
167
168 double inter_mean = mean + (frac-lowerfrac)*(m-mean)/(upperfrac-lowerfrac);
169
170 return inter_mean;
171 }
172 else {
173 s->GetXaxis()->SetRange(1,s->GetNbinsX());
174 return mean;
175 }
176
177}
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
double entries
Definition listroot.cxx:49
const double frac
Definition generate.cxx:295
void getRange(TH1D *s, int imax, double frac, int &lowerbin, int &upperbin, double &lowerfrac, double &upperfrac)
Definition rmsFrac.cxx:32

◆ get_mean()

double generate::get_mean ( const std::vector< double > & m)

mean and rms

Definition at line 47 of file generate.cxx.

47 {
48 double mean = 0;
49 for ( unsigned i=0 ; i<m.size() ; i++ ) mean += m[i];
50 return mean/m.size();
51}

◆ get_rms()

double generate::get_rms ( const std::vector< double > & m)

Definition at line 54 of file generate.cxx.

54 {
55 double mean = get_mean( m );
56 double rms = 0;
57 for ( unsigned i=0 ; i<m.size() ; i++ ) {
58 double r = (m[i]-mean);
59 rms += r*r;
60 }
61 return std::sqrt( rms/m.size() );
62}
int r
Definition globals.cxx:22
double get_mean(const std::vector< double > &m)
mean and rms
Definition generate.cxx:47

◆ GetEntries() [1/2]

double generate::GetEntries ( TH1D * h)

Definition at line 26 of file rmsFrac.cxx.

26 {
27 return GetEntries(h,1,h->GetNbinsX());
28}
Header file for AthHistogramAlgorithm.
double GetEntries(TH1D *h, int ilow, int ihi)
Definition rmsFrac.cxx:20

◆ GetEntries() [2/2]

double generate::GetEntries ( TH1D * h,
int ilow,
int ihi )

Definition at line 20 of file rmsFrac.cxx.

20 {
21 double sum = 0;
22 for ( int i=ilow ; i<=ihi ; i++ ) sum += h->GetBinContent(i);
23 return sum;
24}

◆ getRange()

void generate::getRange ( TH1D * s,
int imax,
double frac,
int & lowerbin,
int & upperbin,
double & lowerfrac,
double & upperfrac )

Definition at line 32 of file rmsFrac.cxx.

33 {
34
35 upperbin = imax;
36 lowerbin = imax;
37
38 upperfrac = 0;
39 lowerfrac = 0;
40
41 double entries = GetEntries(s,0,int(s->GetNbinsX())+1);
42
43 if ( entries>0 ) {
44
45 double sumn = s->GetBinContent(imax);
46 double tsum = s->GetBinContent(imax);
47
48 int i=1;
49 while ( true ) {
50
51 const int upperbin_i = imax+i;
52 const int lowerbin_i = imax-i;
53
54 if ( upperbin_i>s->GetNbinsX() || lowerbin_i<1 ) break;
55
56 tsum += s->GetBinContent(upperbin_i) + s->GetBinContent(lowerbin_i);
57
58 // std::cout << i << " frac: " << lowersum
59 // << "\tx " << s->GetBinCenter(lowerbin)
60 // << "\t " << s->GetBinCenter(upperbin)
61 // << "\ttsum " << tsum
62 // << "\tentries " << entries
63 // << std::endl;
64
65 if ( tsum>=entries*frac ) break;
66
67 sumn = tsum;
68
69 lowerfrac = sumn/entries;
70
71 upperbin = upperbin_i;
72 lowerbin = lowerbin_i;
73
74 i++;
75 }
76
77 upperfrac = tsum/entries;
78 }
79
80}

◆ Gradient()

double generate::Gradient ( TH1D * hin)

Definition at line 121 of file generate.cxx.

121 {
122 double gmax = 0;
123 double imax = 0;
124 for ( int i=2 ; i<hin->GetNbinsX() ; i++ ) {
125 double g = hin->GetBinContent(i+1)-hin->GetBinContent(i-1);
126 if ( std::fabs(gmax)<std::fabs(g) ) {
127 gmax = g;
128 imax = i;
129 }
130 }
131 return gmax/hin->GetBinError(imax);
132}

◆ iGradient()

int generate::iGradient ( TH1D * hin)

Definition at line 136 of file generate.cxx.

136 {
137 double gmax = 0;
138 double imax = 0;
139 for ( int i=2 ; i<hin->GetNbinsX() ; i++ ) {
140 double g = hin->GetBinContent(i+1)-hin->GetBinContent(i-1);
141 if ( std::fabs(gmax)<std::fabs(g) ) {
142 gmax = g;
143 imax = i;
144 }
145 }
146 return imax;
147}

◆ integrate()

TH1D * generate::integrate ( TH1D * hin)

generate an integrated distribution

Definition at line 68 of file generate.cxx.

68 {
69
70 TH1D* hint = (TH1D*)hin->Clone( (std::string(hin->GetName())+"-cumulative").c_str() );
71
72 Zero(hint);
73
74 hint->SetBinContent(0,0);
75 hint->SetBinContent(hin->GetNbinsX()+1,0);
76
77 for ( int i=1 ; i<=hin->GetNbinsX()+1 ; i++ ) {
78 double c = hin->GetBinContent(i) + hint->GetBinContent(i-1);
79 hint->SetBinContent(i, c);
80 hint->SetBinError(i, std::sqrt(c) );
81 }
82
83 hint->SetBinContent(hin->GetNbinsX()+1, hint->GetBinContent(hin->GetNbinsX()) );
84
85 return hint;
86}
void Zero(TH1D *hin)
Definition generate.cxx:32

◆ PDF()

TH1D * generate::PDF ( TH1D * hin)

generate a (normalised) pdf from a distribution

Definition at line 112 of file generate.cxx.

112 {
113 TH1D* h = integrate( hin );
114 h->SetDirectory(0);
115 ScaleIntegral( h );
116 return h;
117}
TH1D * integrate(TH1D *hin)
generate an integrated distribution
Definition generate.cxx:68
void ScaleIntegral(TH1D *hin)
scale a distribution by the integral
Definition generate.cxx:92

◆ rms95() [1/2]

double generate::rms95 ( TH1D * s)

Definition at line 293 of file rmsFrac.cxx.

293 {
294 double rms = rmsFrac(s, 0.95);
295 // printf("rms95 %12.10lf %12.10lf inv %12.10lf\n", rms, rms/0.8711151572, rms*1.1479538518 );
296 // rms *= 1.1479538518;
297 s->GetXaxis()->SetRange(1,s->GetNbinsX());
298 return rms;
299}
double rmsFrac(TH1D *s, double frac, double mean)
Definition rmsFrac.cxx:193

◆ rms95() [2/2]

double generate::rms95 ( TH1D * s,
double mean )

get the fraction of the rms95 value with respect to the rms95 value of a gaussian

Definition at line 285 of file rmsFrac.cxx.

285 {
286 double rms = rmsFrac(s, 0.95, mean);
287 // printf("rms95 %12.10lf %12.10lf inv %12.10lf\n", rms, rms/0.8711151572, rms*1.1479538518 );
288 // rms *= 1.1479538518;
289 s->GetXaxis()->SetRange(1,s->GetNbinsX());
290 return rms;
291}

◆ rmsFrac() [1/2]

double generate::rmsFrac ( TH1D * s,
double frac )

Definition at line 280 of file rmsFrac.cxx.

280 {
281 return rmsFrac(s, frac, findMean( s, frac ) );
282}
double findMean(TH1D *s, double frac=0.95)
Definition rmsFrac.cxx:84

◆ rmsFrac() [2/2]

double generate::rmsFrac ( TH1D * s,
double frac,
double mean )

not enough entries for this fraction, ie we want 0.95 then need 5% to be 1 events, so need at least 20 events

lazyness

technically not correct, since the rms is not definid for only one bin - really need NARROWER BINS!!

Definition at line 193 of file rmsFrac.cxx.

193 {
194
195 double rms = 0;
196 double erms = 0;
197
198 const bool interpolate_flag = true;
199
200 if ( s==0 ) { std::cerr << "rmsFrac() histogram s = " << s << std::endl; return 0; }
201 if ( s->GetXaxis()==0) { std::cerr << "rmsFrac() histogram s->GetXaxis() not definied for histogram " << s->GetName() << std::endl; return 0; }
202
203 // std::cout << "rmsFrac() " << s->GetName() << " " << GetEntries(s) << " " << GetEntries(s, 0, s->GetNbinsX()+1) << std::endl;
204
205 // double entries = GetEntries(s,0,int(s->GetNbinsX())+1);
206 double entries = s->GetEffectiveEntries(); // s,0,int(s->GetNbinsX())+1);
207
211 if ( (1-frac)*entries<1 ) return 0;
212
213 s->GetXaxis()->SetRange(1,s->GetNbinsX());
214
215 double m = mean;
216
217 int imax = s->GetXaxis()->FindBin(m);
218
219 // std::cout << "rmsFrac() mean " << m << " " << imax << " max " << s->GetBinCenter(findMax(s)) << " " << findMax(s) << std::endl;
220 // std::cout << "rmsFrac() entries: " << entries << "\t" << GetEntries(s) << "\t" << GetEntries(s,0,s->GetNbinsX()+1) << std::endl;
221
222 int upperbin = imax;
223 int lowerbin = imax;
224
225 double upperfrac = 0;
226 double lowerfrac = 0;
227
228 if ( entries>0 ) {
229 getRange( s, imax, frac, lowerbin, upperbin, lowerfrac, upperfrac );
230 }
231
232 // std::cout << "rmsFrac() " << s->GetName() << "\tlower bin " << lowerbin << "\t upper bin " << upperbin << std::endl;
233
234 // if ( upperbin!=lowerbin ) {
235 if ( true ) {
236
237 std::vector<double> stats;
238
239 // std::cout << "rmsFrac() GetEntries " << s->GetName() << " " << GetEntries(s) << std::endl;
240
241 if ( interpolate_flag ) { // && upperfrac!=lowerfrac ) {
242
243 double rlower = 0;
244 double erlower = 0;
245 s->GetXaxis()->SetRange(lowerbin, upperbin);
246 rlower = s->GetRMS();
247 // std::cout << "rmsFrac() GetEntries " << s->GetName() << " " << GetEntries(s,lowerbin,upperbin)/sentries << "\trlower " << rlower << std::endl;
248
249 double rupper = 0;
250 double erupper = 0;
251 s->GetXaxis()->SetRange(lowerbin-1, upperbin+1);
252 rupper = s->GetRMS();
253 // std::cout << "rmsFrac() GetEntries " << s->GetName() << " " << GetEntries(s,lowerbin-1,upperbin+1)/sentries << "\trupper " << rupper << std::endl;
254
255 // std::cout << "rmsFrac() Range " << s->GetBinLowEdge(lowerbin) << " - " << s->GetBinLowEdge(upperbin+1) << std::endl;
256
257 if ( upperfrac!=lowerfrac ) {
258 rms = rlower + (frac-lowerfrac)*(rupper-rlower)/(upperfrac-lowerfrac);
259 erms = erlower + (frac-lowerfrac)*(erupper-erlower)/(upperfrac-lowerfrac);
260 }
261 else {
262 rms = erms = 0;
263 }
264
265 }
266 else {
267 s->GetXaxis()->SetRange(lowerbin, upperbin);
268 rms = s->GetRMS();
269 }
270 }
271
272
273 s->GetXaxis()->SetRange(1,s->GetNbinsX());
274
275 return rms;
276
277}

◆ ScaleIntegral()

void generate::ScaleIntegral ( TH1D * hin)

scale a distribution by the integral

Definition at line 92 of file generate.cxx.

92 {
93
94 double N = hin->GetBinContent(hin->GetNbinsX());
95
96 for ( int i=0 ; i<=hin->GetNbinsX()+1 ; i++ ) {
97
98 double n = hin->GetBinContent(i);
99
100 double e = n/N;
101
102 double ee = std::sqrt(e*(1-e)/N);
103
104 hin->SetBinContent(i, e);
105 hin->SetBinError(i, ee);
106 }
107}

◆ smooth()

TH1D * generate::smooth ( TH1D * hin,
bool sym )

smooth a distribution about the maximum NB: maximum is not smoothed

a bit ad hoc, but seems to work

Definition at line 158 of file generate.cxx.

158 {
159
160 TH1D* hsmooth = (TH1D*)hin->Clone( (std::string(hin->GetName())+"-smooth").c_str() );
161 hsmooth->SetDirectory(0);
162
163 Zero( hsmooth );
164
165 // double _w[5] = { -3, 12, 17, 12, -3 };
166 // double w[7] = { -2, 3, 6, 7, 6, 3, -2 };
167 // double _w[11] = { -36, 9, 44, 69, 84, 89, 84, 69, 44, 9, -36 };
168 // double* w = _w+5;
169
170 double mg = 0;
171 int img = 0;
172
173 if ( sym ) {
174
175 TH1D* hint = integrate( hin ) ;
176 hint->SetDirectory(0);
177
178 for ( int i=2 ; i<hsmooth->GetNbinsX() ; i++ ) {
179 double gradient = std::fabs(0.5*(hint->GetBinContent(i+1) - hint->GetBinContent(i-1)));
180 if ( gradient>mg ) {
181 mg = gradient;
182 img = i;
183 }
184 }
185
186 delete hint;
187
188 }
189
190 // double sfac = std::sqrt(Entries(hin)/1000)*0.1;
191 const double sfac = 0.1;
192
193 // std::cout << "max gradient " << img << " " << mg << std::endl;
194
195 for ( int i=1 ; i<=hsmooth->GetNbinsX() ; i++ ) {
196
198 int Nbins = 2*sfac*std::fabs( i-img );
199
200 // Nbins = 0;
201
202 // std::cout << "bin " << i << "\tn " << Nbins << std::endl;
203
204 double sum = 0;
205 int wsum = 0;
206
207 for ( int j=i-Nbins ; j<=i+Nbins ; j++ ) {
208 if ( j<1 ) continue;
209 if ( j>hin->GetNbinsX() ) break;
210 sum += hin->GetBinContent(j);
211 wsum += 1;
212 }
213
214 if ( wsum>0 ) sum /= wsum;
215
216 // std::cout << i << " smoothed " << sum << " " << hin->GetBinContent(i) << std::endl;
217
218 hsmooth->SetBinContent(i, sum);
219 hsmooth->SetBinError(i, 0 );
220
221 }
222
223 return hsmooth;
224
225}

◆ Zero()

void generate::Zero ( TH1D * hin)

Definition at line 32 of file generate.cxx.

32 {
33 for ( int i=0 ; i<=hin->GetNbinsX()+1 ; i++ ) {
34 hin->SetBinContent(i, 0);
35 hin->SetBinError(i, 0);
36 }
37}

Variable Documentation

◆ frac

const double generate::frac = 0.95

Definition at line 295 of file generate.cxx.

◆ Nevent_max

int generate::Nevent_max = 5000

Definition at line 293 of file generate.cxx.

◆ Nevent_min

int generate::Nevent_min = 0

Definition at line 292 of file generate.cxx.

◆ printPad

const bool generate::printPad = false
static

Definition at line 30 of file generate.cxx.