diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index 9fde611b..4c8961aa 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -1382,7 +1382,7 @@ void PMusrCanvas::HandleDataSet(unsigned int plotNo, unsigned int runNo, PRunDat end = data->fDataTimeStart + data->fValue.size()*data->fDataTimeStep + data->fDataTimeStep/2.0; // invoke histo - dataHisto = new TH1F(name, name, data->fValue.size()+2, start, end); + dataHisto = new TH1F(name, name, data->fValue.size()+1, start, end); // fill histogram for (unsigned int i=0; ifValue.size(); i++) { @@ -1422,7 +1422,7 @@ void PMusrCanvas::HandleDataSet(unsigned int plotNo, unsigned int runNo, PRunDat //cout << endl << ">> start = " << start << ", end = " << end << endl; // invoke histo - theoHisto = new TH1F(name, name, data->fTheory.size()+2, start, end); + theoHisto = new TH1F(name, name, data->fTheory.size()+1, start, end); // fill histogram for (unsigned int i=0; ifTheory.size(); i++) { diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index c5890f29..186f4d9f 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -699,7 +699,7 @@ bool PRunAsymmetry::PrepareFitData(PRawRunData* runData, unsigned int histoNo[2] double f, b, ef, eb; // fill data time start, and step // data start at data_start-t0 - fData.fDataTimeStart = fTimeResolution*(((double)start[0]-t0[0])+(double)fRunInfo->fPacking/2.0); + fData.fDataTimeStart = fTimeResolution*((double)start[0]-t0[0]); fData.fDataTimeStep = fTimeResolution*(double)fRunInfo->fPacking; for (unsigned int i=0; ifDataRange[0]-fRunInfo->fPacking*(fRunInfo->fDataRange[0]/fRunInfo->fPacking), - fRunInfo->fDataRange[2]-fRunInfo->fPacking*(fRunInfo->fDataRange[2]/fRunInfo->fPacking)}; + int val = fRunInfo->fDataRange[0]-fRunInfo->fPacking*(fRunInfo->fDataRange[0]/fRunInfo->fPacking); + do { + if (fRunInfo->fDataRange[2] - fRunInfo->fDataRange[0] < 0) + val += fRunInfo->fPacking; + } while (val + fRunInfo->fDataRange[2] - fRunInfo->fDataRange[0] < 0); + + int start[2] = {val, val + fRunInfo->fDataRange[2] - fRunInfo->fDataRange[0]}; int end[2]; double t0[2] = {fT0s[0], fT0s[1]}; +/* +cout << endl << ">> start[0]=" << start[0] << ", end[0]=" << end[0]; +cout << endl << ">> start[1]=" << start[1] << ", end[1]=" << end[1]; +cout << endl; +*/ + // make sure that there are equal number of rebinned bins in forward and backward unsigned int noOfBins0 = (runData->fDataBin[histoNo[0]].size()-start[0])/fRunInfo->fPacking; unsigned int noOfBins1 = (runData->fDataBin[histoNo[1]].size()-start[1])/fRunInfo->fPacking; @@ -806,34 +817,44 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2 double error = 0.0; // forward for (int i=start[0]; ifPacking == 0) && (i != start[0])) { // 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; - forwardPacked.fValue.push_back(value); - if (value == 0.0) - forwardPacked.fError.push_back(1.0); - else - forwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking); - value = 0.0; - error = 0.0; + if (fRunInfo->fPacking == 1) { + forwardPacked.fValue.push_back(fForward[i]); + forwardPacked.fError.push_back(fForwardErr[i]); + } else { // packed data, i.e. fRunInfo->fPacking > 1 + if (((i-start[0]) % fRunInfo->fPacking == 0) && (i != start[0])) { // 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; + forwardPacked.fValue.push_back(value); + if (value == 0.0) + forwardPacked.fError.push_back(1.0); + else + forwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking); + value = 0.0; + error = 0.0; + } } value += fForward[i]; error += fForwardErr[i]*fForwardErr[i]; } // backward for (int i=start[1]; ifPacking == 0) && (i != start[1])) { // 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; - backwardPacked.fValue.push_back(value); - if (value == 0.0) - backwardPacked.fError.push_back(1.0); - else - backwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking); - value = 0.0; - error = 0.0; + if (fRunInfo->fPacking == 1) { + backwardPacked.fValue.push_back(fBackward[i]); + backwardPacked.fError.push_back(fBackwardErr[i]); + } else { // packed data, i.e. fRunInfo->fPacking > 1 + if (((i-start[1]) % fRunInfo->fPacking == 0) && (i != start[1])) { // 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; + backwardPacked.fValue.push_back(value); + if (value == 0.0) + backwardPacked.fError.push_back(1.0); + else + backwardPacked.fError.push_back(TMath::Sqrt(error)/fRunInfo->fPacking); + value = 0.0; + error = 0.0; + } } value += fBackward[i]; error += fBackwardErr[i]*fBackwardErr[i]; @@ -851,9 +872,14 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2 double f, b, ef, eb, alpha = 1.0, beta = 1.0; // fill data time start, and step // data start at data_start-t0 - fData.fDataTimeStart = fTimeResolution*(((double)start[0]-t0[0])+(double)fRunInfo->fPacking/2.0); + fData.fDataTimeStart = fTimeResolution*((double)start[0]-t0[0]); fData.fDataTimeStep = fTimeResolution*(double)fRunInfo->fPacking; +/* +cout << endl << ">> start time = " << fData.fDataTimeStart << ", step = " << fData.fDataTimeStep; +cout << endl << "--------------------------------" << endl; +*/ + // get the proper alpha and beta switch (fAlphaBetaTag) { case 1: // alpha == 1, beta == 1 @@ -877,7 +903,7 @@ bool PRunAsymmetry::PrepareViewData(PRawRunData* runData, unsigned int histoNo[2 } //cout << endl << ">> alpha = " << alpha << ", beta = " << beta; - for (unsigned int i=0; ifDataBin[histoNo[0]].size(); - double startTime = -fT0s[0]*fTimeResolution; - fData.fTheoryTimeStart = startTime; - fData.fTheoryTimeStep = fTimeResolution; + double factor = 1.0; + if (fData.fValue.size() * 10 > runData->fDataBin[histoNo[0]].size()) { + size = fData.fValue.size() * 10; + factor = (double)runData->fDataBin[histoNo[0]].size() / (double)size; + } +//cout << endl << ">> runData->fDataBin[histoNo[0]].size() = " << runData->fDataBin[histoNo[0]].size() << ", fData.fValue.size() * 10 = " << fData.fValue.size() * 10 << ", size = " << size << ", factor = " << factor << endl; + fData.fTheoryTimeStart = fData.fDataTimeStart; + fData.fTheoryTimeStep = fTimeResolution*factor; for (unsigned int i=0; iFunc(time, par, fFuncValues); if (fabs(value) > 10.0) { // dirty hack needs to be fixed!! value = 0.0; diff --git a/src/classes/PRunDataHandler.cpp b/src/classes/PRunDataHandler.cpp index ce6dad26..2866e0cf 100644 --- a/src/classes/PRunDataHandler.cpp +++ b/src/classes/PRunDataHandler.cpp @@ -865,6 +865,7 @@ cout << endl; for (int j=0; jfPacking-1)/2.0); + fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0); fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking; for (int i=start; ifPacking == 1) { value = runData->fDataBin[histoNo][i]; + normalizer = fRunInfo->fPacking * (fTimeResolution * 1e3); // fTimeResolution us->ns + value /= normalizer; fData.fValue.push_back(value); if (value == 0.0) fData.fError.push_back(1.0); @@ -579,8 +581,9 @@ bool PRunSingleHisto::PrepareRawViewData(PRawRunData* runData, const unsigned in double value = 0.0; // 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-(double)t0+(double)(fRunInfo->fPacking-1)/2.0); + fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0); fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking; + /* cout << endl << ">> time resolution = " << fTimeResolution; cout << endl << ">> start = " << start << ", t0 = " << t0 << ", packing = " << fRunInfo->fPacking; @@ -662,11 +665,17 @@ cout << endl << ">> data start time = " << fData.fDataTimeStart; // calculate theory unsigned int size = runData->fDataBin[histoNo].size(); + double factor = 1.0; + if (fData.fValue.size() * 10 > runData->fDataBin[histoNo].size()) { + size = fData.fValue.size() * 10; + factor = (double)runData->fDataBin[histoNo].size() / (double)size; + } +//cout << endl << ">> runData->fDataBin[histoNo].size() = " << runData->fDataBin[histoNo].size() << ", fData.fValue.size() * 10 = " << fData.fValue.size() * 10 << ", size = " << size << ", factor = " << factor << endl; double theoryValue; fData.fTheoryTimeStart = fData.fDataTimeStart; - fData.fTheoryTimeStep = fTimeResolution; + fData.fTheoryTimeStep = fTimeResolution*factor; for (unsigned int i=0; iFunc(time, par, fFuncValues); if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!! theoryValue = 0.0; @@ -773,14 +782,15 @@ bool PRunSingleHisto::PrepareViewData(PRawRunData* runData, const unsigned int h 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-(double)t0+(double)(fRunInfo->fPacking-1)/2.0); + fData.fDataTimeStart = fTimeResolution*((double)start-(double)t0); fData.fDataTimeStep = fTimeResolution*fRunInfo->fPacking; + /* cout << endl << ">> start time = " << fData.fDataTimeStart << ", step = " << fData.fDataTimeStep; cout << endl << ">> start = " << start << ", end = " << end; -cout << endl << "--------------------------------"; +cout << endl << "--------------------------------" << endl; */ + double normalizer = 1.0; for (int i=start; ifPacking == 0) && (i != start)) { // fill data @@ -815,13 +825,19 @@ cout << endl << "--------------------------------"; } // calculate theory - unsigned int size = runData->fDataBin[histoNo].size(); double theoryValue; + unsigned int size = runData->fDataBin[histoNo].size(); + double factor = 1.0; + if (fData.fValue.size() * 10 > runData->fDataBin[histoNo].size()) { + size = fData.fValue.size() * 10; + factor = (double)runData->fDataBin[histoNo].size() / (double)size; + } +//cout << endl << ">> runData->fDataBin[histoNo].size() = " << runData->fDataBin[histoNo].size() << ", fData.fValue.size() * 10 = " << fData.fValue.size() * 10 << ", size = " << size << ", factor = " << factor << endl; fData.fTheoryTimeStart = fData.fDataTimeStart; - fData.fTheoryTimeStep = fTimeResolution; + fData.fTheoryTimeStep = fTimeResolution*factor; //cout << endl << ">> size=" << size << ", startTime=" << startTime << ", fTimeResolution=" << fTimeResolution; for (unsigned int i=0; iFunc(time, par, fFuncValues); if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!! theoryValue = 0.0; diff --git a/src/tests/skewedGaussianTest/fakeData.cpp b/src/tests/skewedGaussianTest/fakeData.cpp index 6d5c2bff..fd40df20 100644 --- a/src/tests/skewedGaussianTest/fakeData.cpp +++ b/src/tests/skewedGaussianTest/fakeData.cpp @@ -83,7 +83,8 @@ void fakeDataSyntax() 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 + // muon life time in (us), see PRL99, 032001 (2007) + const Double_t tau_mu = 2197.019; // muon life time in ns if ((argc != 4) && (argc !=6)) { fakeDataSyntax(); @@ -285,10 +286,11 @@ int main(int argc, char *argv[]) for (UInt_t i=0; iSetBinContent(j, asymmetry[i][j]); + hh->SetBinContent(j+1, asymmetry[i][j]); } hh->Draw("*H HIST"); @@ -380,10 +382,10 @@ int main(int argc, char *argv[]) name += i; title = "histo"; title += i; - hh = new TH1F(name.Data(), title.Data(), noOfChannels, - -timeResolution/2.0, (noOfChannels+0.5)*timeResolution); + hh = new TH1F(name.Data(), title.Data(), noOfChannels+1, + -timeResolution*0.5, (noOfChannels+0.5)*timeResolution); for (Int_t j=0; jSetBinContent(j, histo[i][j]); + hh->SetBinContent(j+1, histo[i][j]); } hh->Draw("*H HIST"); @@ -411,18 +413,18 @@ int main(int argc, char *argv[]) // create histos name = "theoHisto"; name += i; - theoHisto = new TH1F(name.Data(), name.Data(), noOfChannels, - -timeResolution/2.0, (noOfChannels+0.5)*timeResolution); + theoHisto = new TH1F(name.Data(), name.Data(), noOfChannels+1, + -timeResolution*0.5, (noOfChannels+0.5)*timeResolution); if (i < 10) name = "hDecay0"; else name = "hDecay"; name += i; - fakeHisto = new TH1F(name.Data(), name.Data(), noOfChannels, - -timeResolution/2.0, (noOfChannels+0.5)*timeResolution); + fakeHisto = new TH1F(name.Data(), name.Data(), noOfChannels+1, + -timeResolution*0.5, (noOfChannels+0.5)*timeResolution); // fill theoHisto for (Int_t j=0; jSetBinContent(j, histo[i][j]); + theoHisto->SetBinContent(j+1, histo[i][j]); // fill fakeHisto fakeHisto->FillRandom(theoHisto, (Int_t)theoHisto->Integral()); diff --git a/src/tests/skewedGaussianTest/fakeDataNemu.cpp b/src/tests/skewedGaussianTest/fakeDataNemu.cpp index 42c127ec..346f4784 100644 --- a/src/tests/skewedGaussianTest/fakeDataNemu.cpp +++ b/src/tests/skewedGaussianTest/fakeDataNemu.cpp @@ -82,7 +82,8 @@ void fakeDataNemuSyntax() 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 + // muon life time in (us), see PRL99, 032001 (2007) + const Double_t tau_mu = 2197.019; // muon life time in ns if ((argc != 4) && (argc !=6)) { fakeDataNemuSyntax(); diff --git a/src/tests/skewedGaussianTest/paramInput.dat b/src/tests/skewedGaussianTest/paramInput.dat index a7aeec38..9604ba1d 100644 --- a/src/tests/skewedGaussianTest/paramInput.dat +++ b/src/tests/skewedGaussianTest/paramInput.dat @@ -13,7 +13,7 @@ t0, 3300, 3300, 3300, 3300 asym, 0.24, 0.24, 0.24, 0.24 #------------------------------------------------------ # phases in (°) -phase, 10.0, 90.0, 170.0, 270.0 +phase, 0.0, 90.0, 180.0, 270.0 #------------------------------------------------------ # N0's N0, 1100, 1000, 1002, 990 diff --git a/src/tests/skewedGaussianTest/skewedGaussian.C b/src/tests/skewedGaussianTest/skewedGaussian.C index 8329af16..68048f53 100644 --- a/src/tests/skewedGaussianTest/skewedGaussian.C +++ b/src/tests/skewedGaussianTest/skewedGaussian.C @@ -41,18 +41,18 @@ void skewedGaussian() FILE *fp; char fln[256]; - const Double_t w = 0.8; // weight of the skewed Gaussian - const Double_t B0 = 142.0; // skewed Gaussian B0 (G) - const Double_t sm = 2.5; // skewed Gaussian sigma- (G) - const Double_t sp = 4.5; // skewed Gaussian sigma+ (G) + const Double_t w = 1.0; // weight of the skewed Gaussian + const Double_t B0 = 30.0; // skewed Gaussian B0 (G) + const Double_t sm = 2.0; // skewed Gaussian sigma- (G) + const Double_t sp = 2.0; // skewed Gaussian sigma+ (G) - const Double_t B0ext = 150.0; // external field Gaussian B0 (G) - const Double_t sext = 1.2; // external field Gaussian sigma (G) + const Double_t B0ext = 30.0; // external field Gaussian B0 (G) + const Double_t sext = 10; // external field Gaussian sigma (G) sprintf(fln, "skewedGauss-B%0.2lf-sm%0.2lf-sp%0.2lf-w%0.1lf-Bext%0.2lf-sext%0.2lf.dat", B0, sm, sp, w, B0ext, sext); - const Int_t noOfPoints = 2000; + const Int_t noOfPoints = 8000; const Double_t res = 0.1;