newly added for test purposes.
This commit is contained in:
nemu
parent
88c7978ab9
commit
8f817c6cbc
99
src/tests/t0NotEqFirstGoodData/Makefile.PAddPoissonNoise
Normal file
99
src/tests/t0NotEqFirstGoodData/Makefile.PAddPoissonNoise
Normal file
@ -0,0 +1,99 @@
|
||||
#---------------------------------------------------
|
||||
# Makefile.PAddPoissonNoise
|
||||
#
|
||||
# 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 = $(shell uname)
|
||||
|
||||
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 = -Wall -Wno-trigraphs -fPIC
|
||||
INCLUDES = -I../include
|
||||
LD = g++
|
||||
LDFLAGS = -g
|
||||
SOFLAGS = -O -shared
|
||||
endif
|
||||
|
||||
# -- Darwin
|
||||
ifeq ($(OS),DARWIN)
|
||||
CXX = g++
|
||||
CXXFLAGS = -Wall -Wno-trigraphs -fPIC
|
||||
INCLUDES = -I../include
|
||||
LD = g++
|
||||
LDFLAGS = -g
|
||||
SOFLAGS = -dynamic
|
||||
endif
|
||||
|
||||
# the output from the root-config script:
|
||||
CXXFLAGS += $(ROOTCFLAGS)
|
||||
LDFLAGS +=
|
||||
|
||||
# the ROOT libraries (G = graphic)
|
||||
LIBS = $(ROOTLIBS) -lXMLParser
|
||||
GLIBS = $(ROOTGLIBS) -lXMLParser
|
||||
|
||||
# some definitions: headers (used to generate *Dict* stuff), sources, objects,...
|
||||
OBJS =
|
||||
OBJS += PAddPoissonNoise.o PAddPoissonNoiseDict.o
|
||||
|
||||
SHLIB = libPAddPoissonNoise.so
|
||||
|
||||
# make the shared lib:
|
||||
#
|
||||
all: $(SHLIB)
|
||||
|
||||
$(SHLIB): $(OBJS)
|
||||
@echo "---> Building shared library $(SHLIB) ..."
|
||||
/bin/rm -f $(SHLIB)
|
||||
$(LD) $(OBJS) $(SOFLAGS) -o $(SHLIB) $(LIBS)
|
||||
@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) *Dict* core*
|
||||
@echo "---> removing $(OBJS)"
|
||||
|
||||
#
|
||||
$(OBJS): %.o: %.cpp
|
||||
$(CXX) $(INCLUDES) $(CXXFLAGS) -c $<
|
||||
|
||||
PAddPoissonNoiseDict.cpp: PAddPoissonNoise.h PAddPoissonNoiseLinkDef.h
|
||||
@echo "Generating dictionary $@..."
|
||||
rootcint -f $@ -c -p $^
|
||||
|
||||
install: all
|
||||
@echo "Installing shared lib: libPAddPoissonNoise.so ( you must be root ;-) )"
|
||||
ifeq ($(OS),LINUX)
|
||||
cp -pv $(SHLIB) $(ROOTSYS)/lib
|
||||
cp -pv PAddPoissonNoise.h $(ROOTSYS)/include
|
||||
endif
|
||||
173
src/tests/t0NotEqFirstGoodData/PAddPoissonNoise.cpp
Normal file
173
src/tests/t0NotEqFirstGoodData/PAddPoissonNoise.cpp
Normal file
@ -0,0 +1,173 @@
|
||||
/***************************************************************************
|
||||
|
||||
PAddPoissonNoise.cpp
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2009 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 <iostream>
|
||||
using namespace std;
|
||||
|
||||
#include <TMath.h>
|
||||
|
||||
#include "PAddPoissonNoise.h"
|
||||
|
||||
ClassImp(PAddPoissonNoise)
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Constructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p> Constructor.
|
||||
*
|
||||
* \param seed for the random number generator
|
||||
*/
|
||||
PAddPoissonNoise::PAddPoissonNoise(UInt_t seed)
|
||||
{
|
||||
fValid = true;
|
||||
|
||||
fRandom = new TRandom2(seed);
|
||||
if (fRandom == 0) {
|
||||
fValid = false;
|
||||
}
|
||||
|
||||
fSquareRoot = 0.0;
|
||||
fAlxm = 0.0;
|
||||
fG = 0.0;
|
||||
fOldMean = -1.0;
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// Destructor
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p> Destructor.
|
||||
*
|
||||
*/
|
||||
PAddPoissonNoise::~PAddPoissonNoise()
|
||||
{
|
||||
if (fRandom) {
|
||||
delete fRandom;
|
||||
fRandom = 0;
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// SetSeed (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Sets the seed of the random generator.
|
||||
*
|
||||
* \param seed for the random number generator
|
||||
*/
|
||||
void PAddPoissonNoise::SetSeed(UInt_t seed)
|
||||
{
|
||||
if (!fValid)
|
||||
return;
|
||||
|
||||
fRandom->SetSeed(seed);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// AddNoise (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Adds Poisson noise to a single input value.
|
||||
*
|
||||
* \param val input value.
|
||||
*/
|
||||
Double_t PAddPoissonNoise::AddNoise(Double_t val)
|
||||
{
|
||||
return PoiDev(val);
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// AddNoise (public)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Adds Poisson noise to a histogram.
|
||||
*
|
||||
* \param histo histogram to which Poisson noise is added.
|
||||
*/
|
||||
void PAddPoissonNoise::AddNoise(TH1F *histo)
|
||||
{
|
||||
if (histo == 0)
|
||||
return;
|
||||
|
||||
Double_t dval;
|
||||
for (UInt_t i=1; i<histo->GetEntries(); i++) {
|
||||
dval = histo->GetBinContent(i);
|
||||
histo->SetBinContent(i, PoiDev(dval));
|
||||
}
|
||||
}
|
||||
|
||||
//--------------------------------------------------------------------------
|
||||
// PoiDev (private)
|
||||
//--------------------------------------------------------------------------
|
||||
/**
|
||||
* <p>Calculates Poisson noise. Taken from Numerical Recipes: method NR::poidev, chapter 7, p. 298
|
||||
* Numerical Recipes in C++, 2nd Edition 2002.
|
||||
*
|
||||
* Returns as a floating-point number an integer value that is a random deviate drawn from a Poisson
|
||||
* distribution of mean "mean".
|
||||
*
|
||||
* \param mean value to which Poisson noise should be added.
|
||||
*/
|
||||
Double_t PAddPoissonNoise::PoiDev(const Double_t &mean)
|
||||
{
|
||||
Double_t em, t, y;
|
||||
|
||||
if (mean < 12.0) {
|
||||
if (mean != fOldMean) {
|
||||
fOldMean = mean;
|
||||
fG = TMath::Exp(-mean);
|
||||
}
|
||||
em = -1.0;
|
||||
t = 1.0;
|
||||
do {
|
||||
++em;
|
||||
t *= fRandom->Rndm();
|
||||
} while (t > fG);
|
||||
} else {
|
||||
if (mean != fOldMean) {
|
||||
fOldMean = mean;
|
||||
fSquareRoot = TMath::Sqrt(2.0*mean);
|
||||
fAlxm = TMath::Log(mean);
|
||||
fG = mean*fAlxm-TMath::LnGamma(mean+1.0);
|
||||
}
|
||||
do {
|
||||
do {
|
||||
y = TMath::Tan(TMath::Pi()*fRandom->Rndm());
|
||||
em = fSquareRoot*y+mean;
|
||||
} while (em < 0.0);
|
||||
em = TMath::Floor(em);
|
||||
t = 0.9*(1.0+y*y)*TMath::Exp(em*fAlxm-TMath::LnGamma(em+1.0)-fG);
|
||||
} while (fRandom->Rndm() > t);
|
||||
}
|
||||
|
||||
return em;
|
||||
}
|
||||
64
src/tests/t0NotEqFirstGoodData/PAddPoissonNoise.h
Normal file
64
src/tests/t0NotEqFirstGoodData/PAddPoissonNoise.h
Normal file
@ -0,0 +1,64 @@
|
||||
/***************************************************************************
|
||||
|
||||
PAddPoissonNoise.h
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2009 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. *
|
||||
***************************************************************************/
|
||||
|
||||
#ifndef _PADDPOISSONNOISE_H_
|
||||
#define _PADDPOISSONNOISE_H_
|
||||
|
||||
#include <TObject.h>
|
||||
#include <TH1F.h>
|
||||
#include <TRandom2.h>
|
||||
|
||||
class PAddPoissonNoise : public TObject
|
||||
{
|
||||
public:
|
||||
PAddPoissonNoise(UInt_t seed = 0);
|
||||
virtual ~PAddPoissonNoise();
|
||||
|
||||
virtual Bool_t IsValid() { return fValid; }
|
||||
virtual void SetSeed(UInt_t seed);
|
||||
virtual Double_t AddNoise(Double_t val);
|
||||
virtual void AddNoise(TH1F *histo);
|
||||
|
||||
private:
|
||||
Bool_t fValid;
|
||||
TRandom2 *fRandom;
|
||||
|
||||
Double_t fSquareRoot;
|
||||
Double_t fAlxm;
|
||||
Double_t fG;
|
||||
Double_t fOldMean;
|
||||
|
||||
virtual Double_t PoiDev(const Double_t &mean);
|
||||
|
||||
ClassDef(PAddPoissonNoise, 1)
|
||||
};
|
||||
|
||||
#endif // _PADDPOISSONNOISE_H_
|
||||
40
src/tests/t0NotEqFirstGoodData/PAddPoissonNoiseLinkDef.h
Normal file
40
src/tests/t0NotEqFirstGoodData/PAddPoissonNoiseLinkDef.h
Normal file
@ -0,0 +1,40 @@
|
||||
/***************************************************************************
|
||||
|
||||
PAddPoissonNoiseLinkDef.h
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2007 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. *
|
||||
***************************************************************************/
|
||||
|
||||
#ifdef __CINT__
|
||||
|
||||
#pragma link off all globals;
|
||||
#pragma link off all classes;
|
||||
#pragma link off all functions;
|
||||
|
||||
#pragma link C++ class PAddPoissonNoise+;
|
||||
|
||||
#endif
|
||||
151
src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C
Normal file
151
src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C
Normal file
@ -0,0 +1,151 @@
|
||||
/***************************************************************************
|
||||
|
||||
t0NotEqFirstGoodData.C
|
||||
|
||||
Author: Andreas Suter
|
||||
e-mail: andreas.suter@psi.ch
|
||||
|
||||
$Id$
|
||||
|
||||
***************************************************************************/
|
||||
|
||||
/***************************************************************************
|
||||
* Copyright (C) 2009 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. *
|
||||
***************************************************************************/
|
||||
|
||||
void t0NotEqFirstGoodData()
|
||||
{
|
||||
// load library
|
||||
gSystem->Load("$ROOTSYS/lib/libPAddPoissonNoise");
|
||||
|
||||
// generate data
|
||||
TFolder *histosFolder;
|
||||
TFolder *decayAnaModule;
|
||||
TFolder *runInfo;
|
||||
|
||||
histosFolder = gROOT->GetRootFolder()->AddFolder("histos", "Histograms");
|
||||
gROOT->GetListOfBrowsables()->Add(histosFolder, "histos");
|
||||
decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms");
|
||||
|
||||
// feed run info header
|
||||
runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo");
|
||||
gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo");
|
||||
header = new TLemRunHeader();
|
||||
header->SetRunTitle("010101 - test");
|
||||
header->SetLemSetup("trivial");
|
||||
header->SetRunNumber(10101);
|
||||
header->SetStartTime(0);
|
||||
header->SetStopTime(1);
|
||||
header->SetModeratorHV(32.0, 0.01);
|
||||
header->SetSampleHV(0.0, 0.01);
|
||||
header->SetImpEnergy(31.8);
|
||||
header->SetSampleTemperature(0.2, 0.001);
|
||||
header->SetSampleBField(30000, 0.1);
|
||||
header->SetTimeResolution(0.1953125);
|
||||
header->SetNChannels(66601);
|
||||
header->SetNHist(4);
|
||||
header->SetCuts("none");
|
||||
header->SetModerator("none");
|
||||
Double_t tt0[4] = {3419.0, 3419.0, 3419.0, 3419.0};
|
||||
header->SetTimeZero(tt0);
|
||||
runInfo->Add(header); //add header to RunInfo folder
|
||||
|
||||
// create histos
|
||||
UInt_t t0[4] = {3419, 3419, 3419, 3419};
|
||||
UInt_t n0[4] = {200.0, 205.0, 203.0, 198.0};
|
||||
UInt_t bkg[4] = {11.0, 11.0, 5.0, 8.0};
|
||||
const Double_t tau = 2197.019; // ns
|
||||
// asymmetry related stuff
|
||||
const Double_t timeResolution = 0.1953125; // ns
|
||||
Double_t a[4] = {0.16, 0.161, 0.16, 0.159};
|
||||
Double_t phase[4] = {13.0*TMath::Pi()/180.0, 103.0*TMath::Pi()/180.0, 193.0*TMath::Pi()/180.0, 283.0*TMath::Pi()/180.0};
|
||||
const Double_t gamma = 0.00001355; // gamma/(2pi)
|
||||
Double_t bb = 15000.0; // field in Gauss
|
||||
// Double_t bb = 100.0; // field in Gauss
|
||||
Double_t sigma = 1.05/1000.0; // Gaussian sigma in 1/ns
|
||||
|
||||
TH1F *histo[8];
|
||||
char str[128];
|
||||
for (UInt_t i=0; i<4; i++) {
|
||||
sprintf(str, "hDecay0%d", (Int_t)i);
|
||||
histo[i] = new TH1F(str, str, 66601, -0.5, 66600.5);
|
||||
sprintf(str, "hDecay2%d", (Int_t)i);
|
||||
histo[i+4] = new TH1F(str, str, 66601, -0.5, 66600.5);
|
||||
}
|
||||
Double_t time;
|
||||
Double_t dval;
|
||||
for (UInt_t i=0; i<4; i++) {
|
||||
for (UInt_t j=1; j<66601; j++) {
|
||||
if (j<t0[i]) {
|
||||
histo[i]->SetBinContent(j, bkg[i]);
|
||||
} else {
|
||||
time = (Double_t)(j-t0[i])*timeResolution;
|
||||
dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a[i]*TMath::Exp(-0.5*TMath::Power(time*sigma,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb*time+phase[i]))+(Double_t)bkg[i];
|
||||
histo[i]->SetBinContent(j, (UInt_t)dval);
|
||||
}
|
||||
}
|
||||
}
|
||||
|
||||
// add a promp peak
|
||||
Double_t ampl = 500.0;
|
||||
Double_t promptLambda = 500.0/1000.0; // Lorentzain in 1/ns
|
||||
for (UInt_t i=0; i<4; i++) {
|
||||
for (UInt_t j=1; j<66601; j++) {
|
||||
dval = histo[i]->GetBinContent(j);
|
||||
time = (Double_t)(j-t0[i])*timeResolution;
|
||||
dval += ampl*TMath::Exp(-fabs(time)*promptLambda);
|
||||
histo[i]->SetBinContent(j, (UInt_t)dval);
|
||||
histo[i+4]->SetBinContent(j, (UInt_t)dval);
|
||||
}
|
||||
}
|
||||
|
||||
// add Poisson noise
|
||||
PAddPoissonNoise *addNoise = new PAddPoissonNoise();
|
||||
if (!addNoise->IsValid()) {
|
||||
cerr << endl << "**ERROR** while invoking PAddPoissonNoise" << endl;
|
||||
return;
|
||||
}
|
||||
for (UInt_t i=0; i<4; i++) {
|
||||
addNoise->AddNoise(histo[i]);
|
||||
}
|
||||
for (UInt_t i=0; i<4; i++) {
|
||||
for (UInt_t j=1; j<histo[i]->GetEntries(); j++) {
|
||||
histo[i+4]->SetBinContent(j, histo[i]->GetBinContent(j));
|
||||
}
|
||||
}
|
||||
|
||||
for (UInt_t i=0; i<8; i++)
|
||||
decayAnaModule->Add(histo[i]);
|
||||
|
||||
// write file
|
||||
TFile *fout = new TFile("010105.root", "RECREATE", "Midas Fake Histograms");
|
||||
if (fout == 0) {
|
||||
cout << endl << "**ERROR** Couldn't create ROOT file";
|
||||
cout << endl << endl;
|
||||
exit(0);
|
||||
}
|
||||
|
||||
fout->cd();
|
||||
runInfo->Write();
|
||||
histosFolder->Write();
|
||||
fout->Close();
|
||||
delete fout;
|
||||
|
||||
delete [] histo;
|
||||
}
|
||||
Loading…
x
Reference in New Issue
Block a user