diff --git a/src/classes/PMsrHandler.cpp b/src/classes/PMsrHandler.cpp index 17ab8782..4189f6ca 100644 --- a/src/classes/PMsrHandler.cpp +++ b/src/classes/PMsrHandler.cpp @@ -2673,7 +2673,8 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) tokens = 0; } } else if (iter1->fLine.Contains("rrf_phase", TString::kIgnoreCase)) { - // expected entry: rrf_phase value. value given in units of degree + // expected entry: rrf_phase value. value given in units of degree. or + // rrf_phase parX. where X is the parameter number, e.g. par3 tokens = iter1->fLine.Tokenize(" \t"); if (!tokens) { cerr << endl << ">> PMsrHandler::HandlePlotEntry: **SEVERE ERROR** Couldn't tokenize rrf_phase in line " << iter1->fLineNo; @@ -2689,7 +2690,19 @@ Bool_t PMsrHandler::HandlePlotEntry(PMsrLines &lines) if (str.IsFloat()) { param.fRRFPhase = str.Atof(); } else { - error = true; + if (str.BeginsWith("par", TString::kIgnoreCase)) { // parameter value + Int_t no = 0; + if (FilterNumber(str, "par", 0, no)) { + // check that the parameter is in range + if ((Int_t)fParam.size() < no) { + error = true; + } else { + param.fRRFPhase = fParam[no-1].fValue; + } + } + } else { + error = true; + } } } // clean up diff --git a/src/classes/PRunAsymmetry.cpp b/src/classes/PRunAsymmetry.cpp index 8e2703d6..3fc5ddbb 100644 --- a/src/classes/PRunAsymmetry.cpp +++ b/src/classes/PRunAsymmetry.cpp @@ -204,7 +204,7 @@ Double_t PRunAsymmetry::CalcChiSquare(const std::vector& par) */ Double_t PRunAsymmetry::CalcMaxLikelihood(const std::vector& par) { - cout << endl << "PRunSingleHisto::CalcMaxLikelihood(): not implemented yet ..." << endl; + cout << endl << "PRunAsymmetry::CalcMaxLikelihood(): not implemented yet ..." << endl; return 1.0; } @@ -496,7 +496,10 @@ Bool_t PRunAsymmetry::PrepareData() status = PrepareFitData(runData, histoNo); break; case kView: - status = PrepareViewData(runData, histoNo); + if (fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking == 0) + status = PrepareViewData(runData, histoNo); + else + status = PrepareRRFViewData(runData, histoNo); break; default: status = false; @@ -657,25 +660,25 @@ Bool_t PRunAsymmetry::PrepareFitData(PRawRunData* runData, UInt_t histoNo[2]) // check if data range has been provided, and if not try to estimate them if (start[0] < 0) { start[0] = (Int_t)t0[0]+5; - cerr << endl << "PRunSingleHisto::PrepareData(): **WARNING** data range (forward) was not provided, will try data range start = t0+5 = " << start[0] << "."; + cerr << endl << "PRunAsymmetry::PrepareData(): **WARNING** data range (forward) was not provided, will try data range start = t0+5 = " << start[0] << "."; cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; cerr << endl; } if (start[1] < 0) { start[1] = (Int_t)t0[1]+5; - cerr << endl << "PRunSingleHisto::PrepareData(): **WARNING** data range (backward) was not provided, will try data range start = t0+5 = " << start[1] << "."; + cerr << endl << "PRunAsymmetry::PrepareData(): **WARNING** data range (backward) was not provided, will try data range start = t0+5 = " << start[1] << "."; cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; cerr << endl; } if (end[0] < 0) { end[0] = runData->GetDataBin(histoNo[0])->size(); - cerr << endl << "PRunSingleHisto::PrepareData(): **WARNING** data range (forward) was not provided, will try data range end = " << end[0] << "."; + cerr << endl << "PRunAsymmetry::PrepareData(): **WARNING** data range (forward) was not provided, will try data range end = " << end[0] << "."; cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; cerr << endl; } if (end[1] < 0) { end[1] = runData->GetDataBin(histoNo[1])->size(); - cerr << endl << "PRunSingleHisto::PrepareData(): **WARNING** data range (backward) was not provided, will try data range end = " << end[1] << "."; + cerr << endl << "PRunAsymmetry::PrepareData(): **WARNING** data range (backward) was not provided, will try data range end = " << end[1] << "."; cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; cerr << endl; } @@ -839,12 +842,12 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) // check if data range has been provided, and if not try to estimate them if (fRunInfo->GetDataRange(0) < 0) { start[0] = ((Int_t)t0[0]+5) - (((Int_t)t0[0]+5)/packing)*packing; - cerr << endl << "PRunSingleHisto::PrepareData(): **WARNING** data range (forward) was not provided, will try data range start = " << start[0] << "."; + cerr << endl << "PRunAsymmetry::PrepareViewData(): **WARNING** data range (forward) was not provided, will try data range start = " << start[0] << "."; cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; cerr << endl; } else if (fRunInfo->GetDataRange(2) < 0) { start[1] = ((Int_t)t0[1]+5) - (((Int_t)t0[1]+5)/packing)*packing; - cerr << endl << "PRunSingleHisto::PrepareData(): **WARNING** data range (backward) was not provided, will try data range start = " << start[1] << "."; + cerr << endl << "PRunAsymmetry::PrepareViewData(): **WARNING** data range (backward) was not provided, will try data range start = " << start[1] << "."; cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; cerr << endl; } else { @@ -858,12 +861,6 @@ Bool_t PRunAsymmetry::PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]) start[1] = val + fRunInfo->GetDataRange(2) - fRunInfo->GetDataRange(0); } -/* -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 UInt_t noOfBins0 = (runData->GetDataBin(histoNo[0])->size()-start[0])/packing; UInt_t noOfBins1 = (runData->GetDataBin(histoNo[1])->size()-start[1])/packing; @@ -923,9 +920,9 @@ cout << endl; value = 0.0; error = 0.0; } + value += fForward[i]; + error += fForwardErr[i]*fForwardErr[i]; } - value += fForward[i]; - error += fForwardErr[i]*fForwardErr[i]; } // backward for (Int_t i=start[1]; i Prepares the RRF data set for visual representation. This is done the following way: + * -# make all necessary checks + * -# build the asymmetry [ \f[ A(t) \f] ] WITHOUT packing. + * -# \f[ A_R(t) = A(t) \cdot 2 \cos(\omega_R t + \phi_R) \f] + * -# do the packing of \f[ A_R(t) \f] + * -# calculate theory [ \f[ T(t) \f] ] as close as possible to the time resolution [compatible with the RRF frequency] + * -# \f[ T_R(t) = T(t) \cdot 2 \cos(\omega_R t + \phi_R) \f] + * -# do the packing of \f[ T_R(t) \f] + * -# calculate the Kaiser FIR filter coefficients + * -# filter \f[ T_R(t) \f]. + * + * \param runData raw run data needed to perform some crosschecks + * \param histNo array of the histo numbers form which to build the asymmetry + */ +Bool_t PRunAsymmetry::PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]) +{ + // feed the parameter vector + std::vector par; + PMsrParamList *paramList = fMsrInfo->GetMsrParamList(); + for (UInt_t i=0; isize(); i++) + par.push_back((*paramList)[i].fValue); + + // ------------------------------------------------------------ + // 1. make all necessary checks + // ------------------------------------------------------------ + + // first get start data, end data, and t0 + Int_t start[2] = {0, 0}; + Int_t end[2] = {0, 0}; + Double_t t0[2] = {fT0s[0], fT0s[1]}; + UInt_t packing = fMsrInfo->GetMsrPlotList()->at(0).fRRFPacking; + // check if data range has been provided, and if not try to estimate them + if (fRunInfo->GetDataRange(0) < 0) { + start[0] = static_cast(t0[0])+5; + cerr << endl << "PRunAsymmetry::PrepareRRFViewData(): **WARNING** data range (forward) was not provided, will try data range start = " << start[0] << "."; + cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } else if (fRunInfo->GetDataRange(2) < 0) { + start[1] = static_cast(t0[1])+5; + cerr << endl << "PRunAsymmetry::PrepareRRFViewData(): **WARNING** data range (backward) was not provided, will try data range start = " << start[1] << "."; + cerr << endl << "NO WARRANTY THAT THIS DOES MAKE ANY SENSE."; + cerr << endl; + } else { + Int_t val = fRunInfo->GetDataRange(0)-packing*(fRunInfo->GetDataRange(0)/packing); + do { + if (fRunInfo->GetDataRange(2) - fRunInfo->GetDataRange(0) < 0) + val += packing; + } while (val + fRunInfo->GetDataRange(2) - fRunInfo->GetDataRange(0) < 0); + + start[0] = val; + start[1] = val + fRunInfo->GetDataRange(2) - fRunInfo->GetDataRange(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]; + UInt_t noOfBins1 = runData->GetDataBin(histoNo[1])->size()-start[1]; + if (noOfBins0 > noOfBins1) + noOfBins0 = noOfBins1; + 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 + if ((start[i] < 0) || (start[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << "PRunAsymmetry::PrepareRRFViewData(): **ERROR** start data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 3rd check if end is within proper bounds + if ((end[i] < 0) || (end[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << "PRunAsymmetry::PrepareRRFViewData(): **ERROR** end data bin doesn't make any sense!"; + cerr << endl; + return false; + } + // 4th check if t0 is within proper bounds + if ((t0[i] < 0) || (t0[i] > (Int_t)runData->GetDataBin(histoNo[i])->size())) { + cerr << endl << "PRunAsymmetry::PrepareRRFViewData(): **ERROR** t0 data bin doesn't make any sense!"; + cerr << endl; + return false; + } + } + + // ------------------------------------------------------------ + // 2. build the asymmetry [A(t)] WITHOUT packing. + // ------------------------------------------------------------ + + PDoubleVector forward, forwardErr; + PDoubleVector backward, backwardErr; + Double_t error = 0.0; + // forward + for (Int_t i=start[0]; i backward.size()) + noOfBins = backward.size(); + } + + // form asymmetry including error propagation + PDoubleVector asymmetry, asymmetryErr; + Double_t asym; + Double_t f, b, ef, eb, alpha = 1.0, beta = 1.0; + + // get the proper alpha and beta + switch (fAlphaBetaTag) { + case 1: // alpha == 1, beta == 1 + alpha = 1.0; + beta = 1.0; + break; + case 2: // alpha != 1, beta == 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = 1.0; + break; + case 3: // alpha == 1, beta != 1 + alpha = 1.0; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + case 4: // alpha != 1, beta != 1 + alpha = par[fRunInfo->GetAlphaParamNo()-1]; + beta = par[fRunInfo->GetBetaParamNo()-1]; + break; + default: + break; + } +//cout << endl << ">> alpha = " << alpha << ", beta = " << beta; + + for (UInt_t i=0; iGetMsrPlotList()->at(0).fRRFUnit) { + case RRF_UNIT_kHz: + gammaRRF = TMath::TwoPi()*1.0e-3; + break; + case RRF_UNIT_MHz: + gammaRRF = TMath::TwoPi(); + break; + case RRF_UNIT_Mcs: + gammaRRF = 1.0; + break; + case RRF_UNIT_G: + gammaRRF = 0.0135538817*TMath::TwoPi(); + break; + case RRF_UNIT_T: + gammaRRF = 0.0135538817*TMath::TwoPi()*1.0e4; + break; + default: + gammaRRF = TMath::TwoPi(); + break; + } + wRRF = gammaRRF * fMsrInfo->GetMsrPlotList()->at(0).fRRFFreq; + phaseRRF = fMsrInfo->GetMsrPlotList()->at(0).fRRFPhase / 180.0 * TMath::Pi(); + + for (UInt_t i=0; i(start[0])-t0[0]+static_cast(i)); + asymmetry[i] *= 2.0*TMath::Cos(wRRF*time+phaseRRF); + } + + // ------------------------------------------------------------ + // 4. do the packing of A_R(t) + // ------------------------------------------------------------ + Double_t value = 0.0; + error = 0.0; + for (UInt_t i=0; i(start[0])-t0[0]+static_cast(packing-1)/2.0)); + fData.SetDataTimeStep(fTimeResolution*static_cast(packing)); + + // ------------------------------------------------------------ + // 5. calculate theory [T(t)] as close as possible to the time resolution [compatible with the RRF frequency] + // 6. T_R(t) = T(t) * 2 cos(w_R t + phi_R) + // ------------------------------------------------------------ + UInt_t rebinRRF = static_cast((TMath::Pi()/2.0/wRRF)/fTimeResolution); // RRF time resolution / data time resolution + fData.SetTheoryTimeStart(fData.GetDataTimeStart()); + fData.SetTheoryTimeStep(TMath::Pi()/2.0/wRRF/rebinRRF); // = theory time resolution as close as possible to the data time resolution compatible with wRRF + +/* +cout << endl << ">> rebinRRF = " << rebinRRF; +cout << endl << ">> theory time start = " << fData.GetTheoryTimeStart(); +cout << endl << ">> theory time step = " << fData.GetTheoryTimeStep(); +cout << endl; +*/ + + // calculate functions + for (Int_t i=0; iGetNoOfFuncs(); i++) { + fFuncValues[i] = fMsrInfo->EvalFunc(fMsrInfo->GetFuncNo(i), *fRunInfo->GetMap(), par); + } + + Double_t theoryValue; + for (UInt_t i=0; i(i)*fData.GetTheoryTimeStep(); + theoryValue = fTheory->Func(time, par, fFuncValues); + theoryValue *= 2.0*TMath::Cos(wRRF * time + phaseRRF); + + if (fabs(theoryValue) > 10.0) { // dirty hack needs to be fixed!! + theoryValue = 0.0; + } + + fData.AppendTheoryValue(theoryValue); + } + + // ------------------------------------------------------------ + // 7. do the packing of T_R(t) + // ------------------------------------------------------------ + + PDoubleVector theo; + Double_t dval = 0.0; + for (UInt_t i=0; isize(); i++) { + if ((i % rebinRRF == 0) && (i != 0)) { +//cout << endl << "time = " << fData.GetTheoryTimeStart() + i * fData.GetTheoryTimeStep() << ", theory value = " << dval; + theo.push_back(dval/rebinRRF); + dval = 0.0; + } + dval += fData.GetTheory()->at(i); + } + fData.ReplaceTheory(theo); + theo.clear(); + + // set the theory time start and step size + fData.SetTheoryTimeStart(fData.GetTheoryTimeStart()+static_cast(rebinRRF-1)*fData.GetTheoryTimeStep()/2.0); + fData.SetTheoryTimeStep(rebinRRF*fData.GetTheoryTimeStep()); + + // ------------------------------------------------------------ + // 8. calculate the Kaiser FIR filter coefficients + // ------------------------------------------------------------ + CalculateKaiserFilterCoeff(wRRF, 60.0, 0.2); // w_c = wRRF, A = -20 log_10(delta), Delta w / w_c = (w_s - w_p) / (2 w_c) + + // ------------------------------------------------------------ + // 9. filter T_R(t) + // ------------------------------------------------------------ + FilterTheo(); + + // clean up + par.clear(); + forward.clear(); + forwardErr.clear(); + backward.clear(); + backwardErr.clear(); + asymmetry.clear(); + asymmetryErr.clear(); + + return true; +} diff --git a/src/classes/PRunBase.cpp b/src/classes/PRunBase.cpp index e8e9ac3a..37bb1db2 100644 --- a/src/classes/PRunBase.cpp +++ b/src/classes/PRunBase.cpp @@ -200,12 +200,12 @@ cout << endl; } //-------------------------------------------------------------------------- -// FilterData (protected) +// FilterTheo (protected) //-------------------------------------------------------------------------- /** - *

Filters the data with a Kaiser FIR filter. + *

Filters the theory with a Kaiser FIR filter. */ -void PRunBase::FilterData() +void PRunBase::FilterTheo() { PDoubleVector theoFiltered; Double_t dval = 0.0; diff --git a/src/classes/PRunSingleHisto.cpp b/src/classes/PRunSingleHisto.cpp index 5634f294..7594fb67 100644 --- a/src/classes/PRunSingleHisto.cpp +++ b/src/classes/PRunSingleHisto.cpp @@ -1002,9 +1002,9 @@ cout << endl << "--------------------------------" << endl; theo.clear(); } - // filter data + // filter theory CalculateKaiserFilterCoeff(wRRF, 60.0, 0.2); // w_c = wRRF, A = -20 log_10(delta), Delta w / w_c = (w_s - w_p) / (2 w_c) - FilterData(); + FilterTheo(); } // clean up diff --git a/src/include/PRunAsymmetry.h b/src/include/PRunAsymmetry.h index 19719fdf..b3d81ab3 100644 --- a/src/include/PRunAsymmetry.h +++ b/src/include/PRunAsymmetry.h @@ -51,6 +51,7 @@ class PRunAsymmetry : public PRunBase virtual Bool_t PrepareData(); virtual Bool_t PrepareFitData(PRawRunData* runData, UInt_t histoNo[2]); virtual Bool_t PrepareViewData(PRawRunData* runData, UInt_t histoNo[2]); + virtual Bool_t PrepareRRFViewData(PRawRunData* runData, UInt_t histoNo[2]); private: UInt_t fAlphaBetaTag; ///< 1-> alpha = beta = 1; 2-> alpha != 1, beta = 1; 3-> alpha = 1, beta != 1; 4-> alpha != 1, beta != 1 diff --git a/src/include/PRunBase.h b/src/include/PRunBase.h index 87113656..78055fda 100644 --- a/src/include/PRunBase.h +++ b/src/include/PRunBase.h @@ -87,7 +87,7 @@ class PRunBase virtual Bool_t PrepareData() = 0; // pure virtual, i.e. needs to be implemented by the deriving class!! virtual void CalculateKaiserFilterCoeff(Double_t wc, Double_t A, Double_t dw); - virtual void FilterData(); + virtual void FilterTheo(); }; #endif // _PRUNBASE_H_ diff --git a/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C b/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C index f596fa9c..270395fc 100644 --- a/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C +++ b/src/tests/t0NotEqFirstGoodData/t0NotEqFirstGoodData.C @@ -44,19 +44,24 @@ void t0NotEqFirstGoodData() decayAnaModule = histosFolder->AddFolder("DecayAnaModule", "muSR decay histograms"); // feed run info header + UInt_t runNo = 10113; + TString tstr; runInfo = gROOT->GetRootFolder()->AddFolder("RunInfo", "LEM RunInfo"); gROOT->GetListOfBrowsables()->Add(runInfo, "RunInfo"); header = new TLemRunHeader(); - header->SetRunTitle("010101 - test"); + tstr = TString("0"); + tstr += runNo; + tstr += TString(" - test"); + header->SetRunTitle(tstr.Data()); header->SetLemSetup("trivial"); - header->SetRunNumber(10101); + header->SetRunNumber(runNo); header->SetStartTime(0); header->SetStopTime(1); header->SetModeratorHV(32.0, 0.01); header->SetSampleHV(0.0, 0.01); header->SetImpEnergy(31.8); header->SetSampleTemperature(0.2, 0.001); - header->SetSampleBField(30000, 0.1); + header->SetSampleBField(-1.0, 0.1); header->SetTimeResolution(0.1953125); header->SetNChannels(66601); header->SetNHist(4); @@ -68,20 +73,19 @@ void t0NotEqFirstGoodData() // create histos UInt_t t0[4] = {3419, 3419, 3419, 3419}; - UInt_t n0[4] = {2000.0, 2005.0, 2003.0, 1998.0}; + UInt_t n0[4] = {200.0, 205.0, 203.0, 198.0}; UInt_t bkg[4] = {11.0, 11.0, 5.0, 8.0}; const Double_t tau = 2197.019; // ns // asymmetry related stuff const Double_t timeResolution = 0.1953125; // ns - Double_t a0[4] = {0.16, 0.16, 0.16, 0.16}; + Double_t a0[4] = {0.26, 0.26, 0.26, 0.26}; Double_t a1[4] = {0.0, 0.0, 0.0, 0.0}; Double_t phase[4] = {0.0*TMath::Pi()/180.0, 90.0*TMath::Pi()/180.0, 180.0*TMath::Pi()/180.0, 270.0*TMath::Pi()/180.0}; const Double_t gamma = 0.00001355; // gamma/(2pi) Double_t bb0 = 15000.0; // field in Gauss - Double_t bb1 = 15702.0; // field in Gauss -// Double_t bb = 100.0; // field in Gauss + Double_t bb1 = 400.0; // field in Gauss Double_t sigma0 = 0.05/1000.0; // Gaussian sigma in 1/ns - Double_t sigma1 = 0.10/1000.0; // Gaussian sigma in 1/ns + Double_t sigma1 = 0.005/1000.0; // Gaussian sigma in 1/ns TH1F *histo[8]; char str[128]; @@ -99,7 +103,7 @@ void t0NotEqFirstGoodData() histo[i]->SetBinContent(j, bkg[i]); } else { time = (Double_t)(j-t0[i])*timeResolution; - dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-0.5*TMath::Power(time*sigma0,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]) + dval = (Double_t)n0[i]*TMath::Exp(-time/tau)*(1.0+a0[i]*TMath::Exp(-0.5*TMath::Power(time*sigma0,1.2))*TMath::Cos(TMath::TwoPi()*gamma*bb0*time+phase[i]) +a1[i]*TMath::Exp(-0.5*TMath::Power(time*sigma1,2.0))*TMath::Cos(TMath::TwoPi()*gamma*bb1*time+phase[i]))+(Double_t)bkg[i]; histo[i]->SetBinContent(j, (UInt_t)dval); } @@ -107,6 +111,7 @@ void t0NotEqFirstGoodData() } // add a promp peak +/* Double_t ampl = 500.0; Double_t promptLambda = 500.0/1000.0; // Lorentzain in 1/ns for (UInt_t i=0; i<4; i++) { @@ -118,6 +123,7 @@ void t0NotEqFirstGoodData() histo[i+4]->SetBinContent(j, (UInt_t)dval); } } +*/ // add Poisson noise PAddPoissonNoise *addNoise = new PAddPoissonNoise(); @@ -138,7 +144,10 @@ void t0NotEqFirstGoodData() decayAnaModule->Add(histo[i]); // write file - TFile *fout = new TFile("010108.root", "RECREATE", "Midas Fake Histograms"); + tstr = TString("0"); + tstr += runNo; + tstr += TString(".root"); + TFile *fout = new TFile(tstr.Data(), "RECREATE", "Midas Fake Histograms"); if (fout == 0) { cout << endl << "**ERROR** Couldn't create ROOT file"; cout << endl << endl;