Dear Forum,
I am trying to save the results of a toy monte carlo into file, through a TTree with the following code
TFile * fOutTmp;
TTree * treeFitResult;
Double_t Nb_fit,Ns_fit,NbError_fit,NsError_fit;
while(N < Nmax)
{// loop over different values of N == number of signal injected events
cout << "MAIN LOOP: N = " << N << " Nmax = " << Nmax << " Nstep = " << Nstep << endl;
cout << hline << endl;
hB = new TH1D("hB","Log Likelihood Ratio - Background",50000,0.,5000);
hSB = new TH1D("hSB","Log Likelihood Ratio - Signal + Background",50000,0.,5000);
treeFitResult = new TTree("fit_results","fit_results");
treeFitResult->Branch("Nb_injected",&Nbkg,"Nb_injected/I");
treeFitResult->Branch("Ns_injected",&N,"Ns_injected/I");
treeFitResult->Branch("Nb_fit",&Nb_fit,"Nb_fit/D");
treeFitResult->Branch("Ns_fit",&Ns_fit,"Ns_fit/D");
treeFitResult->Branch("NbError_fit",&NbError_fit,"NbError_fit/D");
treeFitResult->Branch("NsError_fit",&NsError_fit,"NsError_fit/D");
for(int a=0; a<Ntoy; a++)
{
Nstatus ++;
if(Nstatus % 10 == 0)cout << "STATUS: done " << 100.*Nstatus/(double)Ntotali << " %" << endl;
h0 = new TH1D("h0","Background Only sample", BackgroundPDF->GetNbinsX() ,BackgroundPDF->GetBinLowEdge( 1 ) , BackgroundPDF->GetBinLowEdge( BackgroundPDF->GetNbinsX() + 1 ) );
h1 = new TH1D("h1","Background + Signal sample", BackgroundPDF->GetNbinsX() ,BackgroundPDF->GetBinLowEdge( 1 ) ,BackgroundPDF->GetBinLowEdge( BackgroundPDF->GetNbinsX() + 1 ) );
cout << "Created histograms for B and S+B." << endl;
cout << "Nbin = \t" << BackgroundPDF->GetNbinsX() << endl;
cout << "Xmin = \t" << BackgroundPDF->GetBinLowEdge( 1 ) << endl;
cout << "Xmax = \t" << BackgroundPDF->GetBinLowEdge( BackgroundPDF->GetNbinsX() + 1 ) << endl;
for(int i=0; i<Ngen; i++)
{
// fill histograms with background
h0->Fill( BackgroundPDF->GetRandom() );
h1->Fill( BackgroundPDF->GetRandom() );
}
for(int j=0; j<N; j++)// inject signal into h1
h1->Fill( SignalPDF->GetRandom() );
// fit the histograms in the BKG case with both models
// 1. --------------------------------------------------------
// BACKGROUND fit with BACKGROUND
cout << "***** FIT B with B" << endl;
fB.SetParameter(0,Nbkg); // setting seed for the fit
fResultB = h0->Fit(&fB,"QRSL");
// some rene brun magic to get likelihood at convergence (amin)
// amin is 2*(log likelihood) (see function H1FitLikelihood in class TH1)
pFitterB = TVirtualFitter::Fitter(h0);
pFitterB->GetStats(aminB,edmB,errdefB,nvparB,nparxB); // amin == minimum likelihood
// BACKGROUND fit with SIGNAL + BACKGROUND
cout << "***** FIT B with S+B" << endl;
fSB.SetParameter(0,Nbkg); // setting seed for the fit
fSB.SetParameter(1,N);
fResultSB = h0->Fit(&fSB,"QRSL");
pFitterSB = TVirtualFitter::Fitter(h0);
pFitterSB->GetStats(aminSB,edmSB,errdefSB,nvparSB,nparxSB); // amin == minimum likelihood
// compute likelihood ratio and add it to histo
LogLikelihoodRatio = aminB - aminSB;
hB->Fill( LogLikelihoodRatio );
cout << "LLR Background = " << LogLikelihoodRatio << endl;
// 2. --------------------------------------------------------
// SIGNAL+BACKGROUND fit with BACKGROUND
cout << "***** FIT S+B with B" << endl;
fB.SetParameter(0,Nbkg); // seed
fResultB = h1->Fit(&fB,"QRSL");
pFitterB = TVirtualFitter::Fitter(h1);
pFitterB->GetStats(aminB,edmB,errdefB,nvparB,nparxB);
cout << "aminB = " << aminB << endl;
// SIGNAL+BACKGROUND fit with SIGNAL + BACKGROUND
cout << "***** FIT S+B with S+B" << endl;
fSB.SetParameter(0,Nbkg); // seed
fSB.SetParameter(1,N);
fResultSB = h1->Fit(&fSB,"QRSL");
pFitterSB = TVirtualFitter::Fitter(h1);
pFitterSB->GetStats(aminSB,edmSB,errdefSB,nvparSB,nparxSB);
cout << "aminSB = " << aminSB << endl;
Nb_fit = pFitterSB->GetParameter(0);
Ns_fit = pFitterSB->GetParameter(1);
NbError_fit = pFitterSB->GetParError(0);
NsError_fit = pFitterSB->GetParError(1);
// compute likelihood ratio and add it to histo
LogLikelihoodRatio = aminB - aminSB;
hSB->Fill( LogLikelihoodRatio );
cout << a << " / " << Ntoy << "\t LLR Signal+Background = " << LogLikelihoodRatio<< endl;
cout << "Filling: Nb_injected = " << Nbkg << endl;
cout << "Filling: Ns_injected = " << N << endl;
cout << "Filling: Nb_fit = " << Nb_fit << endl;
cout << "Filling: Ns_fit = " << Ns_fit << endl;
cout << "Filling: NbError_fit = " << NbError_fit << endl;
cout << "Filling: NsError_fit = " << NsError_fit << endl;
cout << hline << endl;
treeFitResult->Fill();
delete h0;
delete h1;
}
// ONCE draw the distribution of the LLR and save them to file
cout << hline << endl;
if( SingleSimulation == false)cout << "INFO: SingleSimulation = false" << endl;
cout << "N \t = " << N << endl;
fOutTmp = new TFile( Form("%sLLR_Sblue_SBred_N%d.root",outputFolder.c_str(),N), "RECREATE" );
if( fOutTmp->IsOpen() ){
cout << "INFO: File " << Form("%sLLR_Sblue_SBred_N%d.root",outputFolder.c_str(),N) << " opened successfully." << endl;
}else{
cout << "INFO: ERROR opening " << Form("%sLLR_Sblue_SBred_N%d.root",outputFolder.c_str(),N) << endl;
}
// output to file
gDirectory->pwd();
hB->SetDirectory( gDirectory );
hB->Write("hB");
cout << "INFO: hB written to file." << endl;
hSB->SetDirectory( gDirectory );
hSB->Write("hSB");
cout << "INFO: hSB written to file." << endl;
treeFitResult->Show(0);
treeFitResult->SetDirectory( gDirectory );
treeFitResult->Write("treeFitResult");
cout << "INFO: treeFitResult has " << treeFitResult->GetEntries() << " entries." << endl;
cout << "INFO: treeFitResult written to file." << endl;
fOutTmp->Close();
delete fOutTmp;
cout << "INFO: file closed." << endl;
if(Nstep < 0) break; // exit from loop
// DEBUG: Load this into independent macro that makes plots for sensitivity
// compute median of S+B histogram / B histogram
hSB->GetQuantiles(1,&MedianSB,&qSB);
hB->GetQuantiles(1,&MedianB,&qSB);
cout << "INFO: GetQuantiles done." << endl;
// compute p-value of B histogram from median of SB to infty
pValue = hB->Integral( hB->GetXaxis()->FindBin( MedianSB ), hB->GetNbinsX()-1 );
if( pValue < 1. )pValue = 1.;
pValue /= hB->Integral( 0, hB->GetNbinsX()-1 );
cout << "INFO: N \t\t MedianSB \t MedianB \t pValue" << endl;
cout << "INFO: " << N << "\t " << Form("%.1lf",MedianSB) << "\t " << Form("%.1lf",MedianB) << "\t " << pValue << endl;
// FILL A TGRAPH AND OUTPUT IT
pValue_vs_N->SetPoint(counter,GetHalfLife(N),pValue);
counter++;
// clear memory
delete hB;
delete hSB;
N += Nstep;
}
In the log file, everything seems fine. The variables linked to the TTree branches are correctly filled and the ->Show(0) gives the correct output.
But when I open the root file after running this, I get an error:
[gfantini@cuore-login1 output]$ root LLR_Sblue_SBred_N97.root
*******************************************
* *
* W E L C O M E to R O O T *
* *
* Version 5.34/34 2 October 2015 *
* *
* You are welcome to visit our Web site *
* http://root.cern.ch *
* *
*******************************************
ROOT 5.34/34 (v5-34-34@v5-34-34, Oct 02 2015, 16:30:37 on linuxx8664gcc)
CINT/ROOT C/C++ Interpreter version 5.18.00, July 2, 2010
Type ? for help. Commands must be C++ statements.
Enclose multiple statements between { }.
root [0]
Attaching file LLR_Sblue_SBred_N97.root as _file0...
root [1] .ls
TFile** LLR_Sblue_SBred_N97.root
TFile* LLR_Sblue_SBred_N97.root
KEY: TH1D hB;1 Log Likelihood Ratio - Background
KEY: TH1D hSB;1 Log Likelihood Ratio - Signal + Background
KEY: TTree treeFitResult;1 fit_results
root [2] treeFitResult->Show(0)
Error in <TFile::ReadBuffer>: error reading all requested bytes from file LLR_Sblue_SBred_N97.root, got 0 of 277
Error in <TBranch::GetBasket>: File: LLR_Sblue_SBred_N97.root at byte:0, branch:Nb_injected, entry:0, badread=1, nerrors=1, basketnumber=0
Error in <TTree::Show()>: Cannot read entry 0 (I/O error)
root [3]
If anyone has an idea on what can I do to debug this is very much welcome.
I’m at the point of trying to awk-parse the log file to get the information i need!
Cheers
Guido
ROOT Version: 5.34/34
Platform: Not Provided
Compiler: g++ (GCC) 7.1.0