LMU/musrfit
LMU/musrfit
5
0
Fork 0
Code Issues 6 Pull Requests Actions Packages Projects Releases 3 Activity

added a NPP, PPC switch and performed some clean-up

Browse Source
This commit is contained in:
nemu 2008-09-22 06:24:52 +00:00
parent a3756c511b
commit 8ba07a692a
11 changed files with 90 additions and 67 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

12
src/classes/PFitter.cpp
View File

@ -154,9 +154,11 @@ bool PFitter::DoFit()
break;
case PMN_MINOS:
status = ExecuteMinos();
// set positive errors true
for (unsigned int i=0; i<fParams.size(); i++) {
fRunInfo->SetMsrParamPosErrorPresent(i, true);
// set positive errors true if minos has been successfull
if (status) {
for (unsigned int i=0; i<fParams.size(); i++) {
fRunInfo->SetMsrParamPosErrorPresent(i, true);
}
}
break;
case PMN_PLOT:
@ -366,10 +368,12 @@ bool PFitter::ExecuteMigrad()
ROOT::Minuit2::FunctionMinimum min = migrad(maxfcn, tolerance);
if (!min.IsValid()) {
cout << endl << "**WARNING**: PFitter::ExecuteMigrad(): Fit did not converge, sorry ...";
/*
// set flag positive error present to false
for (unsigned int i=0; i<fParams.size(); i++) {
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
*/
return false;
}
@ -435,10 +439,12 @@ bool PFitter::ExecuteMinimize()
ROOT::Minuit2::FunctionMinimum min = minimize(maxfcn, tolerance);
if (!min.IsValid()) {
cout << endl << "**WARNING**: PFitter::ExecuteMinimize(): Fit did not converge, sorry ...";
/*
// set flag positive error present to false
for (unsigned int i=0; i<fParams.size(); i++) {
fRunInfo->SetMsrParamPosErrorPresent(i, false);
}
*/
return false;
}

4
src/classes/PFunction.cpp
View File

@ -59,7 +59,7 @@ using namespace std;
*/
PFunction::PFunction(tree_parse_info<> info)
{
cout << endl << "in PFunction ...";
// cout << endl << "in PFunction ...";
fInfo = info;
@ -528,7 +528,7 @@ void PFunction::EvalTreeForStringExpression(iter_t const& i)
} else if (i->value.id() == PFunctionGrammar::termID) {
termOp++;
if (*i->value.begin() == '*') {
cout << endl << ">> i->children.size() = " << i->children.size() << endl;
//cout << endl << ">> i->children.size() = " << i->children.size() << endl;
assert(i->children.size() == 2);
EvalTreeForStringExpression(i->children.begin());
fFuncString += " * ";

18
src/classes/PFunctionHandler.cpp
View File

@ -47,7 +47,7 @@ PFunctionHandler::PFunctionHandler(PMsrLines lines)
fValid = true;
fLines = lines;
cout << endl << "in PFunctionHandler(PMsrLines lines)";
// cout << endl << "in PFunctionHandler(PMsrLines lines)";
}
//-------------------------------------------------------------
@ -59,7 +59,7 @@ PFunctionHandler::PFunctionHandler(PMsrLines lines)
*/
PFunctionHandler::~PFunctionHandler()
{
cout << endl << "in ~PFunctionHandler()" << endl << endl;
// cout << endl << "in ~PFunctionHandler()" << endl << endl;
fLines.clear();
fFuncs.clear();
@ -74,7 +74,7 @@ PFunctionHandler::~PFunctionHandler()
*/
bool PFunctionHandler::DoParse()
{
cout << endl << "in PFunctionHandler::DoParse() ...";
// cout << endl << "in PFunctionHandler::DoParse() ...";
bool success = true;
PFunctionGrammar function;
@ -82,7 +82,7 @@ bool PFunctionHandler::DoParse()
// feed the function block into the parser. Start with i=1, since i=0 is FUNCTIONS
for (unsigned int i=1; i<fLines.size(); i++) {
cout << endl << "fLines[" << i << "] = '" << fLines[i].fLine.Data() << "'";
// cout << endl << "fLines[" << i << "] = '" << fLines[i].fLine.Data() << "'";
// function line to upper case
line = fLines[i].fLine;
@ -92,7 +92,7 @@ cout << endl << "fLines[" << i << "] = '" << fLines[i].fLine.Data() << "'";
tree_parse_info<> info = ast_parse(line.Data(), function, space_p);
if (info.full) {
cout << endl << "parse successful ..." << endl;
// cout << endl << "parse successful ..." << endl;
PFunction func(info);
fFuncs.push_back(func);
} else {
@ -115,10 +115,10 @@ cout << endl << "fLines[" << i << "] = '" << fLines[i].fLine.Data() << "'";
}
}
if (success) {
for (unsigned int i=0; i<fFuncs.size(); i++)
cout << endl << "func number = " << fFuncs[i].GetFuncNo();
}
// if (success) {
// for (unsigned int i=0; i<fFuncs.size(); i++)
// cout << endl << "func number = " << fFuncs[i].GetFuncNo();
// }
return success;
}

1
src/classes/PMusrCanvas.cpp
View File

@ -612,6 +612,7 @@ void PMusrCanvas::UpdateInfoPad()
// energy if present
tstr += TString("E=");
dval = fRunList->GetEnergy(runs[runNo].fRunName);
cout << endl << ">> dval = " << dval << " (Engery)";
if (dval == -9.9e99) {
tstr += TString("??,");
} else {

72
src/classes/PRunDataHandler.cpp
View File

@ -161,9 +161,9 @@ bool PRunDataHandler::ReadFile()
return false;
// everything looks fine, hence try to read the data file
if (!run_it->fFileFormat.CompareTo("root-npp")) // not post pile up corrected histos
success = ReadRootFile();
success = ReadRootFile(true);
else if (!run_it->fFileFormat.CompareTo("root-ppc")) // post pile up corrected histos
success = ReadRootFile();
success = ReadRootFile(false);
else if (!run_it->fFileFormat.CompareTo("nexus"))
success = ReadNexusFile();
else if (!run_it->fFileFormat.CompareTo("psi-bin"))
@ -342,12 +342,15 @@ cout << endl << ">> generated path: " << str.Data() << endl;
/**
* <p>
*
* \param notPostPileup This flag is used as a switch between "Not Post Pileup Corrected"
* and "Post Pileup Corrected" histogramms.
*/
bool PRunDataHandler::ReadRootFile()
bool PRunDataHandler::ReadRootFile(bool notPostPileup)
{
PDoubleVector histoData;
PRawRunData runData;
//cout << endl << ">> in ReadRootFile() ...";
TFile f(fRunPathName.Data());
if (f.IsZombie()) {
return false;
@ -379,6 +382,7 @@ bool PRunDataHandler::ReadRootFile()
// get implantation energy
runData.fEnergy = runHeader->GetImpEnergy();
//cout << endl << ">> runData.fEnergy = " << runData.fEnergy;
// get setup
runData.fSetup = runHeader->GetLemSetup().GetString();
@ -408,38 +412,42 @@ bool PRunDataHandler::ReadRootFile()
}
// get all the data
char histoName[32];
// read first the data which are NOT post pileup corrected
for (int i=0; i<noOfHistos; i++) {
sprintf(histoName, "hDecay%02d", i);
TH1F *histo = dynamic_cast<TH1F*>(folder->FindObjectAny(histoName));
if (!histo) {
cout << endl << "PRunDataHandler::ReadRootFile: Couldn't get histo " << histoName;
return false;
if (notPostPileup) { // read the data which are NOT post pileup corrected
for (int i=0; i<noOfHistos; i++) {
sprintf(histoName, "hDecay%02d", i);
//cout << endl << ">> histoName = " << histoName;
TH1F *histo = dynamic_cast<TH1F*>(folder->FindObjectAny(histoName));
if (!histo) {
cout << endl << "PRunDataHandler::ReadRootFile: Couldn't get histo " << histoName;
return false;
}
// fill data
for (int j=1; j<histo->GetNbinsX(); j++)
histoData.push_back(histo->GetBinContent(j));
// store them in runData vector
runData.fDataBin.push_back(histoData);
// clear histoData for the next histo
histoData.clear();
}
// fill data
for (int j=1; j<histo->GetNbinsX(); j++)
histoData.push_back(histo->GetBinContent(j));
// store them in runData vector
runData.fDataBin.push_back(histoData);
// clear histoData for the next histo
histoData.clear();
}
// now read the data which ARE post pileup corrected
for (int i=0; i<noOfHistos; i++) {
sprintf(histoName, "hDecay%02d", i+POST_PILEUP_HISTO_OFFSET);
TH1F *histo = dynamic_cast<TH1F*>(folder->FindObjectAny(histoName));
if (!histo) {
cout << endl << "PRunDataHandler::ReadRootFile: Couldn't get histo " << histoName;
return false;
} else { // read the data which ARE post pileup corrected
for (int i=0; i<noOfHistos; i++) {
sprintf(histoName, "hDecay%02d", i+POST_PILEUP_HISTO_OFFSET);
//cout << endl << ">> histoName = " << histoName;
TH1F *histo = dynamic_cast<TH1F*>(folder->FindObjectAny(histoName));
if (!histo) {
cout << endl << "PRunDataHandler::ReadRootFile: Couldn't get histo " << histoName;
return false;
}
// fill data
for (int j=1; j<histo->GetNbinsX(); j++)
histoData.push_back(histo->GetBinContent(j));
// store them in runData vector
runData.fDataBin.push_back(histoData);
// clear histoData for the next histo
histoData.clear();
}
// fill data
for (int j=1; j<histo->GetNbinsX(); j++)
histoData.push_back(histo->GetBinContent(j));
// store them in runData vector
runData.fDataBin.push_back(histoData);
// clear histoData for the next histo
histoData.clear();
}
//cout << endl;
f.Close();

8
src/classes/PRunListCollection.cpp
View File

@ -92,11 +92,11 @@ bool PRunListCollection::Add(int runNo, EPMusrHandleTag tag)
{
bool success = true;
PMsrRunStructure *runList = &(*fMsrInfo->GetMsrRunList())[runNo];
// PMsrRunStructure *runList = &(*fMsrInfo->GetMsrRunList())[runNo];
cout << endl << "PRunListCollection::Add(): will add run no " << runNo;
cout << ", name = " << runList->fRunName.Data();
cout << ", type = " << runList->fFitType;
// cout << endl << "PRunListCollection::Add(): will add run no " << runNo;
// cout << ", name = " << runList->fRunName.Data();
// cout << ", type = " << runList->fFitType;
int fitType = (*fMsrInfo->GetMsrRunList())[runNo].fFitType;

15
src/classes/PTheory.cpp
View File

@ -185,7 +185,7 @@ PTheory::PTheory(PMsrHandler *msrInfo, unsigned int runNo, const bool hasParent)
// if userFcn, the first entry is the function name and needs to be handled specially
if ((fType == THEORY_USER_FCN) && ((i == 1) || (i == 2))) {
cout << endl << ">> userFcn: i=" << i << ", str=" << str.Data() << endl;
//cout << endl << ">> userFcn: i=" << i << ", str=" << str.Data() << endl;
if (i == 1) {
fUserFcnSharedLibName = str;
}
@ -282,7 +282,7 @@ cout << endl << ">> userFcn: i=" << i << ", str=" << str.Data() << endl;
// invoke user function object
fUserFcn = 0;
fUserFcn = (PUserFcnBase*)TClass::GetClass(fUserFcnClassName.Data())->New();
cout << endl << ">> fUserFcn = " << fUserFcn << endl;
//cout << endl << ">> fUserFcn = " << fUserFcn << endl;
if (fUserFcn == 0) {
cout << endl << "**ERROR**: PTheory: user function object could not be invoked. See line no " << line->fLineNo;
fValid = false;
@ -690,7 +690,7 @@ void PTheory::MakeCleanAndTidyTheoryBlock(PMsrLines *fullTheoryBlock)
str.Resize(index);
// tokenize line
tokens = str.Tokenize(" \t");
cout << endl << ">> #tokens=" << tokens->GetEntries() << ", str=" << str.Data();
//cout << endl << ">> #tokens=" << tokens->GetEntries() << ", str=" << str.Data();
// make a handable string out of the asymmetry token
ostr = dynamic_cast<TObjString*>(tokens->At(0));
str = ostr->GetString();
@ -765,7 +765,7 @@ void PTheory::MakeCleanAndTidyPolynom(unsigned int i, PMsrLines *fullTheoryBlock
TObjString *ostr;
char substr[256];
cout << endl << ">> MakeCleanAndTidyPolynom: " << (*fullTheoryBlock)[i].fLine.Data();
//cout << endl << ">> MakeCleanAndTidyPolynom: " << (*fullTheoryBlock)[i].fLine.Data();
// init tidy
tidy = TString("polynom ");
@ -811,7 +811,7 @@ void PTheory::MakeCleanAndTidyUserFcn(unsigned int i, PMsrLines *fullTheoryBlock
TObjArray *tokens = 0;
TObjString *ostr;
cout << endl << ">> MakeCleanAndTidyUserFcn: " << (*fullTheoryBlock)[i].fLine.Data();
//cout << endl << ">> MakeCleanAndTidyUserFcn: " << (*fullTheoryBlock)[i].fLine.Data();
// init tidy
tidy = TString("userFcn ");
@ -1448,6 +1448,7 @@ double PTheory::Polynom(register double t, const PDoubleVector& paramValues, con
*/
double PTheory::UserFcn(register double t, const PDoubleVector& paramValues, const PDoubleVector& funcValues) const
{
/*
static bool first = true;
if (first) {
cout << endl << ">> UserFcn: fParamNo.size()=" << fParamNo.size() << ", fUserParam.size()=" << fUserParam.size();
@ -1456,6 +1457,7 @@ if (first) {
}
cout << endl;
}
*/
// check if FUNCTIONS are used
for (unsigned int i=0; i<fUserParam.size(); i++) {
@ -1466,6 +1468,7 @@ if (first) {
}
}
/*
if (first) {
first = false;
for (unsigned int i=0; i<fUserParam.size(); i++) {
@ -1475,6 +1478,6 @@ if (first) {
cout << endl << ">> t=" << t << ", function value=" << (*fUserFcn)(t, fUserParam);
cout << endl;
}
*/
return (*fUserFcn)(t, fUserParam);
}

2
src/include/PRunDataHandler.h
View File

@ -62,7 +62,7 @@ class PRunDataHandler
virtual bool ReadFile();
virtual bool FileAlreadyRead(PMsrRunStructure &runInfo);
virtual bool FileExistsCheck(PMsrRunStructure &runInfo);
virtual bool ReadRootFile();
virtual bool ReadRootFile(bool notPostPileup);
virtual bool ReadNexusFile();
virtual bool ReadNemuFile();
virtual bool ReadPsiBinFile();

1
src/tests/fourier/PMusrFourier.cpp
View File

@ -505,6 +505,7 @@ void PMusrFourier::ApodizeData(int apodizationTag) {
switch (apodizationTag) {
case F_APODIZATION_NONE:
return;
break;
case F_APODIZATION_WEAK:
c[0] = cweak[0]+cweak[1]+cweak[2];

2
src/tests/fourier/fourier.cpp
View File

@ -185,7 +185,7 @@ cout << endl << "#bins=" << histo->GetNbinsX();
cout << endl << ">> Do you wish to apodize your data (y/n)? ";
cin >> answer;
unsigned int apodizationTag=0;
unsigned int apodizationTag = F_APODIZATION_NONE;
if (strstr(answer, "y")) {
cout << endl << ">> apodization (1=weak, 2=medium, 3=strong, 4=user) = ";
cin >> apodizationTag;

22
src/tests/rootMinuit2Test/minuit2test.C
View File

@ -78,17 +78,17 @@ void minuit2test()
// set parameters
gFitFcn->SetParNames("N0", "asym", "lambda", "B", "phase", "Bkg");
gFitFcn->SetParameter(0, 30.0); // N0
// gFitFcn->SetParLimits(0, 0.0, 1.0e6);
gFitFcn->SetParLimits(0, 0.0, 1.0e6);
gFitFcn->SetParameter(1, 0.26); // asym
// gFitFcn->SetParLimits(1, 0.0, 0.33);
gFitFcn->SetParLimits(1, 0.0, 0.33);
gFitFcn->SetParameter(2, 0.3); // lambda
// gFitFcn->SetParLimits(2, 0.0, 100.0);
gFitFcn->SetParLimits(2, 0.0, 100.0);
gFitFcn->SetParameter(3, 100.0); // B
// gFitFcn->SetParLimits(3, 0.0, 1000.0);
gFitFcn->SetParLimits(3, 0.0, 1000.0);
gFitFcn->SetParameter(4, 0.0); // phase
// gFitFcn->SetParLimits(4, -90.0, 90.0);
gFitFcn->SetParLimits(4, -90.0, 90.0);
gFitFcn->SetParameter(5, 5.0); // Bkg
// gFitFcn->SetParLimits(5, 0.0, 1000.0);
gFitFcn->SetParLimits(5, 0.0, 1000.0);
cout << endl << "gFitFcn->Integral(0.0, 12.0) = " << gFitFcn->Integral(0.0, 12.0);
cout << endl;
@ -100,9 +100,13 @@ void minuit2test()
histo->Draw();
gFitFcn->SetParameter(0, 150.0); // N0
gFitFcn->SetParameter(5, 23.0); // Bkg
gFitFcn->SetParameter(0, 1000.0); // N0
gFitFcn->SetParameter(1, 0.1); // asym
gFitFcn->SetParameter(2, 1.0); // lambda
gFitFcn->SetParameter(3, 50.0); // B
gFitFcn->SetParameter(4, 0.0); // phase
gFitFcn->SetParameter(5, 300.0); // Bkg
TVirtualFitter::SetDefaultFitter("Minuit2");
histo->Fit("gFitFcn", "L"); // L->likleyhood, E->minos
histo->Fit("gFitFcn", ""); // L->likleyhood, E->minos
}