Dear all,
I a piece of code that should be plotting several TGraphs on the same canvas, but doesn’t. Any help would be very much appreciated.
What the code does: in main function, it calls Calibrator::Fit_TH1D which follows the steps below
[ul]
[li]Loop over E bins of THnSparse --> TH2D projected histogram[/li]
[li]Per bin: take TH1D slice and fit Gaussian (if number of events is above threshold)[/li]
[li]Store E value of bin and mean of Gaussian fit, and create TGraph[/li][/ul]
This TGraph can be drawn inside of the Calibrator.cpp code, it can be passed by reference and its properties (axis title, …) can be accessed in the main function. The lines
[quote] cout << “Graph\t” << graph_meas->GetMaximum() << endl;
cout << “Graph\t” << graph_meas->GetMinimum() << endl;
cout << “Graph\t” << graph_meas->GetXaxis()->GetTitle() << endl;[/quote]
return the same values in the Fit_TH1D and the main function, but when I try to draw the graph in the main function, things go awry.
Canvas Test_sector_calibration_%i.C should contain graph graph_meas, but is empty. When I open the macro, I can see that TH1D *Response_1D_bin_49_phi_4 is created there: the histogram is empty, but the name derives from hResponse_1D, which is the fitted histogram halfway through the Fit_TH1D function.
I believe that something is going wrong with the fit, but I have been grinding on this code for days and I am missing something.
MAIN
TFile *_file0 = TFile::Open( filenames[file], "Read");
Calibrator my_first_calibrator(filenames[file], "20150416");
TCanvas *canvas = new TCanvas("canvas", "canvas", 800,800);
TString drawoptions = "ape";
for(int phi_bin = 1; phi_bin <= 4; phi_bin++){
TGraphErrors* graph_true;
TGraphErrors* graph_meas;
cout << "// -- SECTOR\t" << phi_bin << endl;
my_first_calibrator.Fit_TH1D( phi_bin, 50, TString::Format("%i", phi_bin), graph_true, graph_meas);
graph_meas->SetLineColor( getColor(phi_bin) );
graph_meas->SetMarkerColor( getColor(phi_bin) );
graph_meas->Draw( "ape" );
cout << "Graph\t" << graph_meas->GetMaximum() << endl;
cout << "Graph\t" << graph_meas->GetMinimum() << endl;
cout << "Graph\t" << graph_meas->GetXaxis()->GetTitle() << endl;
canvas->SaveAs(TString::Format("Test_sector_calibration_%i.C", phi_bin));
drawoptions = "psame";
}
canvas->SaveAs("Test_sector_calibration.C");
Calibrator::Fit_TH1D
void Calibrator::Fit_TH1D(int phi_bin, int bins, TString label, TGraphErrors* &graph_true, TGraphErrors* &graph_meas){
vector<double> mean_response;
vector<double> mean_response_inverse;
vector<double> mean_eDet;
vector<double> error_eDet;
vector<double> mean_eGen;
vector<double> error_eGen;
vector<double> eGen_center;
vector<double> error_energy;
TH2D* hDetE_selection; Extract_TH2D( phi_bin, hDetE_selection, hDetE_fine_phi_);
TH2D* hGenE_selection; Extract_TH2D( phi_bin, hGenE_selection, hGenE_fine_phi_);
TH2D* hResponse_selection; Extract_TH2D( phi_bin, hResponse_selection, hResponse_);
int rebinner = rebinner_;
int event_threshold_ = 200;
TH1D* hSlice_storage;
TH1D* hEgen_storage;
TH1D* hEdet_storage;
StorageHistograms( hSlice_storage, hEgen_storage, hEdet_storage, hResponse_selection, hGenE_selection, hDetE_selection);
for(int bin_E = 1; bin_E < hResponse_selection->GetNbinsX() ; bin_E++){
// Get 1D projections.
TH1D* hResponse_1D = (TH1D*)hResponse_selection->ProjectionY(TString::Format("Response_1D_bin_%i_phi_%i", bin_E, phi_bin), bin_E, bin_E, "do");
TH1D* hGenE_1D = (TH1D*)hGenE_selection->ProjectionY(TString::Format("eGen_1D_bin_%i_phi_%i", bin_E, phi_bin), bin_E, bin_E, "o");
hGenE_1D->Rebin( rebinner*100 );
TH1D* hDetE_1D = (TH1D*)hDetE_selection->ProjectionY(TString::Format("eDet_1D_bin_%i_phi_%i", bin_E, phi_bin), bin_E, bin_E, "o");
hDetE_1D->Rebin( rebinner*100 );
// Create gaussian.
TF1 *fit_response = new TF1(TString::Format("Fit_Response_1D_bin_%i_phi_%i", bin_E, phi_bin), "gaus");
if( (hResponse_1D->Integral() + hSlice_storage->Integral() )< event_threshold_ ){
hSlice_storage ->Add( hResponse_1D );
hEgen_storage ->Add( hGenE_1D );
hEdet_storage ->Add( hDetE_1D );
continue;
}
else if( ( hResponse_1D->Integral() + hSlice_storage->Integral() ) >= event_threshold_ ){
hResponse_1D ->Add( hSlice_storage );
hGenE_1D ->Add( hEgen_storage );
hDetE_1D ->Add( hEdet_storage );
}
for(int bins_store = 0; bins_store <= hSlice_storage->GetNbinsX(); bins_store++){ hSlice_storage->SetBinContent(bins_store, 0.); hSlice_storage->SetBinError(bins_store, 0.);}
for(int bins_store = 0; bins_store <= hEgen_storage->GetNbinsX(); bins_store++){ hEgen_storage->SetBinContent(bins_store, 0.); hEgen_storage->SetBinError(bins_store, 0.);}
for(int bins_store = 0; bins_store <= hEdet_storage->GetNbinsX(); bins_store++){ hEdet_storage->SetBinContent(bins_store, 0.); hEdet_storage->SetBinError(bins_store, 0.);}
hResponse_1D ->Fit( fit_response , "Q0");
if( fit_response->GetParameter(1) != 0. ){
double R = fit_response->GetParameter(1);
double sR= fit_response->GetParError(1);
mean_response .push_back( R );
mean_response_inverse.push_back( 1./R );
mean_eDet .push_back( GetAverage(hDetE_1D) );
mean_eGen .push_back( GetAverage(hGenE_1D) );
eGen_center .push_back( hGen->GetBinCenter( bin_E ) );
cout << "//--\t" << GetAverage(hGenE_1D) << "\t" << R << "\t\t" << GetAverage(hDetE_1D) << "\t" << 1./R << endl;
if( R == 0.){ error_eGen .push_back( 0. ); }
else{ error_eGen .push_back( sR ); }
if( R == 0.){ error_eDet .push_back( 0. ); }
else{ error_eDet .push_back( sR/(R*R) ); }
error_energy .push_back( 0. );
}
} // Loop over bins.
// -- Transform vectors into arrays.
double * E_gen_axis = &mean_eGen[0],
* E_det_axis = &mean_eDet[0],
* response_val = &mean_response[0],
* response_inverse= &mean_response_inverse[0],
* eGen_bin = &eGen_center[0],
* err_R = &error_eGen[0],
* err_1overR = &error_eDet[0],
* err_energy = &error_energy[0];
graph_meas = new TGraphErrors( mean_eDet.size(), E_det_axis, response_inverse, err_energy, err_1overR);
graph_meas ->GetXaxis()->SetTitle("E_{det}");
graph_meas ->GetXaxis()->SetTitleOffset(0.9);
graph_meas ->GetXaxis()->SetTitleSize(0.06);
graph_meas ->GetXaxis()->SetLabelSize(0.05);
graph_meas ->GetXaxis()->SetNdivisions(505);
graph_meas ->GetYaxis()->SetTitle("< #frac{E_{gen}}{E_{det}} >");
graph_meas ->GetYaxis()->SetRangeUser(0.,2.5);
graph_meas ->GetYaxis()->SetTitleOffset(1.8);
graph_meas ->GetYaxis()->SetTitleSize(0.06);
graph_meas ->GetYaxis()->SetLabelSize(0.07);
graph_meas ->SetMarkerColor( getColor(phi_bin) );
graph_meas ->SetLineColor( getColor(phi_bin) );
graph_meas ->SetMarkerStyle( 20 + phi_bin );
graph_true ->SetName( TString::Format("Meas_response_sector_%i", phi_bin) );
TF1* analytical = new TF1("Analytical", " ( [0] + [1] * log( [2] + x) ) ", 0., 900.);
TCanvas *can_meas = new TCanvas( "can_meas_" + label, "can_meas_" + label, 1000, 1000);
can_meas ->SetLeftMargin(0.25);
can_meas ->SetTopMargin(0.07);
can_meas ->SetBottomMargin(0.14);
graph_meas->Draw("ape");
// graph_meas->Fit( analytical, "");
cout << "Graph\t" << graph_meas->GetMaximum() << endl;
cout << "Graph\t" << graph_meas->GetMinimum() << endl;
cout << "Graph\t" << graph_meas->GetXaxis()->GetTitle() << endl;
// hDetE_distr->Draw("histsame");
can_meas->SaveAs(folder_ + "/can_meas_" + label + ".C");
can_true->SaveAs(folder_ + "/can_true_" + label + ".C");
hResponse_selection->~TH2D();
hDetE_selection->~TH2D();
hGenE_selection->~TH2D();
hSlice_storage->~TH1D();
hEdet_storage->~TH1D();
hEgen_storage->~TH1D();
}
Calibrator::Extract_TH2D
// -- Prepare the TH2D
void Calibrator::Extract_TH2D(int phi_bin, TH2D* &hist_, THnSparseD* &hist_sparse_){
cout << "// -- EXTRACT 2D FROM THN" << endl;
hist_sparse_->GetAxis(2)->SetRange(phi_bin, phi_bin);
hist_ = (TH2D*)hist_sparse_->Projection(1,0);
hist_->RebinX( rebinner_ );
hist_->RebinY( rebinner_ );
}