From 499684fc8bf10b8d909633c5cc1596d7988258c6 Mon Sep 17 00:00:00 2001 From: nemu Date: Wed, 16 Apr 2008 14:32:42 +0000 Subject: [PATCH] some more bits into the direction of musrview --- src/classes/PMsrHandler.cpp | 10 +- src/classes/PMusrCanvas.cpp | 39 +++++- src/classes/PRunSingleHisto.cpp | 203 +++++++++++++++++++++++++++----- src/classes/PStartupHandler.cpp | 26 +++- src/musrfit_startup.xml | 6 + 5 files changed, 246 insertions(+), 38 deletions(-) diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index c4eccec3..b7646b5b 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1846,6 +1846,13 @@ bool PMsrHandler::HandleStatisticEntry(PMsrLines &lines) return false; } + // check if chisq or max.log likelihood + fStatistic.fChisq = true; + for (unsigned int i=0; icd(); if (fData.size() > 0) { - fData[0].data->Draw("pe1"); + fData[0].data->Draw("pe"); // set time range Double_t xmin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmin; Double_t xmax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fTmax; fData[0].data->GetXaxis()->SetRangeUser(xmin, xmax); // check if it is necessary to set the y-axis range + Double_t ymin = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmin; + Double_t ymax = fMsrHandler->GetMsrPlotList()->at(fPlotNumber).fYmax; + if ((ymin != -999.0) && (ymax != -999.0)) { + fData[0].data->GetYaxis()->SetRangeUser(ymin, ymax); + } // plot all remaining data for (unsigned int i=1; iDraw("pe1same"); + fData[i].data->Draw("pesame"); } // plot all the theory for (unsigned int i=0; iSetHeader(tstr); - // get run plot info + // get/set run plot info + unsigned int runNo; + PMsrPlotStructure plotInfo = fMsrHandler->GetMsrPlotList()->at(fPlotNumber); + PMsrRunList runs = *fMsrHandler->GetMsrRunList(); + for (unsigned int i=0; iAddEntry(fData[i].data, tstr.Data(), "p"); + } fInfoPad->Draw(); fMainCanvas->cd(); @@ -709,7 +740,7 @@ void PMusrCanvas::HandleSingleHistoDataSet(unsigned int runNo, PRunData *data) start = data->fTheoryTimeStart; end = data->fTheoryTimeStart + (data->fTheory.size()-1)*data->fTheoryTimeStep; -cout << endl << ">> start = " << start << ", end = " << end << endl; +//cout << endl << ">> start = " << start << ", end = " << end << endl; // invoke histo theoHisto = new TH1F(name, name, data->fTheory.size(), start, end); diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 690e559d..a52e68bd 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -137,28 +137,7 @@ double PRunSingleHisto::CalcChiSquare(const std::vector& par) } } - -// static int counter = 0; -// TString fln=fRunInfo->fRunName+"_"+(Long_t)fRunInfo->fForwardHistoNo+"_data.dat"; -// ofstream f(fln.Data(),ios_base::out); -// for (unsigned int i=0; ifRunName+"_"+(Long_t)fRunInfo->fForwardHistoNo+"_theo.dat"; -// ofstream ft(fln.Data(),ios_base::out); -// for (unsigned int i=0; iFunc(time, par, fFuncValues))+bkg; -// } -// ft.close(); -// counter++; -// if (counter == 4) exit(0); - -//cout << endl << ">> " << chisq*fRunInfo->fPacking; - return chisq*fRunInfo->fPacking; + return chisq; } //-------------------------------------------------------------------------- @@ -407,7 +386,7 @@ bool PRunSingleHisto::PrepareFitData() if (value == 0.0) fData.fError.push_back(1.0); else - fData.fError.push_back(TMath::Sqrt(value)); + fData.fError.push_back(TMath::Sqrt(value/fRunInfo->fPacking)); value = 0.0; } value += runData->fDataBin[histoNo][i]; @@ -506,17 +485,17 @@ bool PRunSingleHisto::PrepareRawViewData() } // 2nd check if start is within proper bounds if ((start < 0) || (start > runData->fDataBin[histoNo].size())) { - cout << endl << "PRunSingleHisto::PrepareData(): **ERROR** start data bin doesn't make any sense!"; + cout << endl << "PRunSingleHisto::PrepareRawViewData(): **ERROR** start data bin doesn't make any sense!"; return false; } // 3rd check if end is within proper bounds if ((end < 0) || (end > runData->fDataBin[histoNo].size())) { - cout << endl << "PRunSingleHisto::PrepareData(): **ERROR** end data bin doesn't make any sense!"; + cout << endl << "PRunSingleHisto::PrepareRawViewData(): **ERROR** end data bin doesn't make any sense!"; return false; } // 4th check if t0 is within proper bounds if ((t0 < 0) || (t0 > runData->fDataBin[histoNo].size())) { - cout << endl << "PRunSingleHisto::PrepareData(): **ERROR** t0 data bin doesn't make any sense!"; + cout << endl << "PRunSingleHisto::PrepareRawViewData(): **ERROR** t0 data bin doesn't make any sense!"; return false; } @@ -535,7 +514,7 @@ bool PRunSingleHisto::PrepareRawViewData() if (value == 0.0) fData.fError.push_back(1.0); else - fData.fError.push_back(TMath::Sqrt(value)); + fData.fError.push_back(TMath::Sqrt(value/fRunInfo->fPacking)); value = 0.0; } value += runData->fDataBin[histoNo][i]; @@ -595,10 +574,176 @@ bool PRunSingleHisto::PrepareRawViewData() // PrepareViewData //-------------------------------------------------------------------------- /** - *

