ATLAS Offline Software
Classes | Functions | Variables
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 More...
 
double get_rms (const std::vector< double > &m)
 
TH1D * integrate (TH1D *hin)
 generate an integrated distribution More...
 
void ScaleIntegral (TH1D *hin)
 scale a distribution by the integral More...
 
TH1D * PDF (TH1D *hin)
 generate a (normalised) pdf from a distribution More...
 
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 More...
 
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 More...
 
double rms95 (TH1D *s)
 

Variables

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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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 }

◆ 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.

AllowedVariables::e
e
Definition: AsgElectronSelectorTool.cxx:37
beamspotman.r
def r
Definition: beamspotman.py:676
python.SystemOfUnits.s
int s
Definition: SystemOfUnits.py:131
mean
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="")
Definition: dependence.cxx:254
python.SystemOfUnits.m
int m
Definition: SystemOfUnits.py:91
python.SystemOfUnits.mg
int mg
Definition: SystemOfUnits.py:171
skel.it
it
Definition: skel.GENtoEVGEN.py:396
GetEntries
TGraphErrors * GetEntries(TH2F *histo)
Definition: TRTCalib_makeplots.cxx:4019
JetTiledMap::N
@ N
Definition: TiledEtaPhiMap.h:44
trigbs_dumpHLTContentInBS.stats
stats
Definition: trigbs_dumpHLTContentInBS.py:91
generate::get_mean
double get_mean(const std::vector< double > &m)
mean and rms
Definition: generate.cxx:47
generate::getRange
void getRange(TH1D *s, int imax, double frac, int &lowerbin, int &upperbin, double &lowerfrac, double &upperfrac)
Definition: rmsFrac.cxx:32
convertTimingResiduals.sum
sum
Definition: convertTimingResiduals.py:55
lumiFormat.i
int i
Definition: lumiFormat.py:85
beamspotman.n
n
Definition: beamspotman.py:731
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
extractSporadic.h
list h
Definition: extractSporadic.py:97
checkxAOD.frac
frac
Definition: Tools/PyUtils/bin/checkxAOD.py:257
imax
int imax(int i, int j)
Definition: TileLaserTimingTool.cxx:33
generate::ScaleIntegral
void ScaleIntegral(TH1D *hin)
scale a distribution by the integral
Definition: generate.cxx:92
h
generate::integrate
TH1D * integrate(TH1D *hin)
generate an integrated distribution
Definition: generate.cxx:68
beamspotnt.rms
rms
Definition: bin/beamspotnt.py:1266
entries
double entries
Definition: listroot.cxx:49
python.compressB64.c
def c
Definition: compressB64.py:93
generate::rmsFrac
double rmsFrac(TH1D *s, double frac)
Definition: rmsFrac.cxx:280
generate::Zero
void Zero(TH1D *hin)
Definition: generate.cxx:32
generate::findMean
double findMean(TH1D *s, double frac=0.95)
Definition: rmsFrac.cxx:84
Generated on Thu Nov 7 2024 22:14:42 for ATLAS Offline Software by   doxygen 1.8.18