diff --git a/src/tests/skewedGaussianTest/README b/src/tests/skewedGaussianTest/README index 259e5704..ea6400fe 100644 --- a/src/tests/skewedGaussianTest/README +++ b/src/tests/skewedGaussianTest/README @@ -5,8 +5,11 @@ $Id$ Generate fake data: -fakeData.pro : qmake input file to generate the Makefile -fakeData.cpp : main +fakeData.pro : qmake input file to generate the Makefile +fakeData.cpp : main -> generates root-file output + +fakeDataNemu.pro : qmake input file to generate the Makefile +fakeDataNemu.cpp : main -> generates nemu-file output paramInput.dat : example input file for setting t0, asym, N0, bkg, time resolution, no of channels @@ -14,6 +17,9 @@ paramInput.dat : example input file for setting t0, asym, N0, bkg, 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). +fakeDataNemu is producing fake data nemu 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 diff --git a/src/tests/skewedGaussianTest/fakeDataNemu.cpp b/src/tests/skewedGaussianTest/fakeDataNemu.cpp new file mode 100644 index 00000000..42c127ec --- /dev/null +++ b/src/tests/skewedGaussianTest/fakeDataNemu.cpp @@ -0,0 +1,501 @@ +/*************************************************************************** + + fakeDataNemu.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 + +typedef vector PDoubleVector; +typedef vector PIntVector; + +void fakeDataNemuSyntax() +{ + cout << endl << "usage: fakeDataNemu [--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)) { + fakeDataNemuSyntax(); + return -1; + } + + TString pBFileName(argv[1]); + TString paramInputFileName(argv[2]); + TString nemuOutputFileName(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 fpB; + + fpB.open(pBFileName.Data(), iostream::in); + if (!fpB.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 (!fpB.eof()) { + // read a line + fpB.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; + fpB.close(); + return -1; + } + + // fill B and pB vector + B.push_back(bb*1.0e-4); // G -> T + pB.push_back(ppb); + + lineNo++; + } + + fpB.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 ..."; + + // read parameter input file + ifstream fparam; + + // open parameter input file and extract the parameters + fparam.open(paramInputFileName.Data(), iostream::in); + if (!fparam.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 (!fparam.eof()) { + // read a line + fparam.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; + } + } + + fparam.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 + 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 = "asym"; + 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; + } + + 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, histo[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 nemu file + ofstream fout(nemuOutputFileName.Data()); + if (!fout.is_open()) { + cout << endl << "**ERROR**: Sorry, couldn't open nemu output file: " << nemuOutputFileName.Data(); + cout << endl; + return -1; + } + + // write header + TTimeStamp *timeStamp = new TTimeStamp(); + fout << endl << "- created by fakeDataNemu (AS - 08/03/03) -"; + fout << endl << "NEMU_Run: " << nemuOutputFileName.Data(); + fout << endl << "nemu_Run: " << nemuOutputFileName.Data(); + fout << endl << "Date: " << timeStamp->AsString("s") << ", " << timeStamp->AsString("s"); + fout << endl << "Title: fakeData"; + fout << endl << "Field: 000000"; + fout << endl << "Setup: ??????"; + fout << endl << "Temp: 000000"; + fout << endl << "TOF(M3S1): nocut"; + fout << endl << "Groups: " << histoData.size(); + fout << endl << "Channels: " << noOfChannels; + fout << endl << "Resolution: " << timeResolution/1000.0; + // clean up + if (timeStamp) + delete timeStamp; + + // write data + for (UInt_t i=0; iGetNbinsX(); j++) { // loop over all data within a histo + if (j % 10 == 0) + fout << endl; + fout.width(8); + fout << (Int_t)histoData[i]->GetBinContent(j); + } + } + + fout.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/fakeDataNemu.pro b/src/tests/skewedGaussianTest/fakeDataNemu.pro new file mode 100644 index 00000000..e156397d --- /dev/null +++ b/src/tests/skewedGaussianTest/fakeDataNemu.pro @@ -0,0 +1,22 @@ +#------------------------------------------------------ +# fakeDataNemu.pro +# qmake file for fakeDataNemu +# +# Andreas Suter, 2008/02/12 +# +# $Id$ +# +#------------------------------------------------------ + +MAKEFILE = Makefile.fakeDataNemu + +CONFIG += warn_on debug + +SOURCES = fakeDataNemu.cpp \ + +INCLUDEPATH += $$(ROOTSYS)/include + +ROOTLIBS = -L$$(ROOTSYS)/lib $(shell $(ROOTSYS)/bin/root-config --libs) +unix:LIBS += -L$$(ROOTSYS)/lib $${ROOTLIBS} -lbsd -lm -ldl -lutil + +TARGET = fakeDataNemu