44 #include <TObjArray.h> 45 #include <TObjString.h> 86 PRunBase(msrInfo, rawData, runNo, tag), fTheoAsData(theoAsData)
93 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no GLOBAL-block present!";
94 std::cerr << std::endl <<
">> For Single Histo RRF the GLOBAL-block is mandatory! Please fix this first.";
95 std::cerr << std::endl;
101 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Frequency found!";
102 std::cerr << std::endl;
109 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: no RRF-Packing found!";
110 std::cerr << std::endl;
123 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PRunSingleHistoRRF(): **SEVERE ERROR**: Couldn't prepare data for fitting!";
124 std::cerr << std::endl <<
">> This is very bad :-(, will quit ...";
125 std::cerr << std::endl;
154 Double_t chisq = 0.0;
177 #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) 201 Double_t chisq = 0.0;
225 #pragma omp parallel for default(shared) private(i,time,diff) schedule(dynamic,chunk) reduction(+:chisq) 231 chisq += diff*diff / theo;
264 std::vector<Double_t> par;
266 for (UInt_t i=0; i<paramList->size(); i++)
267 par.push_back((*paramList)[i].fValue);
279 for (UInt_t i=0; i<size; i++) {
280 time = start +
static_cast<Double_t
>(i)*resolution;
319 TObjArray *tok =
nullptr;
320 TObjString *ostr =
nullptr;
325 tok = fitRange.Tokenize(
" \t");
327 if (tok->GetEntries() == 3) {
329 ostr =
dynamic_cast<TObjString*
>(tok->At(1));
330 str = ostr->GetString();
332 idx = str.First(
"+");
334 str.Remove(0, idx+1);
341 ostr =
dynamic_cast<TObjString*
>(tok->At(2));
342 str = ostr->GetString();
344 idx = str.First(
"-");
346 str.Remove(0, idx+1);
351 }
else if ((tok->GetEntries() > 3) && (tok->GetEntries() % 2 == 1)) {
352 Int_t pos = 2*(
fRunNo+1)-1;
354 if (pos + 1 >= tok->GetEntries()) {
355 std::cerr << std::endl <<
">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange <<
"'";
356 std::cerr << std::endl <<
">> will ignore it. Sorry ..." << std::endl;
359 ostr =
dynamic_cast<TObjString*
>(tok->At(pos));
360 str = ostr->GetString();
362 idx = str.First(
"+");
364 str.Remove(0, idx+1);
371 ostr =
dynamic_cast<TObjString*
>(tok->At(pos+1));
372 str = ostr->GetString();
374 idx = str.First(
"-");
376 str.Remove(0, idx+1);
383 std::cerr << std::endl <<
">> PRunSingleHistoRRF::SetFitRangeBin(): **ERROR** invalid FIT_RANGE command found: '" << fitRange <<
"'";
384 std::cerr << std::endl <<
">> will ignore it. Sorry ..." << std::endl;
433 Bool_t success =
true;
444 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get run " <<
fRunInfo->
GetRunName()->Data() <<
"!";
445 std::cerr << std::endl;
465 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PrepareData(): **PANIC ERROR**:";
466 std::cerr << std::endl <<
">> histoNo found = " << histoNo[i] <<
", which is NOT present in the data file!?!?";
467 std::cerr << std::endl <<
">> Will quit :-(";
468 std::cerr << std::endl;
476 std::cout.precision(10);
477 std::cout << std::endl <<
">> PRunSingleHisto::PrepareData(): time resolution=" << std::fixed << runData->
GetTimeResolution() <<
"(ns)" << std::endl;
485 std::vector<PDoubleVector> forward;
486 forward.resize(histoNo.size());
487 for (UInt_t i=0; i<histoNo.size(); i++) {
488 forward[i].resize(runData->
GetDataBin(histoNo[i])->size());
489 forward[i] = *runData->
GetDataBin(histoNo[i]);
499 if (addRunData ==
nullptr) {
500 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PrepareData(): **ERROR** Couldn't get addrun " <<
fRunInfo->
GetRunName(i)->Data() <<
"!";
501 std::cerr << std::endl;
507 for (UInt_t k=0; k<histoNo.size(); k++) {
508 addRunSize = addRunData->
GetDataBin(histoNo[k])->size();
509 for (UInt_t j=0; j<addRunData->
GetDataBin(histoNo[k])->size(); j++) {
511 if ((static_cast<Int_t>(j)+static_cast<Int_t>(
fAddT0s[i-1][k])-static_cast<Int_t>(
fT0s[k]) >= 0) &&
512 (j+
static_cast<Int_t
>(
fAddT0s[i-1][k])-static_cast<Int_t>(
fT0s[k]) < addRunSize)) {
513 forward[k][j] += addRunData->
GetDataBin(histoNo[k])->at(j+static_cast<Int_t>(
fAddT0s[i-1][k])-static_cast<Int_t>(
fT0s[k]));
522 for (UInt_t i=0; i<
fForward.size(); i++) {
527 for (UInt_t i=1; i<histoNo.size(); i++) {
528 for (UInt_t j=0; j<runData->
GetDataBin(histoNo[i])->size(); j++) {
530 if ((static_cast<Int_t>(j)+static_cast<Int_t>(
fT0s[i])-static_cast<Int_t>(
fT0s[0]) >= 0) &&
531 (j+
static_cast<Int_t
>(
fT0s[i])-static_cast<Int_t>(
fT0s[0]) < runData->
GetDataBin(histoNo[i])->size())) {
532 fForward[j] += forward[i][j+
static_cast<Int_t
>(
fT0s[i])-static_cast<Int_t>(
fT0s[0])];
583 for (UInt_t i=0; i<
fForward.size(); i++) {
587 std::cout <<
"info> freqMax=" << freqMax <<
" (MHz)" << std::endl;
590 Double_t optNoPoints = 8;
606 std::cerr << std::endl <<
">> PRunSingleHistoRRF::PrepareFitData(): **WARNING** Neither fix background nor background bins are given!";
608 std::cerr << std::endl <<
">> NO WARRANTY THAT THIS MAKES ANY SENSE! Better check ...";
609 std::cerr << std::endl;
614 for (UInt_t i=0; i<
fForward.size(); i++)
617 for (UInt_t i=0; i<
fForward.size(); i++) {
624 Int_t t0 =
static_cast<Int_t
>(
fT0s[0]);
629 Double_t time_tau=0.0;
630 Double_t exp_t_tau=0.0;
633 exp_t_tau = exp(time_tau);
640 for (UInt_t i=0; i<
fMerr.size(); i++) {
649 Double_t errN0 = 0.0;
660 exp_t_tau = exp(time_tau);
661 fAerr.push_back(exp_t_tau/n0*sqrt(rawNt[i]+pow(((rawNt[i]-
fBackground)/n0)*errN0,2.0)));
666 Double_t wRRF = globalBlock->
GetRRFFreq(
"Mc");
667 Double_t phaseRRF = globalBlock->
GetRRFPhase()*TMath::TwoPi()/180.0;
671 fForward[i] *= 2.0*cos(wRRF * time + phaseRRF);
741 if (viewPacking > 0) {
743 std::cerr <<
">> PRunSingleHistoRRF::PrepareViewData(): **WARNING** Found View Packing (" << viewPacking <<
") < RRF Packing (" <<
fRRFPacking <<
").";
744 std::cerr <<
">> Will ignore View Packing." << std::endl;
755 std::vector<Double_t> par;
757 for (UInt_t i=0; i<paramList->size(); i++)
758 par.push_back((*paramList)[i].fValue);
779 Double_t theoryValue = 0.0;
780 for (UInt_t i=0; i<size; i++) {
783 if (fabs(theoryValue) > 10.0) {
817 fT0s.resize(histoNo.size());
818 for (UInt_t i=0; i<
fT0s.size(); i++) {
829 if (
fT0s[i] == -1.0) {
835 for (UInt_t i=0; i<histoNo.size(); i++) {
836 if (
fT0s[i] == -1.0) {
837 if (runData->
GetT0Bin(histoNo[i]) > 0.0) {
845 for (UInt_t i=0; i<histoNo.size(); i++) {
846 if (
fT0s[i] == -1.0) {
850 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
852 std::cerr << std::endl <<
">> will try the estimated one: forward t0 = " << runData->
GetT0BinEstimated(histoNo[i]);
853 std::cerr << std::endl <<
">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
854 std::cerr << std::endl;
860 if ((
fT0s[i] < 0.0) || (
fT0s[i] >
static_cast<Int_t
>(runData->
GetDataBin(histoNo[i])->size()))) {
861 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperT0(): **ERROR** t0 data bin (" <<
fT0s[i] <<
") doesn't make any sense!";
862 std::cerr << std::endl;
875 if (addRunData ==
nullptr) {
876 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperT0(): **ERROR** Couldn't get addrun " <<
fRunInfo->
GetRunName(i)->Data() <<
"!";
877 std::cerr << std::endl;
883 fAddT0s[i-1].resize(histoNo.size());
884 for (UInt_t j=0; j<
fAddT0s[i-1].size(); j++) {
894 for (UInt_t j=0; j<histoNo.size(); j++) {
896 if (addRunData->
GetT0Bin(histoNo[j]) > 0.0) {
903 for (UInt_t j=0; j<histoNo.size(); j++) {
908 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperT0(): **WARRNING** NO t0's found, neither in the run data nor in the msr-file!";
910 std::cerr << std::endl <<
">> will try the estimated one: forward t0 = " << addRunData->
GetT0BinEstimated(histoNo[j]);
911 std::cerr << std::endl <<
">> NO WARRANTY THAT THIS OK!! For instance for LEM this is almost for sure rubbish!";
912 std::cerr << std::endl;
919 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperT0(): **ERROR** addt0 data bin (" <<
fAddT0s[i-1][j] <<
") doesn't make any sense!";
920 std::cerr << std::endl;
962 start =
static_cast<Int_t
>(
fT0s[0])+offset;
964 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range start = t0+" << offset <<
"(=10ns) = " << start <<
".";
965 std::cerr << std::endl <<
">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
966 std::cerr << std::endl;
971 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** data range was not provided, will try data range end = " << end <<
".";
972 std::cerr << std::endl <<
">> NO WARRANTY THAT THIS DOES MAKE ANY SENSE.";
973 std::cerr << std::endl;
984 if ((start < 0) || (start > static_cast<Int_t>(
fForward.size()))) {
985 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** start data bin (" << start <<
") doesn't make any sense!";
986 std::cerr << std::endl;
991 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperDataRange(): **ERROR** end data bin (" << end <<
") doesn't make any sense!";
992 std::cerr << std::endl;
995 if (end > static_cast<Int_t>(
fForward.size())) {
996 std::cerr << std::endl <<
">> PRunSingleHistoRRF::GetProperDataRange(): **WARNING** end data bin (" << end <<
") > histo length (" <<
fForward.size() <<
").";
997 std::cerr << std::endl <<
">> Will set end = (histo length - 1). Consider to change it in the msr-file." << std::endl;
998 std::cerr << std::endl;
999 end =
static_cast<Int_t
>(
fForward.size()-1);
1052 std::cerr <<
">> PRunSingleHistoRRF::GetProperFitRange(): **WARNING** Couldn't get fit start/end time!" << std::endl;
1053 std::cerr <<
">> Will set it to fgb/lgb which given in time is: " <<
fFitStartTime <<
"..." <<
fFitEndTime <<
" (usec)" << std::endl;
1067 Double_t freqMax = 0.0;
1074 histo->SetBinContent(i-
fGoodBins[0]+1, data[i]);
1078 std::unique_ptr<PFourier> ft = std::make_unique<PFourier>(histo.get(),
FOURIER_UNIT_FREQ);
1082 TH1F *power = ft->GetPowerFourier();
1083 Double_t maxFreqVal = 0.0;
1084 for (Int_t i=1; i<power->GetNbinsX(); i++) {
1086 if (i<power->GetNbinsX()-1) {
1087 if (power->GetBinContent(i)>power->GetBinContent(i+1))
1091 if (power->GetBinCenter(i) < 10.0)
1094 if (power->GetBinContent(i) > maxFreqVal) {
1095 maxFreqVal = power->GetBinContent(i);
1096 freqMax = power->GetBinCenter(i);
1123 for (Int_t i=0; i<endBin; i++) {
1132 for (Int_t i=0; i<endBin; i++) {
1135 errN0 = sqrt(errN0)/wN;
1137 std::cout <<
"info> PRunSingleHistoRRF::EstimateN0(): N0=" << n0 <<
"(" << errN0 <<
")" << std::endl;
1156 Double_t beamPeriod = 0.0;
1172 std::cout << std::endl <<
"PRunSingleHistoRRF::EstimatBkg(): end = " << end <<
" > start = " << start <<
"! Will swap them!";
1179 if (beamPeriod != 0.0) {
1181 UInt_t fullCycles =
static_cast<UInt_t
>(timeBkg/beamPeriod);
1183 end = start +
static_cast<UInt_t
>((fullCycles*beamPeriod)/
fTimeResolution);
1184 std::cout << std::endl <<
"PRunSingleHistoRRF::EstimatBkg(): Background " << start <<
", " << end;
1191 std::cerr << std::endl <<
">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!";
1192 std::cerr << std::endl <<
">> histo lengths = " <<
fForward.size();
1193 std::cerr << std::endl <<
">> background start = " << start;
1194 std::cerr << std::endl;
1200 std::cerr << std::endl <<
">> PRunSingleHistoRRF::EstimatBkg(): **ERROR** background bin values out of bound!";
1201 std::cerr << std::endl <<
">> histo lengths = " <<
fForward.size();
1202 std::cerr << std::endl <<
">> background end = " << end;
1203 std::cerr << std::endl;
1211 for (UInt_t i=start; i<end; i++)
1213 bkg /=
static_cast<Double_t
>(end - start + 1);
1218 for (UInt_t i=start; i<end; i++)
1220 fBkgErr = sqrt(bkg/(static_cast<Double_t>(end - start)));
1222 std::cout << std::endl <<
"info> fBackground=" <<
fBackground <<
"(" <<
fBkgErr <<
")" << std::endl;
Int_t fGoodBins[2]
keep first/last good bins. 0=fgb, 1=lgb
virtual UInt_t GetT0BinSize()
Bool_t fValid
flag showing if the state of the class is valid
virtual TString GetRRFUnit()
virtual void GetProperFitRange(PMsrGlobalBlock *globalBlock)
virtual void SetDataRange(Int_t ival, Int_t idx)
virtual Double_t CalcMaxLikelihood(const std::vector< Double_t > &par)
virtual Double_t EstimateN0(Double_t &errN0, Double_t freqMax)
virtual const PDoubleVector * GetError()
PDoubleVector fMerr
vector holding the error of M(t): M_err = exp(+t/tau) sqrt(N(t)).
Double_t fN0EstimateEndTime
end time in (us) over which N0 is estimated.
virtual Double_t CalcChiSquare(const std::vector< Double_t > &par)
virtual Int_t GetNoOfFuncs()
virtual Int_t GetDataRange(UInt_t idx)
virtual Double_t GetT0Bin(UInt_t idx=0)
Double_t fBkgErr
estimate error on the estimated background
virtual const Double_t GetTimeResolution()
virtual UInt_t GetNoOfFitBins()
PRunDataHandler * fRawData
holds the raw run data
std::vector< PDoubleVector > fAddT0s
all t0 bins of all addrun's of a run! The derived classes will handle it.
virtual void SetBkgEstimated(Double_t dval, Int_t idx)
virtual Int_t GetFitRangeOffset(UInt_t idx)
virtual Bool_t GetProperDataRange()
virtual Bool_t GetProperT0(PRawRunData *runData, PMsrGlobalBlock *globalBlock, PUIntVector &histoNo)
virtual void SetTheoryTimeStart(Double_t dval)
PMsrHandler * fMsrInfo
msr-file handler
virtual Bool_t IsFitRangeInBin()
virtual const PDoublePairVector * GetTemperature() const
Double_t fTimeResolution
time resolution in (us)
PMetaData fMetaData
keeps the meta data from the data file like field, temperature, energy, ...
#define F_APODIZATION_STRONG
Int_t fEndTimeBin
bin at which the fit ends
virtual PRawRunData * GetRunData(const TString &runName)
PDoubleVector fM
vector holding M(t) = [N(t)-N_bkg] exp(+t/tau). Needed to estimate N0.
virtual Double_t GetAddT0Bin(UInt_t addRunIdx, UInt_t histoIdx)
virtual Int_t GetRRFPacking()
virtual UInt_t GetT0BinSize()
PDoubleVector fAerr
vector holding the errors of estimated A(t)
virtual void CalcTheory()
PDoubleVector fForward
forward histo data
virtual Double_t GetMainFrequency(PDoubleVector &data)
virtual Bool_t EstimateBkg(UInt_t histoNo)
virtual void AppendErrorValue(Double_t dval)
UInt_t fNoOfFitBins
number of bins to be fitted
virtual Bool_t IsFitRangeInBin()
virtual Double_t GetTheoryTimeStep()
std::vector< UInt_t > PUIntVector
PRunData fData
data to be fitted, viewed, i.e. binned data
virtual void SetDataTimeStep(Double_t dval)
virtual Double_t GetTheoryTimeStart()
virtual Bool_t IsPresent()
virtual Double_t GetT0Bin(UInt_t idx=0)
std::vector< PMsrParamStructure > PMsrParamList
std::vector< Double_t > PDoubleVector
virtual Double_t GetRRFFreq(const char *unit)
virtual const Bool_t IsPresent(UInt_t histoNo)
#define ACCEL_PERIOD_TRIUMF
virtual Bool_t PrepareData()
#define FOURIER_UNIT_FREQ
virtual const Double_t GetT0Bin(const UInt_t histoNo)
virtual Bool_t PrepareViewData(PRawRunData *runData, const UInt_t histoNo)
virtual void SetT0Bin(Double_t dval, Int_t idx=-1)
virtual void AppendTheoryValue(Double_t dval)
virtual PMsrGlobalBlock * GetMsrGlobal()
virtual void SetDataTimeStart(Double_t dval)
virtual Int_t GetFitRangeOffset(UInt_t idx)
virtual void AppendValue(Double_t dval)
virtual const PDoubleVector * GetValue()
virtual const Double_t GetField()
virtual void SetAddT0Bin(Double_t dval, UInt_t addRunIdx, UInt_t histoNoIdx)
Int_t fStartTimeBin
bin at which the fit starts
virtual Bool_t PrepareFitData(PRawRunData *runData, const UInt_t histoNo)
virtual Double_t GetBkgFix(UInt_t idx)
virtual void SetFitRangeBin(const TString fitRange)
virtual UInt_t GetFuncNo(Int_t idx)
PDoubleVector fW
vector holding the weight needed to estimate N0, and errN0.
virtual const UInt_t GetNoOfTemperatures()
virtual Int_t GetForwardHistoNo(UInt_t idx=0)
PDoubleVector fT0s
all t0 bins of a run! The derived classes will handle it.
virtual PMsrPlotList * GetMsrPlotList()
virtual Int_t GetBkgRange(UInt_t idx)
EPMusrHandleTag fHandleTag
tag telling whether this is used for fit, view, ...
virtual UInt_t GetForwardHistoNoSize()
virtual Double_t GetRRFPhase()
PMsrRunBlock * fRunInfo
run info used to filter out needed infos of a run
virtual const Double_t GetEnergy()
virtual void SetBkgRange(Int_t ival, Int_t idx)
virtual void SetFitRange(Double_t dval, UInt_t idx)
Double_t fFitEndTime
fit end time
virtual void SetTheoryTimeStep(Double_t dval)
virtual TString * GetInstitute(UInt_t idx=0)
virtual void SetFitRange(Double_t dval, UInt_t idx)
Double_t fFitStartTime
fit start time
virtual PMsrParamList * GetMsrParamList()
virtual PIntVector * GetMap()
virtual Double_t EvalFunc(UInt_t i, std::vector< Int_t > map, std::vector< Double_t > param, PMetaData metaData)
virtual const PDoubleVector * GetDataBin(const UInt_t histoNo)
Bool_t fTheoAsData
true=only calculate the theory points at the data points, false=calculate more points for the theory ...
virtual Double_t GetDataTimeStart()
virtual UInt_t GetRunNameSize()
virtual ~PRunSingleHistoRRF()
Double_t fBackground
needed if background range is given (units: 1/bin)
virtual Int_t GetDataRange(UInt_t idx)
virtual const Double_t GetT0BinEstimated(const UInt_t histoNo)
Int_t fRRFPacking
RRF packing for this particular run. Given in the GLOBAL-block.
virtual void CalcNoOfFitBins()
std::unique_ptr< PTheory > fTheory
theory needed to calculate chi-square
virtual Double_t GetFitRange(UInt_t idx)
virtual Double_t GetDataTimeStep()
virtual Double_t CalcChiSquareExpected(const std::vector< Double_t > &par)
Int_t fRunNo
number of the run within the msr-file
virtual TString * GetRunName(UInt_t idx=0)
virtual Double_t GetFitRange(UInt_t idx)
PDoubleVector fFuncValues
is keeping the values of the functions from the FUNCTIONS block