-
Notifications
You must be signed in to change notification settings - Fork 0
/
Interpolate.h
212 lines (195 loc) · 8.64 KB
/
Interpolate.h
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
#include "TAxis.h"
#include "TH1.h"
//#include "TMath.h"
#include <TH2.h>
//#include <TStyle.h>
//#include <TCanvas.h>
//#include "TFile.h"
//#include "TTree.h"
//#include "TVirtualFitter.h"
//#include "TLegend.h"
#include <iostream>
//#include <sstream>
//#include "TString.h"
//#include <fstream>
//#include "TF1.h"
//#include "TLine.h"
//#include "quicksave.C"
//#include "CMSStyle.C"
#include "utilities.h"
//#include "util.C"
//#include "params_arg.h"
//#include "params_arg.h"
//#include<string.h>
//#include<map>
//#include <algorithm> // std::sort
//#include <vector>
using namespace std;
float bilinear_interpolate(float x1,float val1,float x2, float val2, float x_target);//wad
float bilinear_interpolate( float x1, float x2, float y1, float y2,
float val11, float val12, float val21, float va22,
float x_target,float y_target);
void interpolate_collumn(TH2F* hin, int col, int rowstart, int rowend); //wad
void interpolate_uneven_column(TH2F* hin, int col, int rowstart, int rowend);//wad
void interpolate_collumn_rowwise(TH2F* hin, int col, int rowstart, int rowend); //wad
TH2F* HistOr(TH2F* h1, TH2F* h2, float blankval, float tollorence);
void interpolate_Zywicki(TH2F* shittyhist,bool bigcol = true);
void interpolate_Zywicki(TH2D* shittyhist,bool bigcol = true);
///////////TH2D///////////////
void interpolate_collumn(TH2D* hin, int col, int rowstart, int rowend); //wad
void interpolate_uneven_column(TH2D* hin, int col, int rowstart, int rowend);//wad
void interpolate_collumn_rowwise(TH2D* hin, int col, int rowstart, int rowend); //wad
TH2D* HistOr(TH2D* h1, TH2D* h2, float blankval, float tollorence);
// ##############################
float bilinear_interpolate(float x1,float val1,float x2, float val2, float x_target){
return (x2-x_target)*val1/(x2-x1) + (x_target-x1)*val2/(x2-x1);
}
float bilinear_interpolate( float x1, float x2, float y1, float y2,
float val11, float val12, float val21, float val22,
float x_target,float y_target){
//take four points on the corners of a rectangle, bilinearly interpolate inside the rectangle.
return ((x2-x_target)*(y2-y_target)*val11 + (x_target-x1)*(y2-y_target)*val21 + (x2-x_target)*(y_target-y1)*val12 + (x_target-x1)*(y_target-y1)*val22)/((x2-x1)*(y2-y1));
}
/*float bilinear_interpolate_center(float val11, float val12, float val21, float va22){
//take four points on the edges of a rectangle and interpolate the point in the middle.
//this is just the average of the four points.
return 0.25*(val11+val12+val21+val22);
}*/
void interpolate_collumn(TH2F* hin, int col, int rowstart, int rowend){
//assumes that alternating cells in the are occupied.
//go up the collumn, interpolating the blanks
TAxis* y = hin->GetYaxis();
for(int irow = rowstart; irow <= rowend; irow++){
if(aeq(hin->GetBinContent(col,irow),0)){ //fixed
hin->SetBinContent(col,irow, //fixed
bilinear_interpolate( y->GetBinCenter(irow-1),hin->GetBinContent(col,irow-1), y->GetBinCenter(irow+1), hin->GetBinContent(col,irow+1), y->GetBinCenter(irow) )); //fixed
}
}
}
void interpolate_uneven_column(TH2F* hin, int col, int rowstart, int rowend){
TAxis* y = hin->GetYaxis();
//first we crawl up the column, looking for the first non-zero entry
//then we crawl up to find the next non-zero entry.
//then you fill in all the rows inbetween and
int lowrow = rowstart;
int hirow = rowstart+1;
do{
if(aeq(hin->GetBinContent(col,lowrow),0)){ //fixed
lowrow++;
continue;
}
//now on a non-zero entry
hirow = lowrow+1;
while(aeq(hin->GetBinContent(col,hirow),0) && hirow <= rowend) hirow++;//fixed
if(hirow > rowend) break;
//hirow is now on a no-zero entry.
for(int irow = lowrow+1; irow<hirow; irow++){
hin->SetBinContent(col, irow,//fixed
bilinear_interpolate( y->GetBinCenter(lowrow),hin->GetBinContent(col,lowrow), y->GetBinCenter(hirow), hin->GetBinContent(col,hirow), y->GetBinCenter(irow) ));
}
lowrow = hirow;
}while(lowrow<rowend);
}
void interpolate_collumn_rowwise(TH2F* hin, int col, int rowstart, int rowend){
//assuming a blank column with filled columns on either side.
//interpolate from neighbors in rows.
TAxis* x = hin->GetYaxis();
for(int irow = rowstart; irow <= rowend; irow++){
hin->SetBinContent(col, irow,
bilinear_interpolate( x->GetBinCenter(col-1),hin->GetBinContent(col-1,irow), x->GetBinCenter(col+1), hin->GetBinContent(col+1,irow), x->GetBinCenter(col) ));
}
}
TH2F* HistOr(TH2F* h1, TH2F* h2, float blankval, float tollorence){
//takes two partly interpolated histograms and takes the Or of them.
//if they both have a value, give preference to h1
//if they're not of the same dimension, return h1.
if(h1->GetXaxis()->GetNbins() != h2->GetXaxis()->GetNbins() ||
h1->GetYaxis()->GetNbins() != h2->GetYaxis()->GetNbins() ) return h1;
TH2F* hout = (TH2F*) h1->Clone("hout");
for(int i=0;i<=h1->GetXaxis()->GetNbins();i++){
for(int j=0;j<=h1->GetYaxis()->GetNbins();j++){
if( fabs(h1->GetBinContent(i,j) - blankval) < tollorence) hout->SetBinContent(i,j,h2->GetBinContent(i,j));
}
}
return h1;
}
void interpolate_collumn(TH2D* hin, int col, int rowstart, int rowend){
//assumes that alternating cells in the are occupied.
//go up the collumn, interpolating the blanks
TAxis* y = hin->GetYaxis();
for(int irow = rowstart; irow <= rowend; irow++){
if(aeq(hin->GetBinContent(col,irow),0)){ //fixed
hin->SetBinContent(col,irow, //fixed
bilinear_interpolate( y->GetBinCenter(irow-1),hin->GetBinContent(col,irow-1), y->GetBinCenter(irow+1), hin->GetBinContent(col,irow+1), y->GetBinCenter(irow) )); //fixed
}
}
}
void interpolate_uneven_column(TH2D* hin, int col, int rowstart, int rowend){
TAxis* y = hin->GetYaxis();
//first we crawl up the column, looking for the first non-zero entry
//then we crawl up to find the next non-zero entry.
//then you fill in all the rows inbetween and
int lowrow = rowstart;
int hirow = rowstart+1;
do{
if(aeq(hin->GetBinContent(col,lowrow),0)){ //fixed
lowrow++;
continue;
}
//now on a non-zero entry
hirow = lowrow+1;
while(aeq(hin->GetBinContent(col,hirow),0) && hirow <= rowend) hirow++;//fixed
if(hirow > rowend) break;
//hirow is now on a no-zero entry.
for(int irow = lowrow+1; irow<hirow; irow++){
hin->SetBinContent(col, irow,//fixed
bilinear_interpolate( y->GetBinCenter(lowrow),hin->GetBinContent(col,lowrow), y->GetBinCenter(hirow), hin->GetBinContent(col,hirow), y->GetBinCenter(irow) ));
}
lowrow = hirow;
}while(lowrow<rowend);
}
void interpolate_collumn_rowwise(TH2D* hin, int col, int rowstart, int rowend){
//assuming a blank column with filled columns on either side.
//interpolate from neighbors in rows.
TAxis* x = hin->GetYaxis();
for(int irow = rowstart; irow <= rowend; irow++){
hin->SetBinContent(col, irow,
bilinear_interpolate( x->GetBinCenter(col-1),hin->GetBinContent(col-1,irow), x->GetBinCenter(col+1), hin->GetBinContent(col+1,irow), x->GetBinCenter(col) ));
}
}
TH2D* HistOr(TH2D* h1, TH2D* h2, float blankval, float tollorence){
//takes two partly interpolated histograms and takes the Or of them.
//if they both have a value, give preference to h1
//if they're not of the same dimension, return h1.
if(h1->GetXaxis()->GetNbins() != h2->GetXaxis()->GetNbins() ||
h1->GetYaxis()->GetNbins() != h2->GetYaxis()->GetNbins() ) return h1;
TH2D* hout = (TH2D*) h1->Clone("hout");
for(int i=0;i<=h1->GetXaxis()->GetNbins();i++){
for(int j=0;j<=h1->GetYaxis()->GetNbins();j++){
if( fabs(h1->GetBinContent(i,j) - blankval) < tollorence) hout->SetBinContent(i,j,h2->GetBinContent(i,j));
}
}
return h1;
}
void interpolate_Zywicki(TH2D* shittyhist,bool bigcol){
interpolate_Zywicki((TH2F* )shittyhist,bigcol);
}
void interpolate_Zywicki(TH2F* shittyhist,bool bigcol){
//This assumes a histogram with the following binning:
//float mSBins[nS+1] = {187.5, 212.5, 237.5 ,262.5 ,287.5, 312.5, 337.5, 362.5, 375, 425, 475, 525}
//float mHBins[nH+1] = {130, 140, 160, 190, 210, 220, 230, 235,245, 255, 260, 270, 280, 285,295, 305, 310,320, 330};
//with points filled in according to the spotty-ass grid that PZ sent us.
//very spicific interpolation of
interpolate_uneven_column(shittyhist, 2, 2, 5);
interpolate_uneven_column(shittyhist, 3, 1, 8);
interpolate_uneven_column(shittyhist, 4, 1, 11);
interpolate_uneven_column(shittyhist, 5, 1, 14);
interpolate_uneven_column(shittyhist, 6, 1, 17);
interpolate_uneven_column(shittyhist, 7, 1, 16);
//skip 8
interpolate_uneven_column(shittyhist, 9, 1, 19);
interpolate_uneven_column(shittyhist, 10, 2, 19); //these are regular, so use a simpler algor.
interpolate_uneven_column(shittyhist, 11, 2, 19);
//now do 8
interpolate_collumn_rowwise(shittyhist, 8, 1, 15);
}