analysis_INTT/plot/FitClusPhisize.py

0001 #! /usr/bin/env python
0002 from optparse import OptionParser
0003 import sys
0004 import os
0005 import datetime
0006 from array import *
0007 from ROOT import *
0008 import numpy
0009 import math
0010 import glob
0011 from plotUtil import *
0012
0013 titlesize = 0.05
0014 labelsize = 0.05
0015 pad1scale = 1.05
0016 pad2scale = 2.6
0017
0018 def fitfunction(fstr):
0019     if fstr == 'powerlaw':
0020         return '[0]*x^(-[1])', 'N#timesx^{-#lambda}'
0021     elif fstr == 'expo':
0022         return '[0]*exp(-x/[1])', 'N#timesexp(-x/#lambda)'
0023     elif fstr == 'pwexp':
0024         return '[0]*x^(-[1])*exp(-x/[2])', 'N#timesx^{-#alpha}#timesexp(-x/#lambda)'
0025     else:
0026         print('Invalid function. Returning powerlaw as default')
0027         return '[0]*x^{-[1]}', 'N#timesx^{-#lambda}'
0028
0029 def FillRandom(h, f, n):
0030     for i in range(n):
0031         h.Fill(f.GetRandom())
0032     return h
0033
0034 gROOT.SetBatch(True)
0035
0036 if __name__ == '__main__':
0037     gStyle.SetPalette(kSolar)
0038     os.makedirs('./ClusPhiSize/', exist_ok=True)
0039
0040     hM_ClusPhiSize_all_sim = GetHistogram('./hists/Sim_Ntuple_HIJING_ana443_20241030/Cluster/hists_merged.root', 'hM_ClusPhiSize_all')
0041     hM_ClusPhiSize_all_data = GetHistogram('./hists/Data_CombinedNtuple_Run54280_HotChannel_BCO/Cluster/hists_merged.root', 'hM_ClusPhiSize_all')
0042     # normalize the histograms
0043     hM_ClusPhiSize_all_sim.Scale(1/hM_ClusPhiSize_all_sim.Integral(-1,-1))
0044     hM_ClusPhiSize_all_data.Scale(1/hM_ClusPhiSize_all_data.Integral(-1,-1))
0045     maxbincontent = max(hM_ClusPhiSize_all_sim.GetMaximum(), hM_ClusPhiSize_all_data.GetMaximum())
0046     minbincontent = min(hM_ClusPhiSize_all_sim.GetMinimum(0), hM_ClusPhiSize_all_data.GetMinimum(0))
0047
0048     fitstr_data, fitleg_data = fitfunction('pwexp')
0049     f_data = TF1('f_data', fitstr_data, 1, 100)
0050     f_data.SetParameter(0, 2.3)
0051     f_data.SetParameter(1, 3.1)
0052     f_data.SetParameter(2, 15)
0053     # set parameter range for the fit
0054     # f_data.SetParLimits(0, 0, maxbincontent*10)
0055     # f_data.SetParLimits(1, 1E-6, 10)
0056     # f_data.SetParLimits(2, 10, 100)
0057     f_data.SetParNames('N', '#alpha', '#lambda')
0058     f_data.SetLineColorAlpha(kRed+1,0.5)
0059     f_data.Draw()
0060     # hM_ClusPhiSize_all_data.Fit('f_data', 'R')
0061     # print chi2/ndf
0062     # print('Chi2/ndf for data fit: {}'.format(f_data.GetChisquare()/f_data.GetNDF()))
0063
0064     fitstr_sim, fitleg_sim = fitfunction('pwexp')
0065     f_sim = TF1('f_sim', fitstr_sim, 1, 100)
0066     f_sim.SetParameter(0, 2.9)
0067     f_sim.SetParameter(1, 2.5)
0068     f_sim.SetParameter(2, 4.5)
0069     # set parameter range for the fit
0070     # f_sim.SetParLimits(0, 0, maxbincontent*10)
0071     # f_sim.SetParLimits(1, 1E-6, 50)
0072     # f_sim.SetParLimits(2, 1E-6, 10)
0073     f_sim.SetParNames('N', '#alpha', '#lambda')
0074     f_sim.SetLineColorAlpha(kTMutCyan,0.5)
0075     f_sim.Draw()
0076     # hM_ClusPhiSize_all_sim.Fit('f_sim', 'R')
0077     # print chi2/ndf
0078     # print('Chi2/ndf for sim fit: {}'.format(f_sim.GetChisquare()/f_sim.GetNDF()))
0079
0080     # produce histogram with random numbers from the fit function with various parameters
0081     ref_alpha = 3.0 # for fixed alpha
0082     ref_lambda = 4.5 # for fixed lambda
0083     nentires = 10000000
0084     nhists = 6
0085     l_alpha = [1.+0.5*i for i in range(nhists)]
0086     l_lambda = [4.+2*i for i in range(nhists)]
0087     l_hist_fixed_alpha = []
0088     l_hist_fixed_lambda = []
0089     for i in range(nhists):
0090         f_alpha = TF1('f_alpha', fitstr_sim, 1, 100)
0091         f_alpha.SetParameter(0, 1)
0092         f_alpha.SetParameter(1, l_alpha[i])
0093         f_alpha.SetParameter(2, ref_lambda)
0094         h_alpha = TH1F('h_alpha', 'h_alpha', 100, 1, 100)
0095         h_alpha = FillRandom(h_alpha, f_alpha, nentires)
0096         # normalize the histograms to 1
0097         h_alpha.Scale(1/h_alpha.Integral(-1,-1))
0098         l_hist_fixed_alpha.append(h_alpha)
0099
0100         f_lambda = TF1('f_lambda', fitstr_sim, 1, 100)
0101         f_lambda.SetParameter(0, 1)
0102         f_lambda.SetParameter(1, ref_alpha)
0103         f_lambda.SetParameter(2, l_lambda[i])
0104         h_lambda = TH1F('h_lambda', 'h_lambda', 100, 1, 100)
0105         h_lambda = FillRandom(h_lambda, f_lambda, nentires)
0106         # normalize the histograms to 1
0107         h_lambda.Scale(1/h_lambda.Integral(-1,-1))
0108         l_hist_fixed_lambda.append(h_lambda)
0109
0110     c = TCanvas('c','c',800,600)
0111     c.SetLogy()
0112     hM_ClusPhiSize_all_data.SetLineColor(kBlack)
0113     hM_ClusPhiSize_all_data.SetLineWidth(2)
0114     hM_ClusPhiSize_all_data.SetMarkerStyle(20)
0115     hM_ClusPhiSize_all_data.SetMarkerSize(0.8)
0116     hM_ClusPhiSize_all_data.GetYaxis().SetTitle('Normalized Counts')
0117     hM_ClusPhiSize_all_data.GetYaxis().SetTitleOffset(1.5)
0118     hM_ClusPhiSize_all_data.GetYaxis().SetRangeUser(minbincontent*0.5, maxbincontent*10)
0119     hM_ClusPhiSize_all_data.GetYaxis().SetTitleSize(titlesize)
0120     hM_ClusPhiSize_all_data.GetXaxis().SetTitle('Cluster #phi size')
0121     hM_ClusPhiSize_all_data.GetXaxis().SetRangeUser(1, 100)
0122     hM_ClusPhiSize_all_data.GetXaxis().SetTitleSize(titlesize)
0123     hM_ClusPhiSize_all_data.Draw('PE1')
0124     hM_ClusPhiSize_all_sim.SetLineColor(kTVibBlue)
0125     hM_ClusPhiSize_all_sim.SetFillColorAlpha(kTVibBlue, 0.1)
0126     hM_ClusPhiSize_all_sim.SetLineWidth(2)
0127     hM_ClusPhiSize_all_sim.Draw('hist same')
0128     f_data.Draw('same')
0129     f_sim.Draw('same')
0130     c.RedrawAxis()
0131     leg = TLegend(0.5, 0.6, 0.9, 0.9)
0132     leg.AddEntry(hM_ClusPhiSize_all_data, 'Data', 'PE')
0133     leg.AddEntry(hM_ClusPhiSize_all_sim, 'Simulation', 'lf')
0134     # leg.AddEntry(f_data, 'Function: ' + fitleg_data, 'l')
0135     # for i in range(f_data.GetNpar()):
0136     #     leg.AddEntry(0, f_data.GetParName(i) + ': {:.3e}'.format(f_data.GetParameter(i)) + ' #pm {:.3e}'.format(f_data.GetParError(i)), '')
0137     # leg.AddEntry(f_sim, 'Simulation Fit: ' + fitleg_sim, 'l')
0138     # for i in range(f_sim.GetNpar()):
0139     #     leg.AddEntry(0, f_sim.GetParName(i) + ': {:.3e}'.format(f_sim.GetParameter(i)) + ' #pm {:.3e}'.format(f_sim.GetParError(i)), '')
0140     leg.AddEntry(f_data, 'Function:' + fitleg_data, 'l')
0141     leg.AddEntry(0, f_data.GetParName(1) + '={:.1f}'.format(f_data.GetParameter(1)) + ', ' + f_data.GetParName(2) + '={:.1f}'.format(f_data.GetParameter(2)), '')
0142     leg.AddEntry(f_sim, 'Function:' + fitleg_sim, 'l')
0143     leg.AddEntry(0, f_sim.GetParName(1) + '={:.1f}'.format(f_sim.GetParameter(1)) + ', ' + f_sim.GetParName(2) + '={:.1f}'.format(f_sim.GetParameter(2)), '')
0144     leg.SetBorderSize(0)
0145     leg.SetFillStyle(0)
0146     leg.SetTextSize(0.035)
0147     leg.Draw()
0148     c.SaveAs('./ClusPhiSize/ClusPhiSize_fit.pdf')
0149     c.SaveAs('./ClusPhiSize/ClusPhiSize_fit.png')
0150
0151     # plot the histograms with fixed alpha
0152     c.Clear()
0153     c.SetLogy()
0154     hM_ClusPhiSize_all_data.GetYaxis().SetRangeUser(minbincontent*0.5, maxbincontent*100)
0155     hM_ClusPhiSize_all_data.Draw('PE1')
0156     hM_ClusPhiSize_all_sim.Draw('hist same')
0157     leg_2 = TLegend(0.3, 0.83, 0.6, 0.9)
0158     leg_2.SetNColumns(2)
0159     leg_2.AddEntry(hM_ClusPhiSize_all_data, 'Data', 'PE')
0160     leg_2.AddEntry(hM_ClusPhiSize_all_sim, 'Simulation', 'lf')
0161     leg_2.SetBorderSize(0)
0162     leg_2.SetFillStyle(0)
0163     leg_2.SetTextSize(0.04)
0164     leg_2.Draw()
0165     leg_alpha = TLegend(0.3, 0.66, 0.6, 0.82)
0166     leg_alpha.SetHeader('Toy model N#timesx^{-#alpha}#timesexp(-x/#lambda), fixed'+ ' #lambda={:.1f}'.format(ref_lambda) + ', varying #alpha')
0167     leg_alpha.SetNColumns(2)
0168     for i in range(nhists):
0169         # l_hist_fixed_alpha[i].SetLineColorAlpha(kRed+1+2*i, 0.5)
0170         l_hist_fixed_alpha[i].Draw('same hist PLC')
0171         leg_alpha.AddEntry(l_hist_fixed_alpha[i], '#alpha={:.1f}'.format(l_alpha[i]), 'l')
0172
0173     leg_alpha.SetBorderSize(0)
0174     leg_alpha.SetFillStyle(0)
0175     leg_alpha.SetTextSize(0.04)
0176     leg_alpha.Draw()
0177     c.SaveAs('./ClusPhiSize/ClusPhiSize_toys_fixed_alpha.pdf')
0178     c.SaveAs('./ClusPhiSize/ClusPhiSize_toys_fixed_alpha.png')
0179
0180     c.Clear()
0181     c.SetLogy()
0182     hM_ClusPhiSize_all_data.GetYaxis().SetRangeUser(minbincontent*0.5, maxbincontent*100)
0183     hM_ClusPhiSize_all_data.Draw('PE1')
0184     hM_ClusPhiSize_all_sim.Draw('hist same')
0185     leg_2.Draw()
0186     leg_lambda = TLegend(0.3, 0.66, 0.6, 0.82)
0187     leg_lambda.SetHeader('Toy model N#timesx^{-#alpha}#timesexp(-x/#lambda), fixed' + ' #alpha={:.1f}'.format(ref_alpha) + ', varying #lambda')
0188     leg_lambda.SetNColumns(2)
0189     for i in range(nhists):
0190         # l_hist_fixed_lambda[i].SetLineColorAlpha(kRed+1+2*i, 0.5)
0191         l_hist_fixed_lambda[i].Draw('same hist PLC')
0192         leg_lambda.AddEntry(l_hist_fixed_lambda[i], '#lambda={:.1f}'.format(l_lambda[i]), 'l')
0193
0194     leg_lambda.SetBorderSize(0)
0195     leg_lambda.SetFillStyle(0)
0196     leg_lambda.SetTextSize(0.04)
0197     leg_lambda.Draw()
0198     c.SaveAs('./ClusPhiSize/ClusPhiSize_toys_fixed_lambda.pdf')
0199     c.SaveAs('./ClusPhiSize/ClusPhiSize_toys_fixed_lambda.png')