From 71cac92d235ced77b7448565b7402bef137bae52 Mon Sep 17 00:00:00 2001 From: Suter Andreas Date: Mon, 18 Jan 2016 10:41:27 +0100 Subject: [PATCH] more work toward an asymmetry RRF fit. --- src/classes/PFitterFcn.cpp | 2 + src/classes/PMsrHandler.cpp | 107 +++++++++++++++- src/classes/PMusr.cpp | 2 +- src/classes/PMusrCanvas.cpp | 1 + src/classes/PRunAsymmetryRRF.cpp | 211 ++++++++++++------------------- 5 files changed, 191 insertions(+), 132 deletions(-) diff --git a/src/classes/PFitterFcn.cpp b/src/classes/PFitterFcn.cpp index 5e5b23c6..c6a8566a 100644 --- a/src/classes/PFitterFcn.cpp +++ b/src/classes/PFitterFcn.cpp @@ -79,12 +79,14 @@ Double_t PFitterFcn::operator()(const std::vector& par) const value += fRunListCollection->GetSingleHistoChisq(par); value += fRunListCollection->GetSingleHistoRRFChisq(par); value += fRunListCollection->GetAsymmetryChisq(par); + value += fRunListCollection->GetAsymmetryRRFChisq(par); value += fRunListCollection->GetMuMinusChisq(par); value += fRunListCollection->GetNonMusrChisq(par); } else { // max likelihood value += fRunListCollection->GetSingleHistoMaximumLikelihood(par); value += fRunListCollection->GetSingleHistoRRFMaximumLikelihood(par); value += fRunListCollection->GetAsymmetryMaximumLikelihood(par); + value += fRunListCollection->GetAsymmetryRRFMaximumLikelihood(par); value += fRunListCollection->GetMuMinusMaximumLikelihood(par); value += fRunListCollection->GetNonMusrMaximumLikelihood(par); } diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 97cfa207..81b8095b 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -1151,6 +1151,9 @@ Int_t PMsrHandler::WriteMsrLogFile(const Bool_t messages) case MSR_PLOT_ASYM: fout << "PLOT " << fPlots[plotNo].fPlotType << " (asymmetry plot)" << endl; break; + case MSR_PLOT_ASYM_RRF: + fout << "PLOT " << fPlots[plotNo].fPlotType << " (asymmetry RRF plot)" << endl; + break; case MSR_PLOT_MU_MINUS: fout << "PLOT " << fPlots[plotNo].fPlotType << " (mu minus plot)" << endl; break; @@ -2191,6 +2194,9 @@ Int_t PMsrHandler::WriteMsrFile(const Char_t *filename, map *co case MSR_PLOT_ASYM: fout << "PLOT " << fPlots[i].fPlotType << " (asymmetry plot)" << endl; break; + case MSR_PLOT_ASYM_RRF: + fout << "PLOT " << fPlots[i].fPlotType << " (asymmetry RRF plot)" << endl; + break; case MSR_PLOT_MU_MINUS: fout << "PLOT " << fPlots[i].fPlotType << " (mu minus plot)" << endl; break; @@ -5331,7 +5337,51 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity() } break; case PRUN_SINGLE_HISTO_RRF: - // STILL MISSING + // check that there is a forward parameter number + if (fRuns[i].GetForwardHistoNo() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> forward parameter number not defined. Necessary for single histogram RRF fits." << endl; + return false; + } + if ((fRuns[i].GetNormParamNo() > static_cast(fParam.size())) && !fFourierOnly) { + // check if forward histogram number is a function + if (fRuns[i].GetNormParamNo() - MSR_PARAM_FUN_OFFSET > static_cast(fParam.size())) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> forward histogram number " << fRuns[i].GetNormParamNo() << " is larger than the number of fit parameters (" << fParam.size() << ")."; + cerr << endl << ">> Consider to check the manual ;-)" << endl; + return false; + } + } + // check fit range + if (!fRuns[i].IsFitRangeInBin() && !fFourierOnly) { // fit range given as times in usec (RUN block) + if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in RUN block + if (!fGlobal.IsFitRangeInBin()) { // fit range given as times in usec (GLOBAL block) + if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) { // check fit range in GLOBAL block + cerr << endl << "PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << " Fit range is not defined. Necessary for single histogram fits." << endl; + return false; + } + } + } + } + // check number of T0's provided + if ((fRuns[i].GetT0BinSize() > fRuns[i].GetForwardHistoNoSize()) && + (fGlobal.GetT0BinSize() > fRuns[i].GetForwardHistoNoSize())) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl; + cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl; + return false; + } + // check that RRF frequency is given + if (fGlobal.GetRRFUnitTag() == RRF_UNIT_UNDEF) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF frequency found in the GLOBAL block." << endl; + return false; + } + // check that RRF packing is given + if (fGlobal.GetRRFPacking() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF packing found in the GLOBAL block." << endl; + return false; + } break; case PRUN_ASYMMETRY: // check alpha @@ -5386,7 +5436,60 @@ Bool_t PMsrHandler::CheckRunBlockIntegrity() } break; case PRUN_ASYMMETRY_RRF: - // STILL MISSING + // check alpha + if ((fRuns[i].GetAlphaParamNo() == -1) && !fFourierOnly) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> alpha parameter number missing which is needed for an asymmetry RRF fit."; + cerr << endl << ">> Consider to check the manual ;-)" << endl; + return false; + } + // check that there is a forward parameter number + if (fRuns[i].GetForwardHistoNo() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> forward histogram number not defined. Necessary for asymmetry RRF fits." << endl; + return false; + } + // check that there is a backward parameter number + if (fRuns[i].GetBackwardHistoNo() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> backward histogram number not defined. Necessary for asymmetry RRF fits." << endl; + return false; + } + // check fit range + if (!fRuns[i].IsFitRangeInBin()) { // fit range given as times in usec + if ((fRuns[i].GetFitRange(0) == PMUSR_UNDEFINED) || (fRuns[i].GetFitRange(1) == PMUSR_UNDEFINED)) { + if ((fGlobal.GetFitRange(0) == PMUSR_UNDEFINED) || (fGlobal.GetFitRange(1) == PMUSR_UNDEFINED)) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Fit range is not defined, also NOT present in the GLOBAL block. Necessary for asymmetry RRF fits." << endl; + return false; + } + } + } + // check number of T0's provided + if ((fRuns[i].GetT0BinSize() > 2*fRuns[i].GetForwardHistoNoSize()) && + (fGlobal.GetT0BinSize() > 2*fRuns[i].GetForwardHistoNoSize())) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetForwardHistoNoSize() << " in forward. Needs to be fixed." << endl; + cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetForwardHistoNoSize() << ". Needs to be fixed." << endl; + return false; + } + if ((fRuns[i].GetT0BinSize() > 2*fRuns[i].GetBackwardHistoNoSize()) && + (fGlobal.GetT0BinSize() > 2*fRuns[i].GetBackwardHistoNoSize())) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** in RUN block number " << i+1; + cerr << endl << ">> Found " << fRuns[i].GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetBackwardHistoNoSize() << " in backward. Needs to be fixed." << endl; + cerr << endl << ">> In GLOBAL block: " << fGlobal.GetT0BinSize() << " T0 entries. Expecting only " << 2*fRuns[i].GetBackwardHistoNoSize() << ". Needs to be fixed." << endl; + return false; + } + // check that RRF frequency is given + if (fGlobal.GetRRFUnitTag() == RRF_UNIT_UNDEF) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF frequency found in the GLOBAL block." << endl; + return false; + } + // check that RRF packing is given + if (fGlobal.GetRRFPacking() == -1) { + cerr << endl << ">> PMsrHandler::CheckRunBlockIntegrity(): **ERROR** no RRF packing found in the GLOBAL block." << endl; + return false; + } break; case PRUN_MU_MINUS: // needs eventually to be implemented diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index e0d45c53..8cde496b 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -711,7 +711,7 @@ PMsrGlobalBlock::PMsrGlobalBlock() { fGlobalPresent = false; fRRFFreq = 0.0; // rotating reference frequency in units given by fRRFUnitTag. Only needed for fittype 1 - fRRFUnitTag = RRF_UNIT_MHz; // RRF unit tag. Default: MHz + fRRFUnitTag = RRF_UNIT_UNDEF; // RRF unit tag. Default: undefined fRRFPhase = 0.0; fRRFPacking = -1; // undefined RRF packing/rebinning fFitType = -1; // undefined fit type diff --git a/src/classes/PMusrCanvas.cpp b/src/classes/PMusrCanvas.cpp index ccfa539d..ec99e93a 100644 --- a/src/classes/PMusrCanvas.cpp +++ b/src/classes/PMusrCanvas.cpp @@ -4729,6 +4729,7 @@ void PMusrCanvas::PlotData(Bool_t unzoom) } break; case MSR_PLOT_SINGLE_HISTO_RRF: + case MSR_PLOT_ASYM_RRF: yAxisTitle = "RRF Asymmetry"; break; case MSR_PLOT_ASYM: diff --git a/src/classes/PRunAsymmetryRRF.cpp b/src/classes/PRunAsymmetryRRF.cpp index eef4df0c..7e7f21df 100644 --- a/src/classes/PRunAsymmetryRRF.cpp +++ b/src/classes/PRunAsymmetryRRF.cpp @@ -209,7 +209,7 @@ Double_t PRunAsymmetryRRF::CalcChiSquare(const std::vector& par) chunk = 10; #pragma omp parallel for default(shared) private(i,time,diff,asymFcnValue,a,b,f) schedule(dynamic,chunk) reduction(+:chisq) #endif - for (i=startTimeBin; i < endTimeBin; ++i) { + for (i=startTimeBin; i fData.GetValue()->size()=" << fData.GetValue()->size() << endl; - cout << "debug> fFitStartTime=" << fFitStartTime << endl; - cout << "debug> fData.GetDataTimeStart()=" << fData.GetDataTimeStart() << ", fData.GetDataTimeStep()=" << fData.GetDataTimeStep() << endl; - cout << "debug> -----" << endl; - // In order not having to loop over all bins and to stay consistent with the chisq method, calculate the start and end bins explicitly Int_t startTimeBin = static_cast(ceil((fFitStartTime - fData.GetDataTimeStart())/fData.GetDataTimeStep())); if (startTimeBin < 0) @@ -773,10 +768,10 @@ Bool_t PRunAsymmetryRRF::SubtractEstimatedBkg() // calculate proper background range for (UInt_t i=0; i<2; i++) { if (beamPeriod != 0.0) { - Double_t timeBkg = (Double_t)(end[i]-start[i])*(fTimeResolution*fRRFPacking); // length of the background intervall in time + Double_t timeBkg = (Double_t)(end[i]-start[i])*fTimeResolution; // length of the background intervall in time UInt_t fullCycles = (UInt_t)(timeBkg/beamPeriod); // how many proton beam cylces can be placed within the proposed background intervall // correct the end of the background intervall such that the background is as close as possible to a multiple of the proton cylce - end[i] = start[i] + (UInt_t) ((fullCycles*beamPeriod)/(fTimeResolution*fRRFPacking)); + end[i] = start[i] + (UInt_t) ((fullCycles*beamPeriod)/fTimeResolution); cout << "PRunAsymmetryRRF::SubtractEstimatedBkg(): Background " << start[i] << ", " << end[i] << endl; if (end[i] == start[i]) end[i] = fRunInfo->GetBkgRange(2*i+1); @@ -872,8 +867,6 @@ Bool_t PRunAsymmetryRRF::PrepareFitData() Int_t lgb = fgb + lgb_offset; Int_t dt0 = (Int_t)fT0s[0]-(Int_t)fT0s[1]; - cout << "debug> fgb=" << fgb << ", lgb=" << lgb << endl; - PDoubleVector asym; PDoubleVector asymErr; Double_t asymVal, asymValErr; @@ -889,7 +882,7 @@ Bool_t PRunAsymmetryRRF::PrepareFitData() eff = fForwardErr[i]; ebb = fBackwardErr[i-dt0]; if ((asymVal != 0.0) && (ff+bb) > 0.0) - asymValErr = sqrt(2)/(ff+bb)*sqrt(bb*eff+ff*ebb); + asymValErr = 2.0/pow((ff+bb),2.0)*sqrt(bb*bb*eff*eff+ff*ff*ebb*ebb); else asymValErr = 1.0; asymErr.push_back(asymValErr); @@ -907,8 +900,6 @@ Bool_t PRunAsymmetryRRF::PrepareFitData() asym[i] *= 2.0*cos(wRRF*time+phaseRRF); } - cout << "debug> before packing: startTime=" << startTime << ", endTime=" << startTime+asym.size()*fTimeResolution << endl; - // 3rd: rrf packing PDoubleVector asymRRF; asymVal = 0.0; @@ -925,21 +916,12 @@ Bool_t PRunAsymmetryRRF::PrepareFitData() asymValErr = 0.0; for (UInt_t i=0; i after packing: startTime=" << startTime + fTimeResolution*((Double_t)(fRRFPacking-1)/2.0) << ", endTime=" << startTime + fTimeResolution*((Double_t)(fRRFPacking-1)/2.0)+asymRRF.size()*fTimeResolution*(Double_t)(fRRFPacking) << endl; - - ofstream fout("_data.dat", ofstream::out); - for (UInt_t i=0; iGetMsrPlotList()->at(0).fViewPacking > 0) { - packing = fMsrInfo->GetMsrPlotList()->at(0).fViewPacking; - } - // feed the parameter vector std::vector par; PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); for (UInt_t i=0; isize(); i++) par.push_back((*paramList)[i].fValue); - // transform raw histo data. This is done the following way (for details see the manual): - // first rebin the data, than calculate the asymmetry - // first get start data, end data, and t0 Int_t start[2] = {fGoodBins[0], fGoodBins[2]}; Int_t end[2] = {fGoodBins[1], fGoodBins[3]}; @@ -1018,45 +995,32 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2] fgb[0] = start[0]; fgb[1] = start[1]; } + start[0] = fgb[0]; + start[1] = fgb[1]; - Int_t val = fgb[0]-packing*(fgb[0]/packing); - do { - if (fgb[1] - fgb[0] < 0) - val += packing; - } while (val + fgb[1] - fgb[0] < 0); - - start[0] = val; - start[1] = val + fgb[1] - fgb[0]; - - // make sure that there are equal number of rebinned bins in forward and backward - UInt_t noOfBins0 = (runData->GetDataBin(histoNo[0])->size()-start[0])/packing; - UInt_t noOfBins1 = (runData->GetDataBin(histoNo[1])->size()-start[1])/packing; + // make sure that there are equal number of bins in forward and backward + UInt_t noOfBins0 = runData->GetDataBin(histoNo[0])->size()-start[0]; + UInt_t noOfBins1 = runData->GetDataBin(histoNo[1])->size()-start[1]; if (noOfBins0 > noOfBins1) noOfBins0 = noOfBins1; - end[0] = start[0] + noOfBins0 * packing; - end[1] = start[1] + noOfBins0 * packing; + end[0] = start[0] + noOfBins0; + end[1] = start[1] + noOfBins0; // check if start, end, and t0 make any sense - // 1st check if start and end are in proper order for (UInt_t i=0; i<2; i++) { - if (end[i] < start[i]) { // need to swap them - Int_t keep = end[i]; - end[i] = start[i]; - start[i] = keep; - } - // 2nd check if start is within proper bounds + // 1st check if start is within proper bounds if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **ERROR** start data bin doesn't make any sense!"; cerr << endl; return false; } - // 3rd check if end is within proper bounds + // 2nd check if end is within proper bounds if ((end[i] < 0) || (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **ERROR** end data bin doesn't make any sense!"; cerr << endl; return false; } - // 4th check if t0 is within proper bounds + // 3rd check if t0 is within proper bounds if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { cerr << endl << ">> PRunAsymmetryRRF::PrepareViewData(): **ERROR** t0 data bin doesn't make any sense!"; cerr << endl; @@ -1064,72 +1028,15 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2] } } - // everything looks fine, hence fill packed forward and backward histo - PRunData forwardPacked; - PRunData backwardPacked; - Double_t value = 0.0; - Double_t error = 0.0; - - // forward - for (Int_t i=start[0]; i 1 - if (((i-start[0]) % packing == 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 /= packing; - forwardPacked.AppendValue(value); - if (value == 0.0) - forwardPacked.AppendErrorValue(1.0); - else - forwardPacked.AppendErrorValue(TMath::Sqrt(error)/packing); - value = 0.0; - error = 0.0; - } - value += fForward[i]; - error += fForwardErr[i]*fForwardErr[i]; - } - } - - // backward - for (Int_t i=start[1]; i 1 - if (((i-start[1]) % packing == 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 /= packing; - backwardPacked.AppendValue(value); - if (value == 0.0) - backwardPacked.AppendErrorValue(1.0); - else - backwardPacked.AppendErrorValue(TMath::Sqrt(error)/packing); - value = 0.0; - error = 0.0; - } - value += fBackward[i]; - error += fBackwardErr[i]*fBackwardErr[i]; - } - } - - // check if packed forward and backward hist have the same size, otherwise take the minimum size - UInt_t noOfBins = forwardPacked.GetValue()->size(); - if (forwardPacked.GetValue()->size() != backwardPacked.GetValue()->size()) { - if (forwardPacked.GetValue()->size() > backwardPacked.GetValue()->size()) - noOfBins = backwardPacked.GetValue()->size(); + // check if forward and backward histo have the same size, otherwise take the minimum size + UInt_t noOfBins = fForward.size(); + if (noOfBins > fBackward.size()) { + noOfBins = fBackward.size(); } // form asymmetry including error propagation - Double_t asym; + Double_t asym, error; Double_t f, b, ef, eb, alpha = 1.0, beta = 1.0; - // set data time start, and step - // data start at data_start-t0 - fData.SetDataTimeStart(fTimeResolution*((Double_t)start[0]-t0[0]+(Double_t)(packing-1)/2.0)); - fData.SetDataTimeStep(fTimeResolution*(Double_t)packing); // get the proper alpha and beta switch (fAlphaBetaTag) { @@ -1153,24 +1060,67 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2] break; } - for (UInt_t i=0; isize(); i++) { + PDoubleVector asymVec, asymErr; + Int_t dtBin = start[1]-start[0]; + for (Int_t i=start[0]; iat(i); - b = backwardPacked.GetValue()->at(i); - ef = forwardPacked.GetError()->at(i); - eb = backwardPacked.GetError()->at(i); + f = fForward[i]; + b = fBackward[i+dtBin]; + ef = fForwardErr[i]; + eb = fBackwardErr[i+dtBin]; // check that there are indeed bins if (f+b != 0.0) asym = (alpha*f-b) / (alpha*beta*f+b); else asym = 0.0; - fData.AppendValue(asym); + asymVec.push_back(asym); // calculate the error if (f+b != 0.0) error = 2.0/((f+b)*(f+b))*TMath::Sqrt(b*b*ef*ef+eb*eb*f*f); else error = 1.0; - fData.AppendErrorValue(error); + asymErr.push_back(error); + } + + // RRF transform + // a_rrf = a * 2*cos(w_rrf*t + phi_rrf) + PMsrGlobalBlock *globalBlock = fMsrInfo->GetMsrGlobal(); + Double_t wRRF = globalBlock->GetRRFFreq("Mc"); + Double_t phaseRRF = globalBlock->GetRRFPhase()*TMath::TwoPi()/180.0; + Double_t startTime=fTimeResolution*((Double_t)start[0]-t0[0]); + Double_t time = 0.0; + for (UInt_t i=0; iGetDataBin(histoNo[0])->size(); Double_t factor = 1.0; if (fData.GetValue()->size() * 10 > runData->GetDataBin(histoNo[0])->size()) { @@ -1199,11 +1152,11 @@ Bool_t PRunAsymmetryRRF::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2] fData.SetTheoryTimeStep(fTimeResolution*factor); for (UInt_t i=0; iFunc(time, par, fFuncValues); - if (fabs(value) > 10.0) { // dirty hack needs to be fixed!! - value = 0.0; + dval = fTheory->Func(time, par, fFuncValues); + if (fabs(dval) > 10.0) { // dirty hack needs to be fixed!! + dval = 0.0; } - fData.AppendTheoryValue(value); + fData.AppendTheoryValue(dval); } // clean up