diff --git a/src/classes/PMusr.cpp b/src/classes/PMusr.cpp index a84f4765..6a2142df 100644 --- a/src/classes/PMusr.cpp +++ b/src/classes/PMusr.cpp @@ -1879,6 +1879,19 @@ void PMsrRunBlock::SetMapGlobal(UInt_t idx, Int_t ival) return; } +//-------------------------------------------------------------------------- +// SetEstimatedAlpha (public) +//-------------------------------------------------------------------------- +/** + *

set the value of estimated alpha at position idx + * + * \param alpha is the estimated value + */ +void PMsrRunBlock::SetEstimatedAlpha(Double_t dval) +{ + fAlpha = dval; +} + //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ // implementation PStringNumberList //+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ diff --git a/src/classes/PRunAsymmetryBNMR.cpp b/src/classes/PRunAsymmetryBNMR.cpp index 83d27d2a..8afa4937 100644 --- a/src/classes/PRunAsymmetryBNMR.cpp +++ b/src/classes/PRunAsymmetryBNMR.cpp @@ -106,12 +106,13 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD // check if alpha is given Bool_t alphaFixedToOne = false; + Bool_t alphaSetDefault = false; if (fRunInfo->GetAlphaParamNo() == -1) { // no alpha given // cerr << endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** no alpha parameter given! This is needed for an asymmetry fit!"; // cerr << endl; // fValid = false; // return; - alphaFixedToOne = true; + alphaSetDefault = true; } else if ((fRunInfo->GetAlphaParamNo() < 0) || (fRunInfo->GetAlphaParamNo() > (Int_t)param->size())) { // check if alpha parameter is within proper bounds cerr << endl << ">> PRunAsymmetryBNMR::PRunAsymmetryBNMR(): **ERROR** alpha parameter no = " << fRunInfo->GetAlphaParamNo(); cerr << endl << ">> This is out of bound, since there are only " << param->size() << " parameters."; @@ -141,12 +142,16 @@ PRunAsymmetryBNMR::PRunAsymmetryBNMR(PMsrHandler *msrInfo, PRunDataHandler *rawD } // set fAlphaBetaTag - if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1 + if (alphaFixedToOne && betaFixedToOne) // alpha == 1, beta == 1 fAlphaBetaTag = 1; - else if (!alphaFixedToOne && betaFixedToOne) // alpha != 1, beta == 1 + else if (!alphaFixedToOne && betaFixedToOne & !alphaSetDefault) // alpha != 1, beta == 1 fAlphaBetaTag = 2; - else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1 + else if (alphaFixedToOne && !betaFixedToOne) // alpha == 1, beta != 1 fAlphaBetaTag = 3; + else if (!alphaFixedToOne && betaFixedToOne & alphaSetDefault) // alpha ??, beta == 1 + fAlphaBetaTag = 5; + else if (!alphaFixedToOne && !betaFixedToOne & alphaSetDefault) // alpha ??, beta != 1 + fAlphaBetaTag = 6; else fAlphaBetaTag = 4; @@ -198,7 +203,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) } // calculate chi square - Double_t time(1.0); + Double_t time(1.0),alphaest; Int_t i; // Calculate the theory function once to ensure one function evaluation for the current set of parameters. @@ -206,6 +211,7 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) // for a given set of parameters---which should be done outside of the parallelized loop. // For all other functions it means a tiny and acceptable overhead. asymFcnValue = fTheory->Func(time, par, fFuncValues); + alphaest = fRunInfo->GetEstimatedAlpha(); #ifdef HAVE_GOMP Int_t chunk = (fEndTimeBin - fStartTimeBin)/omp_get_num_procs(); @@ -235,6 +241,17 @@ Double_t PRunAsymmetryBNMR::CalcChiSquare(const std::vector& par) f = fTheory->Func(time, par, fFuncValues)/2.0; asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; + case 5: // alpha ?? , beta == 1 + a = alphaest; + f = fTheory->Func(time, par, fFuncValues)/2.0; + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); + break; + case 6: // alpha ??, beta != 1 + a = alphaest; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues)/2.0; + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); + break; default: asymFcnValue = 0.0; break; @@ -425,8 +442,12 @@ void PRunAsymmetryBNMR::CalcTheory() // calculate asymmetry Double_t asymFcnValue = 0.0; - Double_t a, b, f; + Double_t a, b, f, alphaest; Double_t time; + + // Get estimated alpha + alphaest = fRunInfo->GetEstimatedAlpha(); + for (UInt_t i=0; isize(); i++) { time = fData.GetDataTimeStart() + (Double_t)i*fData.GetDataTimeStep(); switch (fAlphaBetaTag) { @@ -436,18 +457,29 @@ void PRunAsymmetryBNMR::CalcTheory() case 2: // alpha != 1, beta == 1 a = par[fRunInfo->GetAlphaParamNo()-1]; f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); break; case 3: // alpha == 1, beta != 1 b = par[fRunInfo->GetBetaParamNo()-1]; f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0)); + asymFcnValue = f*(b+1.0)/(2.0-f*(b-1.0))-f*(b+1.0)/(2.0+f*(b-1.0)); break; case 4: // alpha != 1, beta != 1 a = par[fRunInfo->GetAlphaParamNo()-1]; b = par[fRunInfo->GetBetaParamNo()-1]; f = fTheory->Func(time, par, fFuncValues); - asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0)); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); + break; + case 5: // alpha ?? , beta == 1 + a = alphaest; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a+1.0)-(a-1.0))/((a+1.0)-f*(a-1.0)) - (-f*(a+1.0)-(a-1.0))/((a+1.0)+f*(a-1.0)); + break; + case 6: // alpha ??, beta != 1 + a = alphaest; + b = par[fRunInfo->GetBetaParamNo()-1]; + f = fTheory->Func(time, par, fFuncValues); + asymFcnValue = (f*(a*b+1.0)-(a-1.0))/((a+1.0)-f*(a*b-1.0))-(-f*(a*b+1.0)-(a-1.0))/((a+1.0)+f*(a*b-1.0)); break; default: asymFcnValue = 0.0; @@ -871,6 +903,14 @@ Bool_t PRunAsymmetryBNMR::SubtractEstimatedBkg() fRunInfo->SetBkgEstimated(bkgm[0], 3); fRunInfo->SetBkgEstimated(bkgm[1], 4); + cout << "I am here 1 " << endl; + // Get estimate for alpha once + Double_t alpha; + alpha = EstimateAlpha(); + cout << "I am here 2 " << alpha << endl; + fRunInfo->SetEstimatedAlpha(alpha); + + return true; } @@ -903,14 +943,8 @@ Bool_t PRunAsymmetryBNMR::PrepareFitData() Double_t valuem = 0.0; Double_t errorm = 0.0; - Double_t SumF[2] = {0.0, 0.0}; - Double_t SumB[2] = {0.0, 0.0}; - Double_t alphaest = 1; - // forward for (Int_t i=fGoodBins[0]; i> PRunAsymmetryBNMR::PrepareFitData(): alpha estimate=" << alphaest << endl; - // check if packed forward and backward hist have the same size, otherwise take the minimum size UInt_t noOfBins = forwardpPacked.GetValue()->size(); if (forwardpPacked.GetValue()->size() != backwardpPacked.GetValue()->size()) { @@ -1249,13 +1277,16 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 // form asymmetry including error propagation Double_t asym; - Double_t fp, bp, efp, ebp, alpha = 1.0, beta = 1.0; + Double_t fp, bp, efp, ebp, alpha, beta = 1.0; Double_t fm, bm, efm, ebm; // 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 estimate for alpha once + alpha = EstimateAlpha(); + // get the proper alpha and beta switch (fAlphaBetaTag) { case 1: // alpha == 1, beta == 1 @@ -1274,6 +1305,14 @@ Bool_t PRunAsymmetryBNMR::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2 alpha = par[fRunInfo->GetAlphaParamNo()-1]; beta = par[fRunInfo->GetBetaParamNo()-1]; break; + case 5: // alpha ?? , beta == 1 + // use estimated value + beta = 1.0; + break; + case 6: // alpha ??, beta != 1 + // use estimated value + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; default: break; } @@ -1701,3 +1740,43 @@ void PRunAsymmetryBNMR::GetProperFitRange(PMsrGlobalBlock *globalBlock) cerr << ">> Will set it to fgb/lgb which given in time is: " << fFitStartTime << "..." << fFitEndTime << " (usec)" << endl; } } + + +//-------------------------------------------------------------------------- +// EstimateAlpha (private) +//-------------------------------------------------------------------------- +/** + *

