Dear experts,
Is there a way to calculate the standard deviation of each bin from a TProfile histogram and then store the values as a new TProfile histogram?
Regards,
Soumyadip
Hi Soumyadip,
Thanks for the post.
You can get the error on the mean of the TProfile bin with the TProfile::GetBinError method. As for your second question, I am not sure I understand what you want to achieve: could you perhaps elaborate on that?
Best,
Danilo
Hi Danilo,
Thanks for your reply.
I have a 2D histogram. I have taken 1D projection with ProfileX to get a TProfile histogram to plot the average of a particular observable.
Now, my question is – is it possible to plot the standard deviation from that 2D histogram so that the standard deviation will be plotted on the y-axis?
I am attaching the picture.
From the 2D histogram, I have got a plot like the upper one where the y-axis represents the mean or average value.
Now, I am trying to get a plot like the bottom one. Where the y-axis is the standard deviation.
How can I plot something like that from the input 2D histogram?
Regards,
Soumyadip
Hi Soumyadip,
Thanks: I think I understand (please correct me if I am wrong).
Isn’t looping over TProfile bins, extract the deviations and replot them as a TGraph?
Cheers,
D
1 Like
Hi,
I have tried something like this but it didn’t worked
for (int id = 0; id < ndef; id++) {
for (int ij = 0; ij < njet; ij++) {
for(int ik=0; ik<nkappa; ik++){
TProfile* Test = (TProfile*)Data_Reco_JCO_PT[id][ij][ik]->Clone();
int ibins = Test->GetNbinsX();
double* xValues = new double[ibins];
double* binErrors = new double[ibins];
for(int ibin=1; ibin<=ibins; ++ibin){
xValues[ibin - 1] = Test->GetBinCenter(ibin);
binErrors[ibin - 1] =Test->GetBinError(ibin);
}
TGraph* graph = new TGraph(ibins, xValues, binErrors);
graph->SetTitle(Form("Bin Errors from Profile [%d][%d][%d];X;Error", id, ij, ik));
graph->SetMarkerStyle(20);
graph->SetMarkerColor(kBlue);
TCanvas* can24 = new TCanvas(Form("canvas_%d_%d_%d", id, ij, ik), "Bin Errors from Profiles", 800, 600);
graph->Draw(); // Draw with axis and points
// Save the canvas
if(id==0 && ij==0 && ik==0){can24->Print("Test_stddev.pdf(","pdf");}
else if(id ==2 && ij==1 && ik==9) {can24->Print("Test_stddev.pdf)","pdf");}
else{can24->Print("Test_stddev.pdf","pdf");};
}
}
}
Regards,
Soumyadip
Hi,
I have been able to solve the issue. I get to know that I have to use option s when taking the projection using ProfileX (ProfileX(histname,0,-1,"s"))
Here is the final code –
for(int id=0; id<ndef; id++){
for(int ij=0; ij<njet; ij++){
for(int ik=0; ik<nkappa; ik++){
DirData->cd();
sprintf(histname, "reco_jco_pt_d%i_j%i_k%i", id,ij,ik);
TProfile *TP_Data_Reco_JCO_PT = Data_Reco_JCO_PT[id][ij][ik]->ProfileX(histname,0,-1,"s");
TP_Data_Reco_JCO_PT->Rebin(3);
int ibins = TP_Data_Reco_JCO_PT->GetNbinsX();
double* xValues = new double[ibins];
double* binErrors = new double[ibins];
for(int ibin=1; ibin<=ibins; ++ibin){
xValues[ibin - 1] = TP_Data_Reco_JCO_PT->GetBinCenter(ibin);
binErrors[ibin - 1] =TP_Data_Reco_JCO_PT->GetBinError(ibin);
}
TGraph* graph = new TGraph(ibins, xValues, binErrors);
sprintf(histname,"reco_jco_pt_StdDev_d%i_j%i_k%i", id, ij, ik);
graph->SetNameTitle(histname,histname);
graph->Write();
}
}
}
But there is one more issue.
When I try to plot those, I want them to plot as points, not as connecting lines.
I have used AP to draw that, but when I am using it, I am getting one ik value in one canvas, but my intention is to plot all the ik values in one canvas. If I remove the AP option, it can plot all the ik values in one canvas, but it’s showing as connecting lines instead of points.
Could you kindly help me with that?
Here is my working code –
for (int id = 0; id < ndef; id++) {
for (int ij = 0; ij < njet; ij++) {
TCanvas* can4 = new TCanvas("can4", "Curve", 800, 900);
can4->SetLeftMargin(0.15);
TLegend* legend4_v1 = new TLegend(.7,0.6,0.85,0.85);
TLegend* legend4_v2 = new TLegend(.25,0.15,0.55,0.3);
legend4_v1->SetBorderSize(0);
legend4_v2->SetBorderSize(0);
for(int ik=0; ik<nkappa; ik++){
TGraph* Data_Reco_JCO_PT_CP = (TGraph*)Data_Reco_JCO_PT[id][ij][ik]->Clone();
//can4->cd();
Data_Reco_JCO_PT_CP->SetTitle("");
Data_Reco_JCO_PT_CP->GetXaxis()->CenterTitle();
Data_Reco_JCO_PT_CP->GetYaxis()->CenterTitle();
sprintf(xtitle, "Jet p_{T,%d} [GeV/c]",ij+1);
Data_Reco_JCO_PT_CP->GetXaxis()->SetTitle(xtitle);
sprintf(ytitle,"Q_{%s,%d}" ,obs_def[id], ij+1);
Data_Reco_JCO_PT_CP->GetYaxis()->SetTitle(ytitle);
Data_Reco_JCO_PT_CP->GetXaxis()->SetTitleSize(0.04); // Adjust the title size
Data_Reco_JCO_PT_CP->GetXaxis()->SetTitleOffset(1.0); // Adjust the title offset
Data_Reco_JCO_PT_CP->GetYaxis()->SetTitleSize(0.04); // Adjust the title size
Data_Reco_JCO_PT_CP->GetYaxis()->SetTitleOffset(1.5); // Adjust the title offset
Data_Reco_JCO_PT_CP->SetLineWidth(2);
Data_Reco_JCO_PT_CP->SetLineColor(color[ik]);
Data_Reco_JCO_PT_CP->SetMarkerColor(color[ik]);
Data_Reco_JCO_PT_CP->SetMarkerStyle(marker[ik]);
Data_Reco_JCO_PT_CP->SetMarkerSize(1.2);
Data_Reco_JCO_PT_CP->GetYaxis()->SetRangeUser(0, 2*Data_Reco_JCO_PT_CP->GetMaximum());
Data_Reco_JCO_PT_CP->GetXaxis()->SetRangeUser(0, Data_Reco_JCO_PT_CP->GetMaximum());
legend4_v1->AddEntry(Data_Reco_JCO_PT_CP, k_fact[ik], "lp");
if(ik==0){
Data_Reco_JCO_PT_CP->Draw("AP");
}else{
Data_Reco_JCO_PT_CP->Draw("same AP");
}
}
legend4_v1->Draw();
CMS_lumi( can4, 0, 0 );
can4->Update();
if(id==0 && ij==0){can4->Print("Data_StdDev.pdf(","pdf");}
else if(id ==2 && ij==1) {can4->Print("Data_StdDev.pdf)","pdf");}
else{can4->Print("Data_StdDev.pdf","pdf");};
}
}
“same” is not an option for drawing TGraphs, so remove that. Also remove the “A” from the “else” part (use it only for the 1st drawing).
if(ik==0){
Data_Reco_JCO_PT_CP->Draw("AP");
}else{
Data_Reco_JCO_PT_CP->Draw("P");
}
2 Likes
Thanks.
It worked fine.
This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.