Muon life time corrected data - * + *

Muon life time corrected data: Starting from + * \f[ N(t) = N_0 e^{-t/\tau} [ 1 + A(t) ] + \mathrm{Bkg} \f] + * it follows that + * \f[ A(t) = (-1) + e^{+t/\tau}\, \frac{N(t)-\mathrm{Bkg}}{N_0}. \f] + * For the error estimate only the statistical error of \f$ N(t) \f$ is used, and hence + * \f[ \Delta A(t) = \frac{e^{t/\tau}}{N_0}\,\sqrt{\frac{N(t)}{p}} \f] + * where \f$ p \f$ is the packing, and \f$ N(t) \f$ are the packed data, i.e. + * \f[ N(t_i) = \frac{1}{p}\, \sum_{j=i}^{i+p} n_j \f] + * with \f$ n_j \f$ the raw histogram data bins. */ bool PRunSingleHisto::PrepareViewData() { + // get the proper run + PRawRunData* runData = fRawData->GetRunData(fRunInfo->fRunName); + if (!runData) { // couldn't get run + cout << endl << "PRunSingleHisto::PrepareViewData(): **ERROR** Couldn't get run " << fRunInfo->fRunName.Data() << "!"; + return false; + } + + // keep the time resolution in (us) + fTimeResolution = runData->fTimeResolution/1.0e3; + + // check if the t0's are given in the msr-file + if (fRunInfo->fT0[0] == -1) { // t0's are NOT in the msr-file + // check if the t0's are in the data file + if (runData->fT0s.size() != 0) { // t0's in the run data + // keep the proper t0's. For single histo runs, forward is holding the histo no + // fForwardHistoNo starts with 1 not with 0 ;-) + fT0s.push_back(runData->fT0s[fRunInfo->fForwardHistoNo-1]); + } else { // t0's are neither in the run data nor in the msr-file -> not acceptable! + cout << endl << "PRunSingleHisto::PrepareViewData(): **ERROR** NO t0's found, neither in the run data nor in the msr-file!"; + return false; + } + } else { // t0's in the msr-file + // check if t0's are given in the data file + if (runData->fT0s.size() != 0) { + // compare t0's of the msr-file with the one in the data file + if (fabs(fRunInfo->fT0[0]-runData->fT0s[fRunInfo->fForwardHistoNo-1])>5.0) { // given in bins!! + cout << endl << "PRunSingleHisto::PrepareViewData(): **WARNING**:"; + cout << endl << " t0 from the msr-file is " << fRunInfo->fT0[0]; + cout << endl << " t0 from the data file is " << runData->fT0s[fRunInfo->fForwardHistoNo-1]; + cout << endl << " This is quite a deviation! Is this done intentionally??"; + cout << endl; + } + } + fT0s.push_back(fRunInfo->fT0[0]); + } + + // check if post pile up data shall be used + unsigned int histoNo; + if (fRunInfo->fFileFormat.Contains("ppc")) { + histoNo = runData->fDataBin.size()/2 + fRunInfo->fForwardHistoNo-1; + } else { + histoNo = fRunInfo->fForwardHistoNo-1; + } + + if ((runData->fDataBin.size() < histoNo) || (histoNo < 0)) { + cout << endl << "PRunSingleHisto::PrepareViewData(): **PANIC ERROR**:"; + cout << endl << " histoNo found = " << histoNo << ", but there are only " << runData->fDataBin.size() << " runs!?!?"; + cout << endl << " Will quite :-("; + cout << endl; + return false; + } + + // transform raw histo data. This is done the following way (for details see the manual): + // for the single histo fit, just the rebinned raw data are copied + // first get start data, end data, and t0 + unsigned int start; + unsigned int end; + double t0 = fT0s[0]; + // raw data, since PMusrCanvas is doing ranging etc. + // start = the first bin which is a multiple of packing backward from t0 + start = (int)t0 - ((int)t0/fRunInfo->fPacking)*fRunInfo->fPacking; + // end = last bin starting from start which is a multipl of packing and still within the data + end = start + ((runData->fDataBin[histoNo].size()-start-1)/fRunInfo->fPacking)*fRunInfo->fPacking; + // check if start, end, and t0 make any sense + // 1st check if start and end are in proper order + if (end < start) { // need to swap them + int keep = end; + end = start; + start = keep; + } + // 2nd check if start is within proper bounds + if ((start < 0) || (start > runData->fDataBin[histoNo].size())) { + cout << endl << "PRunSingleHisto::PrepareViewData(): **ERROR** start data bin doesn't make any sense!"; + return false; + } + // 3rd check if end is within proper bounds + if ((end < 0) || (end > runData->fDataBin[histoNo].size())) { + cout << endl << "PRunSingleHisto::PrepareViewData(): **ERROR** end data bin doesn't make any sense!"; + return false; + } + // 4th check if t0 is within proper bounds + if ((t0 < 0) || (t0 > runData->fDataBin[histoNo].size())) { + cout << endl << "PRunSingleHisto::PrepareViewData(): **ERROR** t0 data bin doesn't make any sense!"; + return false; + } + + // everything looks fine, hence fill data set + + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (unsigned int i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // calculate asymmetry + double N0 = par[fRunInfo->fNormParamNo-1]; +//cout << endl << ">> N0 = " << N0; + + // get tau + double tau; + if (fRunInfo->fLifetimeParamNo != -1) + tau = par[fRunInfo->fLifetimeParamNo-1]; + else + tau = PMUON_LIFETIME; +//cout << endl << ">> tau = " << tau; + + // get background + double bkg = par[fRunInfo->fBkgFitParamNo-1]; +//cout << endl << ">> bkg = " << bkg; + + double value = 0.0; + double expval; + double time; + // data start at data_start-t0 + // time shifted so that packing is included correctly, i.e. t0 == t0 after packing + fData.fDataTimeStart = fTimeResolution*(((double)start-t0)+(double)fRunInfo->fPacking/2.0); + fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking; +//cout << endl << ">> start = " << fData.fDataTimeStart << ", step = " << fData.fDataTimeStep; + for (unsigned i=start; ifPacking == 0) && (i != start)) { // fill data + // in order that after rebinning the fit does not need to be redone (important for plots) + // the value is normalize to per bin + value /= fRunInfo->fPacking; + time = fData.fDataTimeStart+(double)i*fData.fDataTimeStep/fRunInfo->fPacking; + expval = TMath::Exp(+time/tau)/N0; + fData.fValue.push_back(-1.0+expval*(value-bkg)); + fData.fError.push_back(expval*TMath::Sqrt(value/fRunInfo->fPacking)); +//cout << endl << ">> " << time << ", " << expval << ", " << -1.0+expval*(value-bkg) << ", " << expval*TMath::Sqrt(value/fRunInfo->fPacking); + value = 0.0; + } + value += runData->fDataBin[histoNo][i]; + } + + // count the number of bins to be fitted + fNoOfFitBins=0; + for (unsigned int i=0; i= fFitStartTime) && (time <= fFitStopTime)) + fNoOfFitBins++; + } + + // calculate functions + for (int i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), fRunInfo->fMap, par); + } + + // calculate theory + unsigned int size = runData->fDataBin[histoNo].size(); + double startTime = -fT0s[0]*fTimeResolution; + fData.fTheoryTimeStart = startTime; + fData.fTheoryTimeStep = fTimeResolution; + for (unsigned int i=0; iFunc(time, par, fFuncValues)); + } + + // clean up + par.clear(); + return true; } diff --git a/src/classes/PStartupHandler.cpp b/src/classes/PStartupHandler.cpp index d722d882..3f6f1c5f 100644 --- a/src/classes/PStartupHandler.cpp +++ b/src/classes/PStartupHandler.cpp @@ -93,12 +93,34 @@ void PStartupHandler::OnEndDocument() // check if any markers are given if (fMarkerList.size() == 0) { - // to be done yet + fMarkerList.push_back(24); // open circle + fMarkerList.push_back(25); // open square + fMarkerList.push_back(26); // open triangle + fMarkerList.push_back(27); // open diamond + fMarkerList.push_back(28); // open cross + fMarkerList.push_back(29); // full star + fMarkerList.push_back(30); // open star + fMarkerList.push_back(20); // full circle + fMarkerList.push_back(21); // full square + fMarkerList.push_back(22); // full triangle + fMarkerList.push_back(23); // full down triangle + fMarkerList.push_back(2); // thin cross + fMarkerList.push_back(3); // thin star + fMarkerList.push_back(5); // thin cross 45° rotated } // check if any colors are given if (fColorList.size() == 0) { - // to be done yet + fColorList.push_back(TColor::GetColor(0, 0, 0)); // kBlack + fColorList.push_back(TColor::GetColor(255, 0, 0)); // kRed + fColorList.push_back(TColor::GetColor(0, 255, 0)); // kGreen + fColorList.push_back(TColor::GetColor(0, 0, 255)); // kBlue + fColorList.push_back(TColor::GetColor(255, 0, 255)); // kMagneta + fColorList.push_back(TColor::GetColor(0, 255, 255)); // kCyan + fColorList.push_back(TColor::GetColor(156, 0, 255)); // kViolette-3 + fColorList.push_back(TColor::GetColor(99, 101, 49)); // kYellow-1 + fColorList.push_back(TColor::GetColor(49, 101, 49)); // kGreen-1 + fColorList.push_back(TColor::GetColor(156, 48, 0)); // kOrange-4 } } diff --git a/src/musrfit_startup.xml b/src/musrfit_startup.xml index bf3b8264..8d83b4e7 100644 --- a/src/musrfit_startup.xml +++ b/src/musrfit_startup.xml @@ -30,6 +30,12 @@ 255,0,0 0,255,0 0,0,255 + 255,0,255 + 0,255,255 + 156,0,255 + 99,101,49 + 49,101,49 + 156,48,0 \ No newline at end of file