first fourier playground
This commit is contained in:
nemu
parent
b96ac14a18
commit
13c48c72f0
104
src/tests/fourier/Makefile.fourier
Normal file
104
src/tests/fourier/Makefile.fourier
Normal file
@ -0,0 +1,104 @@
|
||||
#---------------------------------------------------
|
||||
# Makefile.fourier
|
||||
#
|
||||
# Author: Andreas Suter
|
||||
# e-mail: andreas.suter@psi.ch
|
||||
#
|
||||
# $Id$
|
||||
#---------------------------------------------------
|
||||
|
||||
#---------------------------------------------------
|
||||
# get compilation and library flags from root-config
|
||||
|
||||
ROOTCFLAGS = $(shell $(ROOTSYS)/bin/root-config --cflags)
|
||||
ROOTLIBS = $(shell $(ROOTSYS)/bin/root-config --libs)
|
||||
ROOTGLIBS = $(shell $(ROOTSYS)/bin/root-config --glibs)
|
||||
|
||||
#---------------------------------------------------
|
||||
# depending on the architecture, choose the compiler,
|
||||
# linker, and the flags to use
|
||||
#
|
||||
|
||||
OSTYPE = linux
|
||||
|
||||
ifeq ($(OSTYPE),linux)
|
||||
OS = LINUX
|
||||
endif
|
||||
ifeq ($(OSTYPE),linux-gnu)
|
||||
OS = LINUX
|
||||
endif
|
||||
ifeq ($(OSTYPE),darwin)
|
||||
OS = DARWIN
|
||||
endif
|
||||
|
||||
# -- Linux
|
||||
ifeq ($(OS),LINUX)
|
||||
CXX = g++
|
||||
CXXFLAGS = -g -Wall -fPIC
|
||||
PMUSRPATH = ./include
|
||||
MNPATH = $(ROOTSYS)/include
|
||||
INCLUDES = -I $(PMUSRPATH) -I $(MNPATH)
|
||||
LD = g++
|
||||
LDFLAGS = -g
|
||||
endif
|
||||
|
||||
# -- Darwin
|
||||
ifeq ($(OS),DARWIN)
|
||||
CXX = g++
|
||||
CXXFLAGS = -g -Wall -fPIC
|
||||
INCLUDES = -I./
|
||||
LD = g++
|
||||
LDFLAGS = -g
|
||||
endif
|
||||
|
||||
# the output from the root-config script:
|
||||
CXXFLAGS += $(ROOTCFLAGS)
|
||||
LDFLAGS +=
|
||||
|
||||
# the ROOT libraries (G = graphic)
|
||||
LIBS = $(ROOTLIBS) -lXMLParser
|
||||
GLIBS = $(ROOTGLIBS) -lXMLParser
|
||||
|
||||
# PSI libs
|
||||
PSILIBS = -lTLemRunHeader -lPMusr
|
||||
# Minuit2 lib
|
||||
MNLIB = -L$(ROOTSYS)/lib -lMinuit2
|
||||
# MathMore lib
|
||||
MMLIB = -L$(ROOTSYS)/lib -lMathMore
|
||||
|
||||
INSTALLPATH = $(HOME)/analysis/bin
|
||||
|
||||
EXEC = fourier
|
||||
|
||||
FOURIER = PMusrFourier
|
||||
|
||||
# some definitions: headers, sources, objects,...
|
||||
OBJS =
|
||||
OBJS += $(EXEC).o
|
||||
OBJS += $(FOURIER).o
|
||||
|
||||
# make the executable:
|
||||
#
|
||||
all: $(EXEC)
|
||||
|
||||
$(EXEC): $(OBJS)
|
||||
@echo "---> Building $(EXEC) ..."
|
||||
/bin/rm -f $(SHLIB)
|
||||
$(LD) $(OBJS) -o $(EXEC) $(GLIBS) $(PSILIBS) $(MNLIB) $(MMLIB)
|
||||
@echo "done"
|
||||
|
||||
# clean up: remove all object file (and core files)
|
||||
# semicolon needed to tell make there is no source
|
||||
# for this target!
|
||||
#
|
||||
clean:; @rm -f $(OBJS)
|
||||
@echo "---> removing $(OBJS)"
|
||||
|
||||
#
|
||||
$(OBJS): %.o: %.cpp
|
||||
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
|
||||
|
||||
install: all
|
||||
@echo "copy $(EXEC) to $(INSTALLPATH)"
|
||||
cp -p $(EXEC) $(INSTALLPATH)
|
||||
|
||||
427
src/tests/fourier/PMusrFourier.cpp
Normal file
427
src/tests/fourier/PMusrFourier.cpp
Normal file
@ -0,0 +1,427 @@
|
||||
/***************************************************************************
|
||||
|
||||
PMusrFourier.cpp
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Suter *
|
||||
* andreas.suter@psi.c *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License for more details. *
|
||||
* *
|
||||
* You should have received a copy of the GNU General Public License *
|
||||
* along with this program; if not, write to the *
|
||||
* Free Software Foundation, Inc., *
|
||||
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||
***************************************************************************/
|
||||
|
||||
#include <cmath>
|
||||
|
||||
#include "TH1F.h"
|
||||
#include "TF1.h"
|
||||
|
||||
#include "TFile.h"
|
||||
|
||||
#include "PMusrFourier.h"
|
||||
|
||||
#define PI 3.14159265358979312
|
||||
#define PI_HALF 1.57079632679489656
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Constructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \dataType tag indicating if data is histogram, asymmetry, ...
|
||||
* \data vector with the real data
|
||||
*/
|
||||
PMusrFourier::PMusrFourier(int dataType, PDoubleVector &data, double timeResolution,
|
||||
double startTime, double endTime, unsigned int rebin,
|
||||
bool estimateN0AndBkg) :
|
||||
fDataType(dataType), fTimeResolution(timeResolution),
|
||||
fStartTime(startTime), fEndTime(endTime), fRebin(rebin)
|
||||
{
|
||||
// init stuff
|
||||
fData = data;
|
||||
if (fData.empty()) {
|
||||
cout << endl << "**ERROR** PMusrFourier::PMusrFourier: no valid data" << endl << endl;
|
||||
fValid = false;
|
||||
return;
|
||||
}
|
||||
|
||||
if ((fStartTime < 0.0) || (fEndTime < 0.0)) {
|
||||
cout << endl << "**ERROR** PMusrFourier::PMusrFourier: no valid start or end time." << endl << endl;
|
||||
fValid = false;
|
||||
return;
|
||||
}
|
||||
|
||||
fValid = true;
|
||||
fIn = 0;
|
||||
fOut = 0;
|
||||
|
||||
// if endTime == 0 set is to the last time slot
|
||||
if (fEndTime == 0.0)
|
||||
fEndTime = (fData.size()-1)*fTimeResolution;
|
||||
else
|
||||
fEndTime *= 1000.0; // us -> ns
|
||||
|
||||
fStartTime *= 1000.0; // us -> ns
|
||||
|
||||
// swap start and end time if necessary
|
||||
if (fStartTime > fEndTime) {
|
||||
double keep = fStartTime;
|
||||
fStartTime = fEndTime;
|
||||
fEndTime = keep;
|
||||
}
|
||||
|
||||
// rebin data if necessary
|
||||
if (fRebin > 1) {
|
||||
fTimeResolution *= fRebin;
|
||||
double dval = 0.0;
|
||||
unsigned int rebinCounter = 0;
|
||||
for (unsigned int i=0; i<fData.size(); i++) {
|
||||
dval += fData[i];
|
||||
rebinCounter++;
|
||||
if (rebinCounter == fRebin) {
|
||||
fDataRebinned.push_back(dval);
|
||||
dval = 0.0;
|
||||
rebinCounter = 0;
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
cout << endl << "dB = " << 1.0/(F_GAMMA_BAR_MUON * (fEndTime-fStartTime)) << " (G), Bmax = " << 1.0/(F_GAMMA_BAR_MUON*fTimeResolution) << " (G)" << endl;
|
||||
|
||||
|
||||
// try to estimate N0 and Bkg just out of the raw data
|
||||
if (estimateN0AndBkg) {
|
||||
EstimateN0AndBkg();
|
||||
}
|
||||
|
||||
// calculate start and end bin
|
||||
unsigned int start = (unsigned int)(fStartTime/fTimeResolution);
|
||||
unsigned int end = (unsigned int)(fEndTime/fTimeResolution);
|
||||
|
||||
// allocate necessary memory
|
||||
fNoOfData = end-start;
|
||||
fIn = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfData);
|
||||
fOut = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfData);
|
||||
|
||||
// check if memory allocation has been successful
|
||||
if ((fIn == 0) || (fOut == 0)) {
|
||||
fValid = false;
|
||||
return;
|
||||
}
|
||||
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
fIn[i][0] = fDataRebinned[i+start];
|
||||
fIn[i][1] = 0.0;
|
||||
}
|
||||
|
||||
// has to be removed after testing
|
||||
TH1F test_raw_in("test_raw_in", "test_raw_in", fNoOfData+1,
|
||||
fStartTime - fTimeResolution/2.0, fEndTime + fTimeResolution/2.0);
|
||||
for (unsigned int i=0; i<fNoOfData; i++)
|
||||
test_raw_in.SetBinContent(i+1, fIn[i][0]);
|
||||
TFile f("test_raw_in.root", "RECREATE");
|
||||
test_raw_in.Write();
|
||||
f.Close();
|
||||
|
||||
fFFTwPlan = fftw_plan_dft_1d(fNoOfData, fIn, fOut, FFTW_FORWARD, FFTW_ESTIMATE);
|
||||
|
||||
if (!fFFTwPlan) {
|
||||
fValid = false;
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Destructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
*/
|
||||
PMusrFourier::~PMusrFourier()
|
||||
{
|
||||
cout << endl << "in ~PMusrFourier() ..." << endl;
|
||||
if (fValid) {
|
||||
fftw_destroy_plan(fFFTwPlan);
|
||||
fftw_free(fIn);
|
||||
fftw_free(fOut);
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Transform
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param apodization
|
||||
* \param filter
|
||||
* \param estimateN0AndBkg
|
||||
*/
|
||||
void PMusrFourier::Transform(int apodization, int filter)
|
||||
{
|
||||
if (!fValid)
|
||||
return;
|
||||
|
||||
if (fDataType == F_SINGLE_HISTO) {
|
||||
PrepareFFTwInputData();
|
||||
}
|
||||
|
||||
// for test only
|
||||
// keep data
|
||||
fftw_complex *data;
|
||||
data = (fftw_complex *)fftw_malloc(sizeof(fftw_complex)*fNoOfData);
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
data[i][0] = fIn[i][0];
|
||||
data[i][1] = 0.0;
|
||||
}
|
||||
|
||||
// loop over the phase
|
||||
double sum;
|
||||
TH1F sumHist("sumHist", "sumHist", 361, -180.5, 180.5);
|
||||
TH1F re("re", "re", fNoOfData+1, -0.5, (double)fNoOfData+0.5);
|
||||
TH1F im("im", "im", fNoOfData+1, -0.5, (double)fNoOfData+0.5);
|
||||
for (int p=-180; p<180; p++) {
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
// recalculate fIn including the phase
|
||||
fIn[i][0] = data[i][0]*cos(p/180.0*PI);
|
||||
fIn[i][1] = data[i][0]*sin(p/180.0*PI);
|
||||
}
|
||||
fftw_execute(fFFTwPlan);
|
||||
|
||||
if (p==51) {
|
||||
for (unsigned int j=0; j<fNoOfData; j++) {
|
||||
re.SetBinContent(j+1, fOut[j][0]);
|
||||
im.SetBinContent(j+1, fOut[j][1]);
|
||||
}
|
||||
}
|
||||
|
||||
// calculate sum of the imaginary part of fOut
|
||||
sum = 0.0;
|
||||
for (unsigned int i=0; i<fNoOfData/2; i++) {
|
||||
sum += fOut[i][1];
|
||||
}
|
||||
sumHist.SetBinContent(p+181, fabs(sum));
|
||||
}
|
||||
TFile f("test_out.root", "RECREATE");
|
||||
re.Write();
|
||||
im.Write();
|
||||
sumHist.Write();
|
||||
f.Close();
|
||||
|
||||
fftw_free(data);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// GetRealFourier
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param realFourier
|
||||
*/
|
||||
void PMusrFourier::GetRealFourier(PDoubleVector &realFourier)
|
||||
{
|
||||
// check if valid flag is set
|
||||
if (!fValid)
|
||||
return;
|
||||
|
||||
// clear realFourier vector first
|
||||
realFourier.clear();
|
||||
|
||||
// fill realFourier vector
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
realFourier.push_back(fOut[i][0]);
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// GetImaginaryFourier
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param imaginaryFourier
|
||||
*/
|
||||
void PMusrFourier::GetImaginaryFourier(PDoubleVector &imaginaryFourier)
|
||||
{
|
||||
// check if valid flag is set
|
||||
if (!fValid)
|
||||
return;
|
||||
|
||||
// clear imaginaryFourier vector first
|
||||
imaginaryFourier.clear();
|
||||
|
||||
// fill imaginaryFourier vector
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
imaginaryFourier.push_back(fOut[i][1]);
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// GetPowerFourier
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param powerFourier
|
||||
*/
|
||||
void PMusrFourier::GetPowerFourier(PDoubleVector &powerFourier)
|
||||
{
|
||||
// check if valid flag is set
|
||||
if (!fValid)
|
||||
return;
|
||||
|
||||
// clear powerFourier vector first
|
||||
powerFourier.clear();
|
||||
|
||||
// fill powerFourier vector
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
powerFourier.push_back(sqrt(fOut[i][0]*fOut[i][0]+fOut[i][0]*fOut[i][0]));
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// GetPhaseFourier
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
* \param powerFourier
|
||||
*/
|
||||
void PMusrFourier::GetPhaseFourier(PDoubleVector &phaseFourier)
|
||||
{
|
||||
// check if valid flag is set
|
||||
if (!fValid)
|
||||
return;
|
||||
|
||||
// clear phaseFourier vector first
|
||||
phaseFourier.clear();
|
||||
|
||||
// fill phaseFourier vector
|
||||
double value = 0.0;
|
||||
for (unsigned int i=0; i<fNoOfData; i++) {
|
||||
// calculate the phase
|
||||
if (fOut[i][0] == 0) {
|
||||
if (fOut[i][1] >= 0.0)
|
||||
value = PI_HALF;
|
||||
else
|
||||
value = -PI_HALF;
|
||||
} else {
|
||||
value = atan(fOut[i][1]/fOut[i][0]);
|
||||
// check sector
|
||||
if (fOut[i][0] < 0.0) {
|
||||
if (fOut[i][1] > 0.0)
|
||||
value += PI_HALF;
|
||||
else
|
||||
value -= PI_HALF;
|
||||
}
|
||||
}
|
||||
|
||||
phaseFourier.push_back(value);
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// EstimateN0AndBkg
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
*/
|
||||
void PMusrFourier::EstimateN0AndBkg()
|
||||
{
|
||||
TH1F summHisto("summHisto", "summHisto", fDataRebinned.size()+1,
|
||||
-fTimeResolution/2.0, (fDataRebinned.size()-1)*fTimeResolution +
|
||||
fTimeResolution/2.0);
|
||||
|
||||
// fill summHisto
|
||||
double value = 0.0;
|
||||
for (unsigned int i=1; i<fDataRebinned.size(); i++) {
|
||||
value += fDataRebinned[i-1];
|
||||
summHisto.SetBinContent(i, value);
|
||||
summHisto.SetBinError(i, sqrt(value));
|
||||
}
|
||||
|
||||
cout << endl << ">> max.summHisto=" << summHisto.GetBinContent(fDataRebinned.size()-1) << endl << endl;
|
||||
|
||||
// define fit function
|
||||
TF1 *func = new TF1("func", "[0]*(1-TMath::Exp(-x/[1]))+[2]*x",
|
||||
-fTimeResolution/2.0, (fDataRebinned.size()-1)*fTimeResolution +
|
||||
fTimeResolution/2.0);
|
||||
// parameter 0 ~ N0 tau
|
||||
func->SetParameter(0, summHisto.GetBinContent(fDataRebinned.size()-1));
|
||||
// parameter 1 == tau
|
||||
func->FixParameter(1, PMUON_LIFETIME*1000.0);
|
||||
// parameter 2 ~ <Bkg>
|
||||
func->SetParameter(2, summHisto.GetBinContent(fDataRebinned.size()-1)/(PMUON_LIFETIME*1000.0)*0.05);
|
||||
|
||||
// do the fit
|
||||
summHisto.Fit(func, "0QR"); // 0->no Graph, Q->quite, R->time range from function
|
||||
|
||||
// get out the parameters
|
||||
double A = func->GetParameter(0);
|
||||
double B = func->GetParameter(2);
|
||||
|
||||
cout << endl << ">> A=" << A << ", B=" << B;
|
||||
cout << endl << ">> N0/per bin=" << A/(PMUON_LIFETIME*1000.0)*fTimeResolution << ", <Bkg> per bin=" << B*fTimeResolution << endl << endl;
|
||||
|
||||
fN0 = A/(PMUON_LIFETIME*1000.0)*fTimeResolution;
|
||||
fBkg = B*fTimeResolution;
|
||||
|
||||
// clean up
|
||||
if (func) {
|
||||
delete func;
|
||||
func = 0;
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// PrepareFFTwInputData
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>
|
||||
*
|
||||
*/
|
||||
void PMusrFourier::PrepareFFTwInputData()
|
||||
{
|
||||
// 1st subtract the Bkg from the data
|
||||
for (unsigned int i=0; i<fNoOfData; i++)
|
||||
fIn[i][0] -= fBkg;
|
||||
|
||||
// 2nd remove the lifetime term
|
||||
unsigned int start = (unsigned int)(fStartTime/fTimeResolution);
|
||||
for (unsigned int i=0; i<fNoOfData; i++)
|
||||
fIn[i][0] *= exp((start+i)*fTimeResolution/(PMUON_LIFETIME*1000.0));
|
||||
|
||||
// 3rd remove the constant N0 term
|
||||
for (unsigned int i=0; i<fNoOfData; i++)
|
||||
fIn[i][0] -= fN0;
|
||||
|
||||
// has to be removed after testing
|
||||
TH1F test_in("test_in", "test_in", fNoOfData,
|
||||
fStartTime - fTimeResolution/2.0, fEndTime + fTimeResolution/2.0);
|
||||
for (unsigned int i=0; i<fNoOfData; i++)
|
||||
test_in.SetBinContent(i, fIn[i][0]);
|
||||
TFile f("test_in.root", "RECREATE");
|
||||
test_in.Write();
|
||||
f.Close();
|
||||
}
|
||||
131
src/tests/fourier/PMusrFourier.h
Normal file
131
src/tests/fourier/PMusrFourier.h
Normal file
@ -0,0 +1,131 @@
|
||||
/***************************************************************************
|
||||
|
||||
PMusrFourier.h
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2008 by Andreas Suter *
|
||||
* andreas.suter@psi.c *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License for more details. *
|
||||
* *
|
||||
* You should have received a copy of the GNU General Public License *
|
||||
* along with this program; if not, write to the *
|
||||
* Free Software Foundation, Inc., *
|
||||
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef _PMUSRFOURIER_H_
|
||||
#define _PMUSRFOURIER_H_
|
||||
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include "fftw3.h"
|
||||
|
||||
#define F_ESTIMATE_N0_AND_BKG true
|
||||
|
||||
#define F_SINGLE_HISTO 0
|
||||
#define F_ASYMMETRY 1
|
||||
|
||||
#define F_APODIZATION_NONE 0
|
||||
#define F_APODIZATION_WEAK 1
|
||||
#define F_APODIZATION_MEDIUM 2
|
||||
#define F_APODIZATION_STRONG 3
|
||||
#define F_APODIZATION_USER 4
|
||||
|
||||
#define F_FILTER_NONE 0
|
||||
#define F_FILTER_LOW_PASS 1
|
||||
#define F_FILTER_LOW_PASS_BESSEL 2
|
||||
#define F_FILTER_LOW_PASS_BUTTERWORTH 3
|
||||
#define F_FILTER_LOW_PASS_CHEBYCHEV 4
|
||||
#define F_FILTER_BAND_PASS 5
|
||||
#define F_FILTER_BAND_PASS_BESSEL 6
|
||||
#define F_FILTER_BAND_PASS_BUTTERWORTH 7
|
||||
#define F_FILTER_BAND_PASS_CHEBYCHEV 8
|
||||
#define F_FILTER_HIGH_PASS 9
|
||||
#define F_FILTER_HIGH_PASS_BESSEL 10
|
||||
#define F_FILTER_HIGH_PASS_BUTTERWORTH 11
|
||||
#define F_FILTER_HIGH_PASS_CHEBYCHEV 12
|
||||
#define F_FILTER_USER 13
|
||||
|
||||
// gamma_muon / (2 pi) = 1.355342e-5 (GHz/G)
|
||||
#define F_GAMMA_BAR_MUON 1.355342e-5
|
||||
|
||||
/**************************************************************
|
||||
!!! START: NEEDS TO BE REMOVED WHEN MERGING WITH MUSRVIEW !!!
|
||||
***************************************************************/
|
||||
// muon life time in (us), see PRL99, 032001 (2007)
|
||||
#define PMUON_LIFETIME 2.197019
|
||||
|
||||
typedef vector<double> PDoubleVector;
|
||||
/**************************************************************
|
||||
!!! END: NEEDS TO BE REMOVED WHEN MERGING WITH MUSRVIEW !!!
|
||||
***************************************************************/
|
||||
|
||||
class PMusrFourier
|
||||
{
|
||||
public:
|
||||
PMusrFourier(int dataType, PDoubleVector &data, double timeResolution,
|
||||
double startTime = 0.0, double endTime = 0.0, unsigned int binning = 1,
|
||||
bool estimateN0AndBkg = false);
|
||||
virtual ~PMusrFourier();
|
||||
|
||||
virtual void SetN0(double n0) { fN0 = n0; }
|
||||
virtual void SetBkg(double bkg) { fBkg = bkg; }
|
||||
|
||||
virtual void Transform(int apodization = -1, int filter = -1);
|
||||
|
||||
virtual double GetFieldResolution() { return fFieldResolution; }
|
||||
virtual double GetPhaseCorrection() { return fPhaseCorrection; }
|
||||
virtual void GetRealFourier(PDoubleVector &realFourier);
|
||||
virtual void GetImaginaryFourier(PDoubleVector &imaginaryFourier);
|
||||
virtual void GetPowerFourier(PDoubleVector &powerFourier);
|
||||
virtual void GetPhaseFourier(PDoubleVector &phaseFourier);
|
||||
|
||||
virtual bool IsValid() { return fValid; }
|
||||
|
||||
private:
|
||||
bool fValid;
|
||||
|
||||
int fDataType; ///< 0=Single Histo, 1=Asymmetry
|
||||
int fUseApodization; ///< 0=no apod., 1=weak apod., 2=medium apod., 3=strong apod., 4=user apod.
|
||||
int fUseFilter; ///< 0=no filter, 1=low pass, 2=band pass, 3=high pass, 4=..., n=user filter
|
||||
|
||||
double fN0;
|
||||
double fBkg;
|
||||
|
||||
double fTimeResolution;
|
||||
double fStartTime;
|
||||
double fEndTime;
|
||||
unsigned int fRebin;
|
||||
double fFieldResolution;
|
||||
double fPhaseCorrection;
|
||||
|
||||
PDoubleVector fData;
|
||||
PDoubleVector fDataRebinned;
|
||||
|
||||
unsigned int fNoOfData;
|
||||
fftw_plan fFFTwPlan;
|
||||
fftw_complex *fIn;
|
||||
fftw_complex *fOut;
|
||||
|
||||
virtual void PrepareFFTwInputData();
|
||||
virtual void EstimateN0AndBkg();
|
||||
};
|
||||
|
||||
#endif // _PMUSRFOURIER_H_
|
||||
165
src/tests/fourier/fourier.cpp
Normal file
165
src/tests/fourier/fourier.cpp
Normal file
@ -0,0 +1,165 @@
|
||||
/***************************************************************************
|
||||
|
||||
fourier.cpp
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2007 by Andreas Suter *
|
||||
* andreas.suter@psi.ch *
|
||||
* *
|
||||
* This program is free software; you can redistribute it and/or modify *
|
||||
* it under the terms of the GNU General Public License as published by *
|
||||
* the Free Software Foundation; either version 2 of the License, or *
|
||||
* (at your option) any later version. *
|
||||
* *
|
||||
* This program is distributed in the hope that it will be useful, *
|
||||
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
|
||||
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
|
||||
* GNU General Public License for more details. *
|
||||
* *
|
||||
* You should have received a copy of the GNU General Public License *
|
||||
* along with this program; if not, write to the *
|
||||
* Free Software Foundation, Inc., *
|
||||
* 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. *
|
||||
***************************************************************************/
|
||||
|
||||
#include <cstdio>
|
||||
|
||||
#include <iostream>
|
||||
#include <vector>
|
||||
using namespace std;
|
||||
|
||||
#include "TFile.h"
|
||||
#include "TFolder.h"
|
||||
#include "TH1F.h"
|
||||
|
||||
#include "TLemRunHeader.h"
|
||||
#include "PMusrFourier.h"
|
||||
|
||||
void fourier_syntax()
|
||||
{
|
||||
cout << endl << "usage: fourier <root-fln>";
|
||||
cout << endl << " <root-fln> root input file";
|
||||
cout << endl << endl;
|
||||
}
|
||||
|
||||
int main(int argc, char *argv[])
|
||||
{
|
||||
PDoubleVector data;
|
||||
double timeResolution = -1.0;
|
||||
double startTime = 0.0;
|
||||
double endTime = 0.0;
|
||||
|
||||
if (argc != 2) {
|
||||
fourier_syntax();
|
||||
return -1;
|
||||
}
|
||||
|
||||
int histoNo = -1;
|
||||
cout << endl << ">> Which histo you would like to transform (0=hDecay00, ...): ";
|
||||
cin >> histoNo;
|
||||
|
||||
// read file
|
||||
TFile f(argv[1], "READ");
|
||||
if (f.IsZombie()) {
|
||||
cout << endl << "**ERROR** Couldn't open file " << argv[1] << endl << endl;
|
||||
return -1;
|
||||
}
|
||||
|
||||
// get header info
|
||||
TFolder *folder;
|
||||
f.GetObject("RunInfo", folder);
|
||||
if (!folder) {
|
||||
cout << endl << "Couldn't obtain RunInfo from " << argv[1] << endl;
|
||||
f.Close();
|
||||
return -1;
|
||||
}
|
||||
|
||||
// read header and check if some missing run info need to be fed
|
||||
TLemRunHeader *runHeader = dynamic_cast<TLemRunHeader*>(folder->FindObjectAny("TLemRunHeader"));
|
||||
|
||||
// check if run header is valid
|
||||
if (!runHeader) {
|
||||
cout << endl << "Couldn't obtain run header info from ROOT file " << argv[1] << endl;
|
||||
f.Close();
|
||||
return -1;
|
||||
}
|
||||
|
||||
// get time resolution
|
||||
timeResolution = runHeader->GetTimeResolution();
|
||||
|
||||
// get number of histogramms
|
||||
int noOfHistos = runHeader->GetNHist();
|
||||
if (noOfHistos < histoNo+1) {
|
||||
cout << endl << "**ERROR** Only " << noOfHistos << " No of histos present (you tried to get histo No " << histoNo << ")" << endl << endl;
|
||||
f.Close();
|
||||
return -1;
|
||||
}
|
||||
|
||||
// get t0's
|
||||
Double_t *t0 = runHeader->GetTimeZero();
|
||||
char answer[128];
|
||||
cout << endl << "t0 of histo No " << histoNo << " = " << t0[histoNo] << ", is this ok (y/n)? ";
|
||||
cin >> answer;
|
||||
if (strstr(answer, "n")) {
|
||||
cout << endl << "please enter the new t0: ";
|
||||
cin >> t0[histoNo];
|
||||
}
|
||||
|
||||
|
||||
// read data ---------------------------------------------------------
|
||||
// check if histos folder is found
|
||||
f.GetObject("histos", folder);
|
||||
if (!folder) {
|
||||
cout << endl << "Couldn't obtain histos from " << argv[1] << endl;
|
||||
f.Close();
|
||||
return -1;
|
||||
}
|
||||
// get data histo
|
||||
char histoName[32];
|
||||
// read first the data which are NOT post pileup corrected
|
||||
sprintf(histoName, "hDecay%02d", histoNo);
|
||||
TH1F *histo = dynamic_cast<TH1F*>(folder->FindObjectAny(histoName));
|
||||
if (!histo) {
|
||||
cout << endl << "PRunDataHandler::ReadRootFile: Couldn't get histo " << histoName;
|
||||
f.Close();
|
||||
return -1;
|
||||
}
|
||||
// rebinning?
|
||||
unsigned int rebin;
|
||||
cout << endl << ">> Rebinning: ";
|
||||
cin >> rebin;
|
||||
// fill data
|
||||
int start = (int)t0[histoNo];
|
||||
|
||||
cout << endl << "t0=" << start;
|
||||
cout << endl << "#bins=" << histo->GetNbinsX();
|
||||
for (int i=start; i<histo->GetNbinsX(); i++) {
|
||||
data.push_back(histo->GetBinContent(i));
|
||||
}
|
||||
|
||||
f.Close();
|
||||
|
||||
cout << endl << ">> Start Time in (us): ";
|
||||
cin >> startTime;
|
||||
cout << endl << ">> End Time in (us): ";
|
||||
cin >> endTime;
|
||||
|
||||
PMusrFourier fourier(F_SINGLE_HISTO, data, timeResolution, startTime, endTime,
|
||||
rebin, F_ESTIMATE_N0_AND_BKG);
|
||||
|
||||
if (fourier.IsValid()) {
|
||||
fourier.Transform(F_APODIZATION_NONE, F_FILTER_NONE);
|
||||
} else {
|
||||
cout << endl << "**ERROR** something went wrong when invoking the PMusrFourier object ...";
|
||||
cout << endl << endl;
|
||||
}
|
||||
|
||||
return 0;
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user