//----------Author: Markus D. Oldenburg
//----------Last Modified: 05.03.1999
//----------Copyright: © MDO Production 1997
#include "MWavelet.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TAxis.h"
#include "TMath.h"
#include "TProfile.h"
#include "TObject.h"
////////////////////////////////////////////////////////////////////////////
// MWavelet //
// //
// The MWavelet class calculates factorial wavelet moments of 2nd, 3rd //
// and 4th order and for one, two and three dimensional distributions. //
// It is also able to display the results in a fancy way. It can do the //
// inverse wavelet transformation and some histogramming of the facorial //
// wavelet moments of a whole set of distributions //
// The class design itself isn't very good yet - actually it is no good //
// at all!
//
// 20.05.1999 I had to remove all calculations for dimension = 3 due to //
// problems with cint. The code is still there but it is //
// commented out. Also I commented out the concerning //
// variabels. //
// In the constructor the setting of all pointers to NULL was //
// commented out as well to avoid warnings like //
// "unused variable". //
////////////////////////////////////////////////////////////////////////////
ClassImp(MWavelet)
MWavelet::MWavelet(Char_t *na) {
// MWavelet constructor. Gives the new MWavelet object its name
// and sets several pointers to NULL.
name = na;
/*
TH1D *wav_mom_ord2 = NULL;
TH1D *wav_mom_ord3 = NULL;
TH1D *wav_mom_ord4 = NULL;
TH1D *coeff = NULL;
TH1D *coeff_scaled = NULL;
TH1D *revert = NULL;
TProfile *prof_ord2 = NULL;
TProfile *prof_ord3 = NULL;
TProfile *prof_ord4 = NULL;
TCanvas *canvas = NULL;
TPad *original = NULL;
TPad *wav2 = NULL;
TPad *wav3 = NULL;
TPad *wav4 = NULL;
TNtuple *ord2 = NULL;
TNtuple *ord3 = NULL;
TNtuple *ord4 = NULL;
*/
}
MWavelet::~MWavelet() {
// MWavelet destructor. Does nothing at all.
}
void MWavelet::CalcHaar(TH1 *input_histo, Int_t jm, Int_t ord_input, Int_t option) {
// Calculates the factorial wavelet moments of the given
// input distribution (input_histo) up to the scale j (jm)
// and for the given order ord_input. By option=1 is specified if
// the output is drawn immediatly.
Int_t j, jend;
Int_t i1, i1end, i2 /*, i3*/;
Int_t l1, l1end, l2 /*, l3*/;
Int_t k1, k1end, k2 /*, k3*/;
// Int_t x, y, z;
Int_t ord, ord_end = 0, min_maxbin;
Int_t dim;
Int_t bin_x, bin_y/*, bin_z */;
Int_t maxbin;
Int_t maxbin_histo_x, maxbin_histo_y/*, maxbin_histo_z */;
Double_t eintrag1, eintrag2;
// Double_t phi[2], psi1d[2];
for (j=0;j<10; j++) {
max_coeff[j] = 0.;
max_coeff_scaled[j] = 0;
}
if (!coeff) {
TH1D *a = new TH1D("coeff", "Wavelet Koeffizienten", 512, 0., 512.);
TH1D *b = new TH1D("coeff_scaled", "Wavelet Koeffizienten (skaliert)", 512, 0., 512.);
coeff = a;
coeff_scaled = b;
}
else {
coeff->Reset();
coeff_scaled->Reset();
}
// Evaluation of the dimension of the histogram
dim = input_histo->GetDimension();
// Evaluation of the order of the histogram
if (ord_input < 0 || ord_input == 1 || ord_input > 4) {
printf("Only wavelet moments of range 2 to 4 supported!n");
return;
}
// Dimension = 1 -----
if (dim==1) {
// Set jmax. If it is not set, calculate maximal possible value.
maxbin_histo_x = input_histo->GetNbinsX();
if (jm == 0) {
while (((Int_t) pow(2, jm))<=maxbin_histo_x) jm++;
jm--;
}
jmax = jm;
maxbin = (Int_t) pow(2, jmax);
// Restrict jmax to 512 (which equals 9 orders).
if (maxbin > 512) {
maxbin = 512;
jmax = 9;
}
if (maxbin_histo_x < maxbin) {
printf("pow(2, jmax = %d) > #bins of histogram!n", jmax);
return;
}
else if (maxbin_histo_x > maxbin) {
printf("Take into account only the first %d bins!n", maxbin);
}
for (bin_x=1; bin_x<=maxbin; bin_x++) {
matrix_1dim[jmax][bin_x-1] = input_histo->GetBinContent(bin_x);
}
// Allocate psi.
psi_1dim[0] = 1.;
psi_1dim[1] = -1.;
// Loop over different orders.
if (ord_input==0) {
ord_end = 4;
}
else {
ord_end=ord_input;
}
for (ord=2; ord<=ord_end; ord++) {
// Set correlations rho to zero.
for (j=0; j<jmax; j++) {
i1end = (Int_t) pow(2, j);
for (i1=0; i1<i1end; i1++) {
rho_1dim[j][i1] = 0.;
}
}
// Loop over different scales.
for (j=jmax; j>0; j--) {
this->Mitteln_1dim(matrix_1dim, j);
l1end = (Int_t) pow(2, j);
for (l1=0; l1<l1end; l1+=2) {
for (i1=0; i1<=1; i1++) {
ma_1dim[i1] = matrix_1dim[j][i1 + l1];
}
if (ord==2) {
eintrag1 = (matrix_1dim[j][l1]-matrix_1dim[j][l1+1])/
pow(2,jmax-j+1.);
eintrag2 = (matrix_1dim[j][l1]-matrix_1dim[j][l1+1])/2.;
if (TMath::Abs(eintrag1) > max_coeff[j])
max_coeff[j] = TMath::Abs(eintrag1);
if (TMath::Abs(eintrag2) > max_coeff_scaled[j])
max_coeff_scaled[j] = TMath::Abs(eintrag2);
coeff->SetBinContent(Int_t(coeff->GetBin(Int_t(pow(2., Double_t(j)-1.) + Double_t(l1/2.+1.)))), eintrag1);
coeff_scaled->SetBinContent(Int_t(coeff_scaled->GetBin(Int_t(pow(2., Double_t(j)-1.) + Double_t(l1/2.+1.)))), eintrag2);
}
wertges=0.;
if (ord==2) {
this->Ord2_1dim(ma_1dim, psi_1dim);
}
else if (ord==3) {
this->Ord3_1dim(ma_1dim, psi_1dim);
}
else if (ord==4) {
this->Ord4_1dim(ma_1dim, psi_1dim);
}
rho_1dim[j-1][l1/2] += wertges;
eintrag1 = (matrix_1dim[1][0]+matrix_1dim[1][1])/
pow(2., jmax);
eintrag2 = (matrix_1dim[1][0]+matrix_1dim[1][1])/2.;
coeff->SetBinContent(1, eintrag1);
coeff_scaled->SetBinContent(1, eintrag2);
max_coeff[0] = eintrag1;
max_coeff_scaled[0] = eintrag2;
}
}
// Calculation of wavelet moments.
for (j=0; j<jmax; j++) {
mo[ord][j] = 0.;
k1end = (Int_t) pow(2, j);
for (k1=0; k1<k1end; k1++) {
mo[ord][j] += (rho_1dim[j][k1] / k1end);
}
}
}
}
// Dimesion = 2 --------
else if (dim==2) {
// Set jmax. If it is not set, calculate maximal possible value.
maxbin_histo_x = input_histo->GetNbinsX();
maxbin_histo_y = input_histo->GetNbinsY();
min_maxbin = (maxbin_histo_x < maxbin_histo_y)
? maxbin_histo_x : maxbin_histo_y;
if (jm == 0) {
while (((Int_t) pow(2, jm))<=min_maxbin) jm++;
jm--;
}
jmax = jm;
maxbin = (Int_t) pow(2, jmax);
// Restrict jmax to 64 (which equals 6 orders).
if (maxbin > 64) {
maxbin = 64;
jmax = 6;
}
if (maxbin_histo_x < maxbin || maxbin_histo_y < maxbin) {
printf("pow(2, jmax = %d) > #bins of histogram in one direction!n", jmax);
return;
}
if (maxbin_histo_x > maxbin) {
printf("Take into account only the first %d bins of x!n", maxbin);
}
if (maxbin_histo_y > maxbin) {
printf("Take into account only the first %d bins of y!n", maxbin);
}
for (bin_x=1; bin_x<=maxbin; bin_x++) {
for (bin_y=1; bin_y<=maxbin; bin_y++) {
matrix_2dim[jmax][bin_x-1][bin_y-1] = input_histo->GetBinContent(input_histo->GetBin(bin_x, bin_y));
}
}
// Allocate psi.
psi_2dim[0][0][0] = 1.;
psi_2dim[0][0][1] = 1.;
psi_2dim[0][1][0] = -1.;
psi_2dim[0][1][1] = -1.;
psi_2dim[1][0][0] = 1.;
psi_2dim[1][0][1] = -1.;
psi_2dim[1][1][0] = 1.;
psi_2dim[1][1][1] = -1.;
psi_2dim[2][0][0] = 1.;
psi_2dim[2][0][1] = -1.;
psi_2dim[2][1][0] = -1.;
psi_2dim[2][1][1] = 1.;
// Loop over different orders.
if (ord_input==0) ord_end = 4;
else ord_end=ord_input;
for (ord=2; ord<=ord_end; ord++) {
// Set correlations rho to zero.
for (j=0; j<jmax; j++) {
jend = (Int_t) pow(2 ,j);
for (i1=0; i1<jend; i1++) {
for (i2=0; i2<jend; i2++) {
rho_2dim[j][i1][i2] = 0.;
}
}
}
// Loop over different scales.
for (j=jmax; j>0; j--) {
this->Mitteln_2dim(matrix_2dim, j);
l1end = (Int_t) pow(2, jmax)-1;
for (l1=0; l1<l1end; l1+=2) {
for (l2=0; l2<l1end; l2+=2) {
for (i1=0; i1<=1; i1++) {
for (i2=0; i2<=1; i2++) {
ma_2dim[i1][i2] = matrix_2dim[j][i1+l1][i2+l2];
}
}
wertges=0.;
if (ord==2) {
this->Ord2_2dim(ma_2dim, psi_2dim);
}
else if (ord==3) {
this->Ord3_2dim(ma_2dim, psi_2dim);
}
else if (ord==4) {
this->Ord4_2dim(ma_2dim, psi_2dim);
}
rho_2dim[j-1][l1/2][l2/2] += wertges;
}
}
}
// Calculation of wavelet moments.
for (j=0; j<jmax; j++) {
mo[ord][j] = 0.;
k1end = (Int_t) pow(2,j);
for (k1=0; k1<k1end; k1++) {
for (k2=0; k2<k1end; k2++) {
mo[ord][j] += (rho_2dim[j][k1][k2] / (3 * pow(2, 2*j)));
}
}
}
}
}
/*
// Dimension = 3 ------
else if (dim==3) {
// Set jmax. If it is not set, calculate maximal possible value.
maxbin_histo_x = input_histo->GetNbinsX();
maxbin_histo_y = input_histo->GetNbinsY();
maxbin_histo_z = input_histo->GetNbinsZ();
min_maxbin = (maxbin_histo_x < maxbin_histo_y) ? maxbin_histo_x : maxbin_histo_y;
min_maxbin = (min_maxbin < maxbin_histo_z) ? min_maxbin : maxbin_histo_z;
if (jm == 0) {
while (((Int_t) pow(2, jmax))<=min_maxbin) jm++;
jm--;
}
jmax = jm;
maxbin = (Int_t) pow(2, jmax);
// Restrict jmax to 64 (which equals 6 orders).
if (maxbin > 64 ) {
maxbin = 64;
jmax = 6;
}
if (maxbin_histo_x < maxbin || maxbin_histo_y < maxbin || maxbin_histo_z < maxbin) {
printf("pow(2, jmax = %d) > #bins of histogram in one direction!n", jmax);
return;
}
if (maxbin_histo_x > maxbin) {
printf("Take into account only the first %d bins of x!n", maxbin);
}
if (maxbin_histo_y > maxbin) {
printf("Take into account only the first %d bins of y!n", maxbin);
}
if (maxbin_histo_z > maxbin) {
printf("Take into account only the first %d bins of z!n", maxbin);
}
for (bin_x=1; bin_x<=maxbin; bin_x++) {
for (bin_y=1; bin_y<=maxbin; bin_y++) {
for (bin_z=1; bin_z<=maxbin; bin_z++) {
matrix_3dim[jmax][bin_x-1][bin_y-1][bin_z-1] = input_histo->GetBinContent(input_histo->GetBin(bin_x, bin_y, bin_z));
}
}
}
// Allocate psi.
phi[0] = 1;
phi[1] = 1;
for (x=0; x<=1; x++) {
for (y=0; y<=1; y++) {
for (z=0; z<=1; z++) {
psi_3dim[0][x][y][z] = phi[x] * phi[y] * psi1d[z];
psi_3dim[1][x][y][z] = phi[x] * psi1d[y] * phi[z];
psi_3dim[2][x][y][z] = psi1d[x] * phi[y] * phi[z];
psi_3dim[3][x][y][z] = phi[x] * psi1d[y] * psi1d[z];
psi_3dim[4][x][y][z] = psi1d[x] * phi[y] * psi1d[z];
psi_3dim[5][x][y][z] = psi1d[x] * psi1d[y] * phi[z];
psi_3dim[6][x][y][z] = psi1d[x] * psi1d[y] * psi1d[z];
}
}
}
// Loop over different orders.
if (ord_input==0) ord_end = 4;
else ord_end=ord_input;
for (ord=2; ord<=ord_end; ord++) {
// Set correlations rho to zero.
for (j=0; j<jmax; j++) {
i1end = (Int_t) pow(2, j);
for (i1=0; i1<i1end; i1++) {
for (i2=0; i2<i1end; i2++) {
for (i3=0; i3<i1end; i3++) {
rho_3dim[j][i1][i1][i3] = 0.;
}
}
}
}
// Loop over different scales.
for (j=jmax; j>=1; j--) {
this->Mitteln_3dim(matrix_3dim, j);
l1end = (Int_t) pow(2, j) - 1;
for (l1=0; l1<l1end; l1+=2) {
for (l2=0; l2<l1end; l2+=2) {
for (l3=0; l3<l1end; l3+=2) {
for (i1=0; i1<=1; i1++) {
for (i2=0; i2<=1; i2++) {
for (i3=0; i3<=1; i3++) {
ma_3dim[i1][i2][i3] = matrix_3dim[j][i1+l1][i2+l2][i3+l3];
}
}
}
wertges=0.;
if (ord==2) {
this->Ord2_3dim(ma_3dim, psi_3dim);
}
else if (ord==3) {
this->Ord3_3dim(ma_3dim, psi_3dim);
}
else if (ord==4) {
this->Ord4_3dim(ma_3dim, psi_3dim);
}
rho_3dim[j-1][l1/2][l2/2][l3/2] += wertges;
}
}
}
}
// Calculation of wavelet moments.
for (j=0; j<jmax; j++) {
mo[ord][j] = 0.;
k1end = (Int_t) pow(2, j);
for (k1=0; k1<k1end; k1++) {
for (k2=0; k2<k1end; k2++) {
for (k3=0; k3<k1end; k3++) {
mo[ord][j] += (rho_3dim[j][k1][k2][k3] / (7 * pow(2, 3*j)));
}
}
}
}
}
}
*/
// Generates output.
for (ord=2; ord<=ord_end; ord++) {
if (ord==2) {
if(wav_mom_ord2) {
delete wav_mom_ord2;
}
sprintf(histo, "%s__wav_mom_ord2", name);
wav_mom_ord2 = new TH1D(histo, "Wavelet Momente 2. Ordnung", jmax, 0., jmax);
for (j=0; j<jmax; j++) {
wav_mom_ord2->SetBinContent(j+1, (Stat_t)mo[ord][j]);
}
}
if (ord==3) {
if(wav_mom_ord3) {
delete wav_mom_ord3;
}
sprintf(histo, "%s__wav_mom_ord3", name);
wav_mom_ord3 = new TH1D(histo, "Wavelet Momente 3. Ordnung", jmax, 0., jmax);
for (j=0; j<jmax; j++) {
wav_mom_ord3->SetBinContent(j+1, (Stat_t)mo[ord][j]);
}
}
if (ord==4) {
if(wav_mom_ord4) {
delete wav_mom_ord4;
}
sprintf(histo, "%s__wav_mom_ord4", name);
wav_mom_ord4 = new TH1D(histo, "Wavelet Momente 4. Ordnung", jmax, 0., jmax);
for (j=0; j<jmax; j++) {
wav_mom_ord4->SetBinContent(j+1, (Stat_t)mo[ord][j]);
}
}
}
if (option == 1 && ord_input == 0) {
this->Draw(input_histo, option);
}
return;
}
Double_t MWavelet::GetWavMom(Int_t i, Int_t ord) {
// Returns the factorial wavelet moment of scale i and order ord.
return mo[ord][i];
}
Double_t MWavelet::GetMaxCoeff(Int_t ord) {
// Returns the maximal coefficient of the factorial wavelet transform
// for a given order ord.
return max_coeff[ord];
}
Double_t MWavelet::GetMaxCoeffSc(Int_t ord) {
// Returns the scaled maximal wavelet coefficients of the wavelet transform
// for a given order ord.
return max_coeff_scaled[ord];
}
void MWavelet::Mitteln_1dim(Double_t matrix_1dim[10][512], Int_t j) {
// Calculates the mean of the given 1-dimensional distribution of
// wavelet coefficients of scale j.
Int_t k1, m1, end;
Double_t mamittel;
end = (Int_t) pow(2, j-1);
for (k1=0; k1<end; k1++) {
mamittel=0;
for (m1=0; m1<=1; m1++) {
mamittel += matrix_1dim[j][m1+2*k1];
}
matrix_1dim[j-1][k1] = mamittel;
}
return;
}
void MWavelet::Mitteln_2dim(Double_t matrix[7][64][64], Int_t j) {
// Calculates the mean of the given 2-dimensional distribution of
// wavelet coefficients of scale j.
Int_t k1, k2, m1, m2, end;
Double_t mamittel;
end = (Int_t) pow(2, j-1);
for (k1=0; k1<end; k1++) {
for (k2=0; k2<end; k2++) {
mamittel = 0.;
for (m1=0; m1<=1; m1++) {
for (m2=0; m2<=1; m2++) {
mamittel += matrix_2dim[j][m1+2*k1][m2+2*k2];
}
}
matrix[j-1][k1][k2] = mamittel;
}
}
return;
}
/*
void MWavelet::Mitteln_3dim(Double_t matrix[7][64][64][64], Int_t j) {
// Calculates the mean of the given 3-dimensional distribution of
// wavelet coefficients of scale j.
Int_t k1, k1end, k2, k3, m1, m2, m3;
Double_t mamittel;
k1end = (Int_t) pow(2, j-1);
for (k1=0; k1<k1end; k1++) {
for (k2=0; k2<k1end; k2++) {
for (k3=0; k3<k1end; k3++) {
mamittel = 0.;
for (m1=0; m1<=1; m1++) {
for (m2=0; m2<=1; m2++) {
for (m3=0; m3<=1; m3++) {
mamittel += matrix_3dim[j][m1+2*k1][m2+2*k2][m3+2*k3];
}
}
}
}
}
}
return;
}
*/
void MWavelet::Ord2_1dim(Double_t ma_1dim[2], Double_t psi_1dim[2]) {
// Calculation for 2nd order (1-dimensional).
Double_t wert=0., rest=0.;
Int_t n1;
for (n1=0; n1<=1; n1++) {
wert += ma_1dim[n1] * psi_1dim[n1];
rest += ma_1dim[n1];
}
wertges = wert * wert - rest;
return;
}
void MWavelet::Ord2_2dim(Double_t ma_2dim[2][2], Double_t psi_2dim[3][2][2]) {
// Calculation for 2nd order (2-dimensional).
Int_t n1, n2, n;
Double_t wert, rest;
for (n=0; n<3; n++) {
wert = rest =0.;
for (n1=0; n1<=1; n1++) {
for (n2=0; n2<=1; n2++) {
wert += ma_2dim[n1][n2] * psi_2dim[n][n1][n2];
rest += ma_2dim[n1][n2];
}
}
wert = wert * wert - rest;
wertges += wert;
}
return;
}
/*
void MWavelet::Ord2_3dim(Double_t ma[2][2][2], Double_t psi_3dim[7][2][2][2]) {
// Calculation for 2nd order (3-dimensional).
Int_t n, n1, n2, n3;
Double_t rest, wert;
for (n=0; n<7; n++) {
wert = rest = 0.;
for (n1=0; n1<=1; n1++) {
for (n2=0; n2<=1; n2++) {
for (n3=0; n3<=1; n3++) {
wert += ma_3dim[n1][n2][n3] * psi_3dim[n][n1][n2][n3];
rest += ma_3dim[n1][n2][n3];
}
}
}
wert = wert*wert - rest;
wertges += wert;
}
return;
}
*/
void MWavelet::Ord3_1dim(Double_t ma_1dim[2], Double_t psi_1dim[2]) {
// Calculation for 3rd order (1-dimensional).
Double_t wert, wert1, wert2=0;
Int_t n1, k1, m1;
wert2=0.;
for (n1=0; n1<=1; n1++) {
ma_1dim[n1]--;
wert1=0.;
for (k1=0; k1<=1; k1++) {
ma_1dim[k1]--;
wert=0.;
for (m1=0; m1<=1; m1++) {
wert += ma_1dim[m1] * psi_1dim[m1];
}
ma_1dim[k1]++;
wert *= ma_1dim[k1] * psi_1dim[k1];
wert1 += wert;
}
ma_1dim[n1]++;
wert1 *= ma_1dim[n1] * psi_1dim[n1];
wert2 += wert1;
}
wertges = wert2;
return;
}
void MWavelet::Ord3_2dim(Double_t ma_2dim[2][2], Double_t psi_2dim[3][2][2]) {
// Calculation for 3rd order (2-dimensional).
Int_t n1, n2, k1, k2, m1, m2, n;
Double_t wert, wert1, wert2;
for (n=0; n<3; n++) {
wert2 = 0.;
for (n1=0; n1<=1; n1++) {
for (n2=0; n2<=1; n2++) {
ma_2dim[n1][n2]--;
wert1 = 0.;
for (k1=0; k1<=1; k1++) {
for (k2=0; k2<=1; k2++) {
ma_2dim[k1][k2]--;
wert = 0.;
for (m1=0; m1<=1; m1++) {
for (m2=0; m2<=1; m2++) {
wert += ma_2dim[m1][m2] * psi_2dim[n][m1][m2];
}
}
ma_2dim[k1][k2]++;
wert *= ma_2dim[k1][k2] * psi_2dim[n][k1][k2];
wert1 += wert;
}
}
ma_2dim[n1][n2]++;
wert1 *= ma_2dim[n1][n2] * psi_2dim[n][n1][n2];
wert2 += wert1;
}
}
wertges += wert2;
}
return;
}
/*
void MWavelet::Ord3_3dim(Double_t ma_3dim[2][2][2], Double_t psi_3dim[7][2][2][2]) {
// Calculation for 3rd order (3-dimensional).
Int_t n, n1, n2, n3, k1, k2, k3, m1, m2, m3;
Double_t wert, wert1, wert2;
wertges = 0.;
for (n=0; n<7; n++) {
wert2 = 0.;
for (n1=0; n1<=1; n1++) {
for (n2=0; n2<=1; n2++) {
for (n3=0; n3<=1; n3++) {
ma_3dim[n1][n2][n3]--;
wert1 = 0.;
for (k1=0; k1<=1; k1++) {
for (k2=0; k2<=1; k2++) {
for (k3=0; k3<=1; k3++) {
ma_3dim[k1][k2][k3]--;
wert = 0.;
for (m1=0; m1<=1; m1++) {
for (m2=0; m2<=1; m2++) {
for (m3=0; m3<=1; m3++) {
wert += ma_3dim[m1][m2][m3] * psi_3dim[n][m1][m2][m3];
}
}
}
ma_3dim[k1][k2][k3]++;
wert *= ma_3dim[m1][m2][m3] * psi_3dim[n][k1][k2][k3];
wert1 += wert;
}
}
}
ma_3dim[n1][n2][n3]++;
wert1 *= ma_3dim[n1][n2][n3] * psi_3dim[n][n1][n2][n3];
}
}
}
wertges += wert2;
}
return;
}
*/
void MWavelet::Ord4_1dim(Double_t ma_1dim[2], Double_t psi_1dim[2]) {
// Calculation for 4th order (1-dimensional).
Int_t n1, k1, l1, m1;
Double_t wert, wert1, wert2 ,wert3;
wert3 = 0.;
for (l1=0; l1<=1; l1++) {
ma_1dim[l1]--;
wert2 = 0.;
for (n1=0; n1<=1; n1++) {
ma_1dim[n1]--;
wert1=0.;
for (k1=0; k1<=1; k1++) {
ma_1dim[k1]--;
wert=0.;
for (m1=0; m1<=1; m1++) {
wert += ma_1dim[m1] * psi_1dim[m1];
}
ma_1dim[k1]++;
wert *= ma_1dim[k1] * psi_1dim[k1];
wert1 += wert;
}
ma_1dim[n1]++;
wert1 *= ma_1dim[n1] * psi_1dim[n1];
wert2 += wert1;
}
ma_1dim[l1]++;
wert2 *= ma_1dim[l1] * psi_1dim[l1];
wert3 += wert2;
}
wertges = wert3;
return;
}
void MWavelet::Ord4_2dim(Double_t ma_2dim[2][2], Double_t psi_2dim[3][2][2]) {
// Calculation for 4th order (2-dimensional).
Int_t k1, k2, l1, l2, m1, m2, n1, n2, n;
Double_t wert, wert1, wert2, wert3;
for (n=0; n<3; n++) {
wert3 = 0.;
for (l1=0; l1<=1; l1++) {
for (l2=0; l2<=1; l2++) {
ma_2dim[l1][l2]--;
wert2 = 0.;
for (n1=0; n1<=1; n1++) {
for (n2=0; n2<=1; n2++) {
ma_2dim[n1][n2]--;
wert1 = 0.;
for (k1=0; k1<=1; k1++) {
for (k2=0; k2<=1; k2++) {
ma_2dim[k1][k2]--;
wert = 0.;
for (m1=0; m1<=1; m1++) {
for (m2=0; m2<=1; m2++) {
wert += ma_2dim[m1][m2] * psi_2dim[n][m1][m2];
}
}
ma_2dim[k1][k2]++;
wert *= ma_2dim[k1][k2] * psi_2dim[n][k1][k2];
wert1 += wert;
}
}
ma_2dim[n1][n2]++;
wert1 *= ma_2dim[n1][n2] * psi_2dim[n][n1][n2];
wert2 += wert1;
}
}
ma_2dim[l1][l2]++;
wert2 *= ma_2dim[l1][l2] * psi_2dim[n][l1][l2];
wert3 += wert2;
}
}
wertges += wert3;
}
return;
}
/*
void MWavelet::Ord4_3dim(Double_t ma_3dim[2][2][2], Double_t psi_3dim[7][2][2][2]) {
// Calculation for 4th order (3-dimensional).
Int_t n, n1, n2, n3, k1, k2, k3, l1, l2, l3, m1, m2 ,m3;
Double_t wert, wert1, wert2, wert3;
wertges = 0.;
for (n=0; n<7; n++) {
wert3 = 0.;
for (l1=0; l1<=1; l1++) {
for (l2=0; l2<=1; l2++) {
for (l3=0; l3<=1; l3++) {
ma_3dim[l1][l2][l3]--;
wert2 = 0.;
for (n1=0; n1<=1; n1++) {
for (n2=0; n2<=1; n2++) {
for (n3=0; n3<=1; n3++) {
ma_3dim[n1][n2][n3]--;
wert1 = 0.;
for (k1=0; k1<=1; k1++) {
for (k2=0; k2<=1; k2++) {
for (k3=0; k3<=1; k3++) {
ma_3dim[k1][k2][k3]--;
wert = 0.;
for (m1=0; m1<=1; m1++) {
for (m2=0; m2<=1; m2++) {
for (m3=0; m3<=1; m3++) {
wert += ma_3dim[m1][m2][m3] * psi_3dim[n][m1][m2][m3];
}
}
}
ma_3dim[k1][k2][k3]++;
wert *= ma_3dim[k1][k2][k3] * psi_3dim[n][k1][k2][k3];
wert1 += wert;
}
}
}
ma_3dim[n1][n2][n3]++;
wert1 *= ma_3dim[n1][n2][n3] * psi_3dim[n][n1][n2][n3];
wert2 += wert1;
}
}
}
ma_3dim[l1][l2][l3]++;
wert2 *= ma_3dim[l1][l2][l3] * psi_3dim[n][l1][l2][l3];
wert3 += wert2;
}
}
}
wertges = wert3;
}
return;
}
*/
// Umkehrung der Transformation (nur fuer 512 bins!)
void MWavelet::Reverse(TH1 *in) {
// Calculates inverse wavelet transform. Generates out of a given
// wavelet coefficient distribution a histogram. Works only for 512 bins!
Int_t durchlauf, vorzeichenwechsler, lesebin, schreibbin, teil;
Double_t wert;
if (revert) {
delete revert;
revert = NULL;
}
if (!revert) {
revert = new TH1D("revert", "berechnete Originalfunktion", 512, 0., 1.);
}
// Fill global mean.
wert = in->GetBinContent(1);
for (schreibbin=1; schreibbin<=512; schreibbin++) {
revert->SetBinContent(schreibbin, wert);
}
// Fill scales.
for (durchlauf=1; durchlauf<=9; durchlauf++) {
teil=-1;
for (lesebin=Int_t(pow(2, durchlauf-1)+1); lesebin<=pow(2, durchlauf); lesebin++) {
for (vorzeichenwechsler=0; vorzeichenwechsler<=1; vorzeichenwechsler++) {
teil++;
for (schreibbin=1; schreibbin<=pow(2, 9-durchlauf); schreibbin++) {
wert = in->GetBinContent(lesebin);
revert->AddBinContent(Int_t(schreibbin+teil*pow(2, 9-durchlauf)), pow(-1., vorzeichenwechsler)*wert);
}
}
}
}
return;
}
void MWavelet::Draw(TH1 *input_histo, Int_t option) {
// Draw routine for the calculated wavelet moments.
if (original) delete original;
if (wav2) delete wav2;
if (wav3) delete wav3;
if (wav4) delete wav4;
if (canvas) canvas->Clear();
else {
canvas = new TCanvas("Wavelets");
canvas->SetFillColor(0);
}
original = new TPad("original", "original", 0.025, 0.525, 0.475, 0.975);
original->SetFillColor(0);
wav2 = new TPad("wav2", "wav2", 0.525, 0.525, 0.975, 0.975);
wav2->SetFillColor(0);
wav_mom_ord2->SetMarkerStyle(8);
wav_mom_ord2->SetMarkerSize(0.5);
wav3 = new TPad("wav3", "wav3", 0.025, 0.025, 0.475, 0.475);
wav3->SetFillColor(0);
wav_mom_ord3->SetMarkerStyle(8);
wav_mom_ord3->SetMarkerSize(0.5);
wav4 = new TPad("wav4", "wav4", 0.525, 0.025, 0.975, 0.475);
wav4->SetFillColor(0);
wav_mom_ord4->SetMarkerStyle(8);
wav_mom_ord4->SetMarkerSize(0.5);
original->Draw();
wav2->Draw();
wav3->Draw();
wav4->Draw();
original->cd();
input_histo->Draw();
wav2->cd();
wav_mom_ord2->Draw("p");
wav3->cd();
wav_mom_ord3->Draw("p");
wav4->cd();
wav_mom_ord4->Draw("p");
canvas->Update();
return;
}
void MWavelet::FillProf() {
// Creates profile histograms for different orders.
Int_t j;
if (!prof_ord2 /* accordingly other profiles shouldn't be there, either */) {
sprintf(histo, "%s__prof_ord2", name);
prof_ord2 = new TProfile(histo, "Wavelet Momente Profile 2. Ordnung", jmax, 0., (Double_t)jmax);
sprintf(histo, "%s__prof_ord3", name);
prof_ord3 = new TProfile(histo, "Wavelet Momente Profile 3. Ordnung", jmax, 0., (Double_t)jmax);
sprintf(histo, "%s__prof_ord4", name);
prof_ord4 = new TProfile(histo, "Wavelet Momente Profile 4. Ordnung", jmax, 0., (Double_t)jmax);
sprintf(histo, "%s__ord2", name);
ord2 = new TNtuple(histo, "ord2", "scale1:scale2:scale3:scale4:scale5:scale6:scale7:scale8:scale9");
sprintf(histo, "%s__ord3", name);
ord3 = new TNtuple(histo, "ord3", "scale1:scale2:scale3:scale4:scale5:scale6:scale7:scale8:scale9");
sprintf(histo, "%s__ord4", name);
ord4 = new TNtuple(histo, "ord4", "scale1:scale2:scale3:scale4:scale5:scale6:scale7:scale8:scale9");
}
ord2->Fill(mo[2][0], mo[2][1], mo[2][2], mo[2][3], mo[2][4], mo[2][5], mo[2][6], mo[2][7], mo[2][8]);
ord3->Fill(mo[3][0], mo[3][1], mo[3][2], mo[3][3], mo[3][4], mo[3][5], mo[3][6], mo[3][7], mo[3][8]);
ord4->Fill(mo[4][0], mo[4][1], mo[4][2], mo[4][3], mo[4][4], mo[4][5], mo[4][6], mo[4][7], mo[4][8]);
for (j = 0; j < jmax; j++) {
prof_ord2->Fill((Axis_t)j, (Axis_t)mo[2][j]);
prof_ord3->Fill((Axis_t)j, (Axis_t)mo[3][j]);
prof_ord4->Fill((Axis_t)j, (Axis_t)mo[4][j]);
}
return;
}
void MWavelet::SaveMom(TH1 *input_histo, Text_t *name, TObject *obj1, TObject *obj2, TObject *obj3, TObject *obj4, TObject *obj5) {
// Saves the input histogram, the different wavelet moments, and
// any additionally given objects (up to 5) under the given name into a ROOT-file.
TFile *file_mom = new TFile(name, "RECREATE", "Wavelet Momente");
input_histo->Write();
wav_mom_ord2->Write();
wav_mom_ord3->Write();
wav_mom_ord4->Write();
if (obj1) obj1->Write();
if (obj2) obj2->Write();
if (obj3) obj3->Write();
if (obj4) obj4->Write();
if (obj5) obj5->Write();
file_mom->Close();
return;
}
void MWavelet::SaveProf(Text_t *name, TObject *obj1, TObject *obj2, TObject *obj3, TObject *obj4, TObject *obj5) {
// Saves the different profile histograms, the ntuples of different orders, and
// any additionally given objects (up to 5) under the given name into a ROOT-file.
TFile *file_prof = new TFile(name, "RECREATE", "Wavelet Momente (Profiles, Bins)");
file_prof->cd();
prof_ord2->Write();
prof_ord3->Write();
prof_ord4->Write();
ord2->Write();
ord3->Write();
ord4->Write();
if (obj1) obj1->Write();
if (obj2) obj2->Write();
if (obj3) obj3->Write();
if (obj4) obj4->Write();
if (obj5) obj5->Write();
file_prof->Close();
return;
}
TProfile* MWavelet::GetProf_Ord2(void) {
// Returns the pointer to the profile histogram prof_ord2.
return prof_ord2;
}
TProfile* MWavelet::GetProf_Ord3(void) {
// Returns the pointer to the profile histogram prof_ord3.
return prof_ord3;
}
TProfile* MWavelet::GetProf_Ord4(void) {
// Returns the pointer to the profile histogram prof_ord4.
return prof_ord4;
}
TNtuple* MWavelet::GetOrd2(void) {
// Returns the pointer to the ntuple ord2.
return ord2;
}
TNtuple* MWavelet::GetOrd3(void) {
// Returns the pointer to the ntuple ord3.
return ord3;
}
TNtuple* MWavelet::GetOrd4(void) {
// Returns the pointer to the ntuple ord4.
return ord4;
}
ROOT page - Home page - Class index - Top of the page
This page has been automatically generated. If you have any comments or suggestions about the page layout send a mail to ROOT support, or contact the developers with any questions or problems regarding ROOT.