Get an estimate for alpha from the forward and backward histograms + * + * \param globalBlock pointer to the GLOBAL block information form the msr-file. + */ +Double_t PRunAsymmetryBNMR::EstimateAlpha() +{ + + Double_t SumF[2] = {0.0, 0.0}; + Double_t SumB[2] = {0.0, 0.0}; + Double_t alpha = 1; + + + // forward + for (Int_t i=0; i> PRunAsymmetryBNMR::EstimateAlpha(): alpha estimate=" << alpha << endl; + + return alpha; +} diff --git a/src/include/PMusr.h b/src/include/PMusr.h index cd9d6375..357cf6b4 100644 --- a/src/include/PMusr.h +++ b/src/include/PMusr.h @@ -570,6 +570,7 @@ class PMsrGlobalBlock { virtual Double_t GetFitRange(UInt_t idx); virtual Int_t GetFitRangeOffset(UInt_t idx); virtual Int_t GetPacking() { return fPacking; } + virtual Double_t GetEstimatedAlpha() { return fAlpha; } virtual void SetGlobalPresent(Bool_t bval) { fGlobalPresent = bval; } virtual void SetRRFFreq(Double_t freq, const char *unit); @@ -598,6 +599,7 @@ class PMsrGlobalBlock { Double_t fFitRange[2]; ///< fit range in (us) Int_t fFitRangeOffset[2]; ///< if fit range is given in bins it can have the form fit fgb+n0 lgb-n1. This variable holds the n0 and n1. Int_t fPacking; ///< packing/rebinning + Double_t fAlpha; ///< estimated alpha value from F/B counts }; //------------------------------------------------------------- @@ -645,6 +647,7 @@ class PMsrRunBlock { virtual Double_t GetFitRange(UInt_t idx); virtual Int_t GetFitRangeOffset(UInt_t idx); virtual Int_t GetPacking() { return fPacking; } + virtual Double_t GetEstimatedAlpha() { return fAlpha; } virtual Int_t GetXDataIndex() { return fXYDataIndex[0]; } virtual Int_t GetYDataIndex() { return fXYDataIndex[1]; } virtual TString* GetXDataLabel() { return &fXYDataLabel[0]; } @@ -667,6 +670,7 @@ class PMsrRunBlock { virtual void SetForwardHistoNo(Int_t histoNo, Int_t idx=-1); virtual void SetBackwardHistoNo(Int_t histoNo, Int_t idx=-1); virtual void SetBkgEstimated(Double_t dval, Int_t idx); + virtual void SetEstimatedAlpha(Double_t dval); virtual void SetBkgFix(Double_t dval, Int_t idx); virtual void SetBkgRange(Int_t ival, Int_t idx); virtual void SetDataRange(Int_t ival, Int_t idx); @@ -707,6 +711,7 @@ class PMsrRunBlock { Bool_t fFitRangeInBins; ///< flag telling if fit range is given in time or in bins Double_t fFitRange[2]; ///< fit range in (us) Int_t fFitRangeOffset[2]; ///< if fit range is given in bins it can have the form fit fgb+n0 lgb-n1. This variable holds the n0 and n1. + Double_t fAlpha; ///< estimated alpha value from F/B counts Int_t fPacking; ///< packing/rebinning Int_t fXYDataIndex[2]; ///< used to get the data indices when using db-files (fit type 8) TString fXYDataLabel[2]; ///< used to get the indices via labels when using db-files (fit type 8) diff --git a/src/include/PRunAsymmetryBNMR.h b/src/include/PRunAsymmetryBNMR.h index a1fc5377..5ed85c9f 100644 --- a/src/include/PRunAsymmetryBNMR.h +++ b/src/include/PRunAsymmetryBNMR.h @@ -87,6 +87,7 @@ class PRunAsymmetryBNMR : public PRunBase virtual Bool_t GetProperT0(PRawRunData* runData, PMsrGlobalBlock *globalBlock, PUIntVector &forwardHisto, PUIntVector &backwardHistoNo); virtual Bool_t GetProperDataRange(PRawRunData* runData, UInt_t histoNo[2]); virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock); + virtual Double_t EstimateAlpha(); }; #endif // _PRUNASYMMETRYBNMR_H_