From 7d384793af54b3445146b4dfecc18ae8425d5e71 Mon Sep 17 00:00:00 2001 From: nemu Date: Wed, 13 Feb 2008 14:18:46 +0000 Subject: [PATCH] newly added for testing --- src/tests/skewedGaussianTest/README | 17 + src/tests/skewedGaussianTest/fakeData.cpp | 481 ++++++++++++++++++ src/tests/skewedGaussianTest/fakeData.pro | 23 + src/tests/skewedGaussianTest/paramInput.dat | 25 + src/tests/skewedGaussianTest/skewedGaussian.C | 56 ++ src/tests/spirit/fcnInput.txt | 2 +- 6 files changed, 603 insertions(+), 1 deletion(-) create mode 100644 src/tests/skewedGaussianTest/README create mode 100644 src/tests/skewedGaussianTest/fakeData.cpp create mode 100644 src/tests/skewedGaussianTest/fakeData.pro create mode 100644 src/tests/skewedGaussianTest/paramInput.dat create mode 100644 src/tests/skewedGaussianTest/skewedGaussian.C diff --git a/src/tests/skewedGaussianTest/README b/src/tests/skewedGaussianTest/README new file mode 100644 index 00000000..1df56fcd --- /dev/null +++ b/src/tests/skewedGaussianTest/README @@ -0,0 +1,17 @@ +------------------------------------------------------------------------- +Author: Andreas Suter, 2008/02/13 +$Id$ +------------------------------------------------------------------------- + +Generate fake data: + +fakeData.pro +fakeData.cpp + +fakeData is producing fake data root files which can be used in analysis +software like maxent and wkm (wkm only after converting to nemu files). + +------------------------------------------------------------------------- + +skewedGaussian.C is a simple root macro which is generating a skewed +Gaussian input file for fakeData diff --git a/src/tests/skewedGaussianTest/fakeData.cpp b/src/tests/skewedGaussianTest/fakeData.cpp new file mode 100644 index 00000000..afc250f1 --- /dev/null +++ b/src/tests/skewedGaussianTest/fakeData.cpp @@ -0,0 +1,481 @@ +/*************************************************************************** + + fakeData.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 +#include +#include +using namespace std; + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "TLemRunHeader.h" + +typedef vector PDoubleVector; +typedef vector PIntVector; + +void fakeDataSyntax() +{ + cout << endl << "usage: fakeData [--dump ]"; + cout << endl << " : file holding p(B) in columns B, pB, separated by ','"; + cout << endl << " comment lines start with a '#'"; + cout << endl << " : holding parameters needed to generate the histograms"; + cout << endl << " (see below)"; + cout << endl << " : output file name of the fake data histo file"; + cout << endl << " --dump : option to dump p(B) and A_i(t) in a separate root file"; + cout << endl; + cout << endl << " structure:"; + cout << endl << " comment lines start with a '#'"; + cout << endl << " time resolution : \"resolution\", followed by time resolution in (us)"; + cout << endl << " No of channles : \"channels\", followed by number of channels in a histogram"; + cout << endl << " t0 data line : \"t0\", followed by the t0's for each histogram, separated by ','"; + cout << endl << " asym data line : \"asym\", followed by the asymmetries for each histogram"; + cout << endl << " phases data line: \"phase\", followed by the phases for each histogram"; + cout << endl << " N0 data line : \"N0\", followed by the N0's for each histogram"; + cout << endl << " Bkg data line : \"bkg\", followed by the background for each histogram"; + cout << endl << " Example: "; + cout << endl << " # example"; + cout << endl << " resolution, 0.0001953125"; + cout << endl << " channels, 66601"; + cout << endl << " t0, 3300, 3300, 3300, 3300"; + cout << endl << " asym, 0.26, 0.26, 0.26, 0.26"; + cout << endl << " phase, 0.0, 90.0, 180.0, 270.0"; + cout << endl << " N0, 1000, 1100, 980, 1020"; + cout << endl << " bkg, 8, 11, 6, 15"; + cout << endl; +} + +int main(int argc, char *argv[]) +{ + const Double_t gamma_mu = TMath::TwoPi() * 0.1355; // in (rad/ns/T) + const Double_t tau_mu = 2197.147; // muon life time in ns + + if ((argc != 4) && (argc !=6)) { + fakeDataSyntax(); + return -1; + } + + TString pBFileName(argv[1]); + TString paramInputFileName(argv[2]); + TString rootOutputFileName(argv[3]); + TString dumpFileName(""); + if (argc == 6) { + dumpFileName = TString(argv[5]); + } + + vector B, pB; + + cout << endl << ">> read p(B) file ..."; + + // read pB input file and fill B and pB + ifstream f; + + f.open(pBFileName.Data(), iostream::in); + if (!f.is_open()) { + cout << endl << "**ERROR**: Sorry, couldn't open input file (p(B) file): " << pBFileName.Data(); + cout << endl << " Will quit"; + cout << endl; + return -1; + } + + char str[256]; + Double_t bb, ppb; + int status, lineNo=0; + while (!f.eof()) { + // read a line + f.getline(str, sizeof(str)); + + // ignore comments or empty lines + if ((str[0] == '#') || (strlen(str)==0)) + continue; + + // get data + status = sscanf(str, "%lf, %lf", &bb, &ppb); + if (status != 2) { + cout << endl << "**ERROR**: Problems while reading the input file (line no " << lineNo << "), will quit."; + cout << endl << " status = " << status; + cout << endl; + f.close(); + return -1; + } + + // fill B and pB vector + B.push_back(bb*1.0e-4); // G -> T + pB.push_back(ppb); + + lineNo++; + } + + f.close(); + + // normalize p(B) + Double_t sum = 0.0; + for (unsigned int i=0; i t0; + vector asym; + vector phase; + vector N0; + vector bkg; + + cout << endl << ">> read parameter input file ..."; + + // open parameter input file and extract the parameters + f.open(paramInputFileName.Data(), iostream::in); + if (!f.is_open()) { + cout << endl << "**ERROR**: Sorry, couldn't open parameter input file: " << paramInputFileName.Data(); + cout << endl << " Will quit"; + cout << endl; + return -1; + } + + TObjArray *tokens; + TObjString *ostr; + TString tstr; + Double_t dval; + Int_t ival; + + lineNo = 0; + while (!f.eof()) { + // read a line + f.getline(str, sizeof(str)); + lineNo++; + + // ignore comments or empty lines + if ((str[0] == '#') || (strlen(str)==0)) + continue; + + // tokenize string + tokens = TString(str).Tokenize(", \t"); + + // find to proper tag and fill the parameters + if (tokens->GetEntries() < 2) { + cout << endl << "**ERROR**: found '" << str << "'"; + cout << endl << " This is not compatible with the parameter input file, will quit"; + cout << endl; + return -1; + } + + ostr = dynamic_cast(tokens->At(0)); + tstr = ostr->GetString(); + tstr.ToLower(); + if (tstr.Contains("resolution")) { // handle time resolution + ostr = dynamic_cast(tokens->At(1)); + timeResolution = ostr->GetString().Atof()*1000.0; // us->ns + if (timeResolution <= 0.0) { + cout << endl << "**ERROR** time resolution must be > 0.0!"; + cout << endl; + return -1; + } + } else if (tstr.Contains("channels")) { // handle number of channels + ostr = dynamic_cast(tokens->At(1)); + noOfChannels = ostr->GetString().Atoi(); + if (noOfChannels <= 0.0) { + cout << endl << "**ERROR** number of channels must be > 0.0!"; + cout << endl; + return -1; + } + } else if (tstr.Contains("t0")) { // handle t0's + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tokens->At(i)); + ival = ostr->GetString().Atoi(); + t0.push_back(ival); + } + } else if (tstr.Contains("asym")) { // handle asymmetries + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tokens->At(i)); + dval = ostr->GetString().Atof(); + asym.push_back(dval); + } + } else if (tstr.Contains("phase")) { // handle phases + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tokens->At(i)); + dval = ostr->GetString().Atof(); + phase.push_back(dval/180.0*TMath::Pi()); + } + } else if (tstr.Contains("n0")) { // handle N0's + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tokens->At(i)); + ival = ostr->GetString().Atoi(); + N0.push_back(ival); + } + } else if (tstr.Contains("bkg")) { // handle bkg's + for (Int_t i=1; iGetEntries(); i++) { + ostr = dynamic_cast(tokens->At(i)); + ival = ostr->GetString().Atoi(); + bkg.push_back(ival); + } + } else { + cout << endl << "**ERROR** while reading the input file in line " << lineNo; + cout << endl; + return -1; + } + + // clean up + if (tokens) { + delete tokens; + tokens = 0; + } + } + + f.close(); + + cout << endl << ">> number of pB's = " << pB.size(); + cout << endl << ">> number of channels = " << noOfChannels; + + // check that all the parameter vectors have the same length + if ((t0.size() != asym.size()) || (t0.size() != phase.size()) || + (t0.size() != N0.size()) || (t0.size() != bkg.size())) { + cout << endl << "**ERROR** the number of elements of t0's etc. in the input parameter file must be the same!"; + cout << endl; + return -1; + } + + cout << endl << ">> calculate all the A_i(t) ..." << endl; + + PDoubleVector ddata; + vector asymmetry; + // calculate the asymmetry + for (UInt_t i=0; i> asymmetry " << i+1 << "/" << asym.size() << " done ..."; + } + + cout << endl << ">> calculate all the histo_i(t) ..." << endl; + + // calculate the histograms + PIntVector idata; + vector histo; + for (UInt_t i=0; i> histo " << i+1 << "/" << asym.size() << " done ..."; + } + + + if (dumpFileName.Length() > 0) { + // dump some test stuff + TFile tf(dumpFileName.Data(), "recreate"); + + TCanvas *c1 = new TCanvas("c1", "p(B)", 10, 10, 800, 600); + c1->Show(); + + // fill in the TH1F's TGraph's etc. + TGraph *graph_pB = new TGraph(B.size()); + for (UInt_t i=0; iSetPoint(i, B[i], pB[i]); + } + graph_pB->SetMarkerStyle(21); + graph_pB->DrawClone("apl"); + + tf.WriteTObject(c1); + + if (graph_pB) + delete graph_pB; + if (c1) + delete c1; + + TCanvas *chist; + TH1F *hh; + TString name, title; + for (UInt_t i=0; iShow(); + + name = "hh"; + name += i; + title = "histo"; + title += i; + hh = new TH1F(name.Data(), title.Data(), noOfChannels, + -timeResolution/2.0, (noOfChannels+0.5)*timeResolution); + for (Int_t j=0; jSetBinContent(j, asymmetry[i][j]); + } + hh->Draw("*H HIST"); + + tf.WriteTObject(chist); + + if (hh) + delete hh; + if (chist) + delete chist; + } + + tf.Close(); + } + + // add Poisson noise to the histograms + + cout << endl << ">> add Poisson noise ..." << endl; + + TH1F* theoHisto; + TH1F* fakeHisto; + vector histoData; + + TString name; + for (UInt_t i=0; iSetBinContent(j, histo[i][j]); + // fill fakeHisto + fakeHisto->FillRandom(theoHisto, (Int_t)theoHisto->Integral()); + + // keep fake data + histoData.push_back(fakeHisto); + + // cleanup + if (theoHisto) { + delete theoHisto; + theoHisto = 0; + } + } + + // save the histograms as root files + + // create run info folder and content + TFolder *runInfoFolder = new TFolder("RunInfo", "Run Info"); + TLemRunHeader *runHeader = new TLemRunHeader(); + sprintf(str, "Fake Data generated from %s", pBFileName.Data()); + runHeader->SetRunTitle(str); + Float_t fval = timeResolution; + runHeader->SetTimeResolution(fval); + runHeader->SetNChannels(noOfChannels); + runHeader->SetNHist(histoData.size()); + Int_t *t0array = new Int_t[histoData.size()]; + for (UInt_t i=0; iSetTimeZero(t0array); + if (t0array) + delete t0array; + runInfoFolder->Add(runHeader); + + // create decay histo folder and content + TFolder *histoFolder = new TFolder("histos", "histos"); + TFolder *decayAnaModule = new TFolder("DecayAnaModule", "DecayAnaModule"); + histoFolder->Add(decayAnaModule); + // no post pileup corrected (NPP) + for (UInt_t i=0; iAdd(histoData[i]); + // post pileup corrected (PPC) + vector histoDataPPC; + for (UInt_t i=0; i(histoData[i]->Clone())); + if (i < 10) + name = "hDecay2"; + else + name = "hDecay"; + name += i; + histoDataPPC[i]->SetNameTitle(name.Data(), name.Data()); + decayAnaModule->Add(histoDataPPC[i]); + } + + // write file + TFile fdf(rootOutputFileName.Data(), "recreate"); + runInfoFolder->Write("RunInfo", TObject::kSingleKey); + histoFolder->Write(); + fdf.Close(); + + // clean up + B.clear(); + pB.clear(); + t0.clear(); + asym.clear(); + phase.clear(); + N0.clear(); + bkg.clear(); + + for (UInt_t i=0; i> DONE." << endl; + + return 0; +} diff --git a/src/tests/skewedGaussianTest/fakeData.pro b/src/tests/skewedGaussianTest/fakeData.pro new file mode 100644 index 00000000..e36c4b1e --- /dev/null +++ b/src/tests/skewedGaussianTest/fakeData.pro @@ -0,0 +1,23 @@ +#------------------------------------------------------ +# fakeData.pro +# qmake file for fakeData +# +# Andreas Suter, 2008/02/12 +# +# $Id$ +# +#------------------------------------------------------ + +MAKEFILE = Makefile.fakeData + +CONFIG += warn_on debug + +SOURCES = fakeData.cpp \ + +INCLUDEPATH += $$(ROOTSYS)/include + +ROOTLIBS = -L$$(ROOTSYS)/lib $(shell $(ROOTSYS)/bin/root-config --libs) +NEMULIBS = -lTLemRunHeader +unix:LIBS += -L$$(ROOTSYS)/lib $${ROOTLIBS} $${NEMULIBS} -lbsd -lm -ldl -lutil + +TARGET = fakeData diff --git a/src/tests/skewedGaussianTest/paramInput.dat b/src/tests/skewedGaussianTest/paramInput.dat new file mode 100644 index 00000000..b8827c11 --- /dev/null +++ b/src/tests/skewedGaussianTest/paramInput.dat @@ -0,0 +1,25 @@ +# parameter input file to generate fake data histograms +#------------------------------------------------------ +# time resolution in (us) +resolution, 0.0001953125 +#------------------------------------------------------ +# number of channels +#channels, 66601 +channels, 66601 +#------------------------------------------------------ +# t0's +t0, 3300, 3300, 3300, 3300 +#------------------------------------------------------ +# asym's +asym, 0.26, 0.26, 0.26, 0.26 +#------------------------------------------------------ +# phases in (°) +phase, 0.0, 90.0, 180.0, 270.0 +#------------------------------------------------------ +# N0's +N0, 1100, 1000, 1020, 990 +#------------------------------------------------------ +# bkg's +bkg, 12, 10, 15, 8 +#------------------------------------------------------ +# end diff --git a/src/tests/skewedGaussianTest/skewedGaussian.C b/src/tests/skewedGaussianTest/skewedGaussian.C new file mode 100644 index 00000000..f6f3cced --- /dev/null +++ b/src/tests/skewedGaussianTest/skewedGaussian.C @@ -0,0 +1,56 @@ +/*************************************************************************** + + skewedGaussian.C + + Author: Andreas Suter + e-mail: andreas.suter@psi.ch + + $Id$ + +***************************************************************************/ + +//------------------------------------------------------------------------- +Double_t PofB(Double_t B, Double_t B0, Double_t sm, Double_t sp) +{ + Double_t val = 0.0; + + if (B < B0) { + val = TMath::Exp(-0.5*TMath::Power((B-B0)/sm, 2.0)); + } else { + val = TMath::Exp(-0.5*TMath::Power((B-B0)/sp, 2.0)); + } + + return val; +} + +//------------------------------------------------------------------------- +void skewedGaussian() +{ + FILE *fp; + char fln[256]; + + strcpy(fln, "skewedGauss.dat"); + + const Double_t B0 = 100.0; // (G) + const Double_t sm = 2.0; // (G) + const Double_t sp = 10.0; // (G) + + const Int_t noOfPoints = 1000; + + const Double_t res = 10.0*(sp+sm)/(Double_t)noOfPoints; + + fp = fopen(fln, "w"); + + fprintf(fp, "# title: skewed Gaussian\n"); + fprintf(fp, "# parameter: B0 = %lf, sm = %lf, sp = %lf\n", B0, sm, sp); + fprintf(fp, "# -----------------------------------------------------\n"); + fprintf(fp, "# B, pB\n"); + Double_t b, pb; + for (Int_t i=0; i