Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add pT weights #330

Open
wants to merge 7 commits into
base: master
Choose a base branch
from
Open

Add pT weights #330

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
7 commits
Select commit Hold shift + click to select a range
e205393
Adds code for computing pT shape weights
wuctlby Dec 20, 2024
036e579
relative path
wuctlby Dec 20, 2024
90bc76b
change the way how to add reflection
wuctlby Jan 1, 2025
3f75e41
Upgrade the plot of v2VSfrac, and consider the case with only one pT …
wuctlby Jan 1, 2025
8211f7b
debug: suffix --> suffix_pdf
wuctlby Jan 6, 2025
df51030
copy the flow config
wuctlby Jan 6, 2025
3b703ae
debug the copuying of flow config
wuctlby Jan 6, 2025
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
518 changes: 518 additions & 0 deletions run3/flow/BDT/ComputePtWeights.py
View file
Open in desktop

Large diffs are not rendered by default.

24 changes: 16 additions & 8 deletions run3/flow/BDT/ComputeV2vsFDFrac.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -8,15 +8,15 @@
import yaml
import sys
from ROOT import TFile, TCanvas, TLegend, TLatex, TGraphErrors, TF1, TH1D, TVirtualFitter, Double_t
from ROOT import kBlack, kAzure, kRed
from ROOT import kBlack, kAzure, kCyan, kOrange
from ROOT import kFullCircle, kOpenCircle
sys.path.append('../../../')
sys.path.append('../')
from flow_analysis_utils import get_particle_info
from utils.StyleFormatter import SetGlobalStyle, SetObjectStyle, GetROOTColor

def set_frame_style(canv, Title, particleTit):
hFrame = canv.DrawFrame(0.0, 0.0001, 1, 0.35, f"{Title};Non-prompt {particleTit} fraction; #it{{v}}_{{2}}^{{#it{{obs}}}}")
hFrame = canv.DrawFrame(0.0, -0.2, 1, 0.35, f"{Title};Non-prompt {particleTit} fraction; #it{{v}}_{{2}}^{{#it{{obs}}}}")
hFrame.GetYaxis().SetDecimals()
hFrame.GetYaxis().SetNoExponent()
hFrame.GetXaxis().SetMoreLogLabels()
Expand Down Expand Up @@ -74,7 +74,6 @@ def v2_vs_frac(config, inputdir, outputdir, suffix):
hV2[-1].SetDirectory(0)
hFracFD[-1].SetDirectory(0)
hFracPrompt[-1].SetDirectory(0)
SetObjectStyle(gV2[-1], linecolor=1, linewidth=2, markerstyle=20, markersize=1, markercolor=1)

gFracVsV2, hV2VsFrac = [], [] # gFracVsV2 used for fitting, hV2VsFrac used for plotting
hV2VsPtFD = hV2[0].Clone("hV2VsPtFD")
Expand All @@ -88,6 +87,8 @@ def v2_vs_frac(config, inputdir, outputdir, suffix):
gFracVsV2.append(TGraphErrors(-1))
hV2VsFrac.append(TH1D(f"hV2VsFrac_{iPt}", "", 1000, 0.0, 1.0))
hV2VsFrac[-1].SetDirectory(0)
SetObjectStyle(hV2VsFrac[-1], markerstyle=kFullCircle, markersize=0)
SetObjectStyle(gFracVsV2[-1], linecolor=kAzure+4, linewidth=2, markerstyle=kFullCircle, markersize=1, markercolor=kAzure+4)

ptMin = hFracFD[0].GetBinLowEdge(iPt + 1)
ptMax = ptMin + hFracFD[0].GetBinWidth(iPt + 1)
Expand All @@ -107,15 +108,17 @@ def v2_vs_frac(config, inputdir, outputdir, suffix):
gFracVsV2[iPt].SetPoint(iSet, fracFD, v2)
gFracVsV2[iPt].SetPointError(iSet, fracFDUnc, v2Unc)

# gFracVsV2Fit = TGraphErrors(gFracVsV2[-1])
linFunc = TF1("linear", "pol1", 0, 1)
SetObjectStyle(linFunc, color=kOrange+1, linestyle=9, linewidth=2)
gFracVsV2[-1].Fit("linear", "", "", 0, 1)
chi2 = linFunc.GetChisquare()
ndf = linFunc.GetNDF()

# get the confidence intervals 0.683
fitter = TVirtualFitter.GetFitter()
fitter.GetConfidenceIntervals(hV2VsFrac[-1], 0.683)
hV2VsFrac[-1].SetLineColorAlpha(4, 0.15)
hV2VsFrac[-1].SetLineColorAlpha(kAzure+5, 0.15)

# get the v2 value at the FD fraction = 1, and it is not the last bin?
hV2VsPtFD.SetBinContent(iPt + 1,
Expand Down Expand Up @@ -151,8 +154,10 @@ def v2_vs_frac(config, inputdir, outputdir, suffix):
suffix_pdf = ')'
else:
suffix_pdf = ''
if nPtBins == 1:
suffix_pdf = ''

cFrac.append(TCanvas(f"cFrac_{ptStrings[iPt]}", "", 500, 500))
cFrac.append(TCanvas(f"cFrac_{ptStrings[iPt]}", "", 1200, 1200))
set_frame_style(cFrac[-1], ptStrings[iPt], particleTit)

t.SetTextSize(0.04)
Expand All @@ -161,12 +166,15 @@ def v2_vs_frac(config, inputdir, outputdir, suffix):
t.SetTextSize(0.035)
t.DrawLatex(0.250, 0.23, f'{chi2Strings[iPt]}')

gV2[iPt].Draw("pZ")
hV2VsFrac[iPt].DrawCopy("same")
hV2VsFrac[iPt].Draw("same pZ")
gFracVsV2[iPt].Draw("same pZ")
input("Press Enter to continue...")

gV2[iPt].Write()
gFracVsV2[iPt].Write()
hV2VsFrac[iPt].Write()

cFrac[-1].Update()

cFrac[iPt].SaveAs(f"{outputdir}/V2VsFrac/FracV2_{suffix}.pdf{suffix_pdf}")

PtTit = "#it{p}_{T} GeV/#it{c}"
Expand Down
33 changes: 18 additions & 15 deletions run3/flow/BDT/proj_thn_mc.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -13,13 +13,12 @@
from ROOT import gROOT
from alive_progress import alive_bar
from scipy.interpolate import InterpolatedUnivariateSpline
### please fill your path of DmeasonAnalysis
sys.path.append('../../..')
from utils.TaskFileLoader import LoadSparseFromTask



def proj_MC(config, cutsetConfig, outputdir, suffix):
def proj_MC(config, cutsetConfig, ptweights, ptweightsB, outputdir, suffix):

with open(config, 'r') as ymlCfgFile:
config = yaml.load(ymlCfgFile, yaml.FullLoader)
Expand All @@ -31,17 +30,15 @@ def proj_MC(config, cutsetConfig, outputdir, suffix):
particleName = config['Dmeson']
enableRef = config['enableRef']
enableSecPeak = config['enableSecPeak']
ptweights = [config['ptweightPath'], config['ptweightName']]
ptweightsB = [config['ptweightBPath'], config['ptweightBName']]
Bspeciesweights = config['Bspeciesweights'] if 'Bspeciesweights' in config else None

#TODO: safety checks for Dmeson reflecton and secondary peak

# check the arguments
if '' in ptweights:
if ptweights == []:
print('\033[91m WARNING: pt weights will not be provided! \033[0m')
ptweights = None
if '' in ptweightsB:
if ptweightsB == []:
print('\033[91m WARNING: B weights will not not be provided! \033[0m')
ptweightsB = None
if Bspeciesweights == None:
Expand Down Expand Up @@ -115,22 +112,22 @@ def proj_MC(config, cutsetConfig, outputdir, suffix):
sparseReco['RecoAll'].GetAxis(axisNum).SetRange(binMin, binMax)

if particleName == 'Dzero':
sparseReco['RecoPrompt'].GetAxis(6).SetRange(2, 2) # make sure it is prompt
sparseReco['RecoFD'].GetAxis(6).SetRange(3, 3) # make sure it is non-prompt
sparseReco['RecoPrompt'].GetAxis(8).SetRange(1, 2) # make sure it is signal
sparseReco['RecoFD'].GetAxis(8).SetRange(1, 2) # make sure it is signal
sparseReco['RecoPrompt'].GetAxis(8).SetRange(2, 2) # make sure it is prompt
sparseReco['RecoFD'].GetAxis(8).SetRange(3, 3) # make sure it is non-prompt
sparseReco['RecoPrompt'].GetAxis(6).SetRange(1, 2) # make sure it is signal
sparseReco['RecoFD'].GetAxis(6).SetRange(1, 2) # make sure it is signal
#TODO: add other particles
sparseReco['RecoPrompt'].GetAxis(axisNum).SetRange(binMin, binMax)
sparseReco['RecoFD'].GetAxis(axisNum).SetRange(binMin, binMax)

if enableRef:
sparseReco['RecoRefl'].GetAxis(8).SetRange(3, 4) # make sure it is reflection
sparseReco['RecoRefl'].GetAxis(6).SetRange(3, 4) # make sure it is reflection
sparseReco['RecoRefl'].GetAxis(axisNum).SetRange(binMin, binMax)
sparseReco['RecoReflPrompt'].GetAxis(8).SetRange(3, 4) # make sure it is reflection
sparseReco['RecoReflPrompt'].GetAxis(6).SetRange(2, 2) # make sure it is prompt reflection
sparseReco['RecoReflPrompt'].GetAxis(6).SetRange(3, 4) # make sure it is reflection
sparseReco['RecoReflPrompt'].GetAxis(8).SetRange(2, 2) # make sure it is prompt reflection
sparseReco['RecoReflPrompt'].GetAxis(axisNum).SetRange(binMin, binMax)
sparseReco['RecoReflFD'].GetAxis(8).SetRange(3, 4) # make sure it is reflection
sparseReco['RecoReflFD'].GetAxis(6).SetRange(3, 3) # make sure it is non-prompt reflection
sparseReco['RecoReflFD'].GetAxis(6).SetRange(3, 4) # make sure it is reflection
sparseReco['RecoReflFD'].GetAxis(8).SetRange(3, 3) # make sure it is non-prompt reflection
sparseReco['RecoReflFD'].GetAxis(axisNum).SetRange(binMin, binMax)

if enableSecPeak:
Expand Down Expand Up @@ -406,6 +403,10 @@ def proj_MC(config, cutsetConfig, outputdir, suffix):
default="config.yaml", help="flow configuration file")
parser.add_argument('cutsetConfig', metavar='text',
default='cutsetConfig.yaml', help='cutset configuration file')
parser.add_argument("--ptweights", "-w", metavar="text", nargs=2, required=False,
default=[], help="path to pt weights file and histogram name")
parser.add_argument("--ptweightsB", "-wb", metavar="text", nargs=2, required=False,
default=[], help="path to pt weightsB file and histogram name")
parser.add_argument("--outputdir", "-o", metavar="text",
default=".", help="output directory")
parser.add_argument("--suffix", "-s", metavar="text",
Expand All @@ -414,5 +415,7 @@ def proj_MC(config, cutsetConfig, outputdir, suffix):

proj_MC(config=args.config,
cutsetConfig=args.cutsetConfig,
ptweights=args.ptweights,
ptweightsB=args.ptweightsB,
outputdir=args.outputdir,
suffix=args.suffix)
54 changes: 47 additions & 7 deletions run3/flow/BDT/run_cutvar.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -20,6 +20,7 @@ def check_dir(dir):
return

def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffix, vn_method,
skip_calc_weights=False,
skip_make_yaml=False,
skip_cut_variation=False,
skip_proj_mc=False,
Expand All @@ -30,7 +31,7 @@ def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffi
skip_v2_vs_frac=False):

#___________________________________________________________________________________________________________________________
# Load the configuration file
# Load and copy the configuration file
with open(config_flow, 'r') as cfgFlow:
config = yaml.safe_load(cfgFlow)

Expand All @@ -51,6 +52,29 @@ def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffi
print(f"\033[32mNumber of cutsets: {nCutSets}\033[0m")

output_dir = f"{output}/cutvar_{suffix}"

os.system(f"mkdir -p {output_dir}")

# the pT weights histograms
PtWeightsDHistoName = 'hPtWeightsFONLLtimesTAMUDcent'
PtWeightsBHistoName = 'hPtWeightsFONLLtimesTAMUBcent'

# copy the configuration file
config_suffix = 1
while os.path.exists(f'{output_dir}/config_flow_{suffix}_{config_suffix}.yml'):
config_suffix = config_suffix + 1
os.system(f'cp {config_flow} {output_dir}/config_flow_{suffix}_{config_suffix}.yml')

#___________________________________________________________________________________________________________________________
# calculate the pT weights
if not skip_calc_weights:
check_dir(f"{output_dir}/ptweights")
CalcWeiPath = "./ComputePtWeights.py"

print(f"\033[32mpython3 {CalcWeiPath} {config_flow} -o {output_dir} -s {suffix}\033[0m")
os.system(f"python3 {CalcWeiPath} {config_flow} -o {output_dir} -s {suffix}")
else:
print("\033[33mWARNING: Calculation of weights will not be performed\033[0m")

#___________________________________________________________________________________________________________________________
# make yaml file
Expand Down Expand Up @@ -80,11 +104,25 @@ def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffi
if not skip_proj_mc:
check_dir(f"{output_dir}/proj_mc")
ProjMcPath = "./proj_thn_mc.py"

if not os.path.exists(f'{output_dir}/ptweights/pTweight_{suffix}.root'):
for i in range(nCutSets):
iCutSets = f"{i:02d}"
print(f"\033[32mpython3 {ProjMcPath} {config_flow} {output_dir}/config/cutset_{suffix}_{iCutSets}.yml -o {output_dir} -s {suffix}_{iCutSets}\033[0m")
os.system(f"python3 {ProjMcPath} {config_flow} {output_dir}/config/cutset_{suffix}_{iCutSets}.yml -o {output_dir} -s {suffix}_{iCutSets}")
else:
for i in range(nCutSets):
iCutSets = f"{i:02d}"
print(
f"\033[32mpython3 {ProjMcPath} {config_flow} {output_dir}/config/cutset_{suffix}_{iCutSets}.yml "
f"-w {output_dir}/ptweights/pTweight_{suffix}.root hPtWeightsFONLLtimesTAMUDcent "
f"-wb {output_dir}/ptweights/pTweight_{suffix}.root hPtWeightsFONLLtimesTAMUBcent "
f"-o {output_dir} -s {suffix}_{iCutSets} \033[0m"
)
os.system(f"python3 {ProjMcPath} {config_flow} {output_dir}/config/cutset_{suffix}_{iCutSets}.yml \
-w {output_dir}/ptweights/pTweight_{suffix}.root hPtWeightsFONLLtimesTAMUDcent \
-wb {output_dir}/ptweights/pTweight_{suffix}.root hPtWeightsFONLLtimesTAMUBcent -o {output_dir} -s {suffix}_{iCutSets}")

for i in range(nCutSets):
iCutSets = f"{i:02d}"
print(f"\033[32mpython3 {ProjMcPath} {config_flow} {output_dir}/config/cutset_{suffix}_{iCutSets}.yml -o {output_dir} -s {suffix}_{iCutSets}\033[0m")
os.system(f"python3 {ProjMcPath} {config_flow} {output_dir}/config/cutset_{suffix}_{iCutSets}.yml -o {output_dir} -s {suffix}_{iCutSets}")
else:
print("\033[33mWARNING: Projection for MC will not be performed\033[0m")

Expand All @@ -110,9 +148,9 @@ def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffi

for i in range(nCutSets):
iCutSets = f"{i:02d}"
print(f"\033[32mpython3 {SimFitPath} {config_flow} {cent} {output_dir}/proj/proj_{suffix}.root -o {output_dir}/ry -s {suffix} -vn {vn_method}\033[0m")
print(f"\033[32mpython3 {SimFitPath} {config_flow} {cent} {output_dir}/proj/proj_{suffix}.root -o {output_dir}/ry -s _{suffix}_{iCutSets} -vn {vn_method}\033[0m")
print(f"\033[32mProcessing cutset {iCutSets}\033[0m")
os.system(f"python3 {SimFitPath} {config_flow} {cent} {output_dir}/proj/proj_{suffix}_{iCutSets}.root -o {output_dir}/ry -s {suffix}_{iCutSets} -vn {vn_method} --batch")
os.system(f"python3 {SimFitPath} {config_flow} {cent} {output_dir}/proj/proj_{suffix}_{iCutSets}.root -o {output_dir}/ry -s _{suffix}_{iCutSets} -vn {vn_method} --batch")
else:
print("\033[33mWARNING: vn extraction will not be performed\033[0m")

Expand Down Expand Up @@ -161,6 +199,7 @@ def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffi
parser.add_argument("--outputdir", "-o", metavar="text", default=".", help="output directory")
parser.add_argument("--suffix", "-s", metavar="text", default="", help="suffix for output files")
parser.add_argument("--vn_method", "-vn", metavar="text", default="sp", help="vn technique (sp, ep, deltaphi)")
parser.add_argument("--skip_calc_weights", "-scw", action="store_true", help="skip calculation of weights")
parser.add_argument("--skip_make_yaml", "-smy", action="store_true", help="skip make yaml")
parser.add_argument("--skip_cut_variation", "-scv", action="store_true", help="skip cut variation")
parser.add_argument("--skip_proj_mc", "-spm", action="store_true", help="skip projection for MC")
Expand All @@ -172,6 +211,7 @@ def run_full_cut_variation(config_flow, anres_dir, cent, res_file, output, suffi
args = parser.parse_args()

run_full_cut_variation(args.flow_config, args.anres_dir, args.centrality, args.resolution, args.outputdir, args.suffix, args.vn_method,
args.skip_calc_weights,
args.skip_make_yaml,
args.skip_cut_variation,
args.skip_proj_mc,
Expand Down
35 changes: 22 additions & 13 deletions run3/flow/BDT/run_cutvar.sh
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,6 +9,7 @@
#- outputdir (str): output directory
#- suffix (str): suffix for output files
#- vn_method (str): vn technique (sp, ep, deltaphi)
#- skip_calc_weights (bool): skip calculation of weights
#- skip_make_yaml (bool): skip make yaml
#- skip_cut_variation (bool): skip cut variation
#- skip_proj_mc (bool): skip projection for MC
Expand All @@ -19,22 +20,29 @@
#- skip_v2_vs_frac (bool): skip v2 vs FD fraction
#----------

export config_flow="path/to/config_flow.yml"
export anres_dir="/media/wuct/wulby/ALICE/AnRes/D0_flow/pass4/ML/Results/306077/AnalysisResults_data.root"
export output_dir="/home/wuct/ALICE/local/DmesonAnalysis/run3/flow/BDT/results/test"
export config_flow="/home/wuct/ALICE/local/DmesonAnalysis/run3/flow/config_flow_d0_test_ml.yml"
export anres_dir="/media/wuct/wulby/ALICE/AnRes/D0_flow/pass4/ML/Results/CombPID_3SigCut_AND/AnalysisResults_data_medium.root"
export output_dir="/home/wuct/ALICE/local/DmesonAnalysis/run3/flow/BDT/results/medium"
export cent="k3050"
export vn_method="sp"
export res_file="/home/wuct/ALICE/local/DmesonAnalysis/run3/flow/Results/2060/k3050/small/sp/resolution/resosp3050s_291131_inte_gain_Reso.root"
export suffix="small_pt2_3"
export res_file="/media/wuct/wulby/ALICE/AnRes/resolution/output_reso/resospk3050_inte.root"
export suffix="pt2_3_full"

export smy=False # True or False
export scv=False # True or False
export spm=False # True or False
export seff=False # True or False
export svn=False # True or False
export sfcv=False # True or False
export sddf=False # True or False
export sv2vf=False # True or False
export spw=False # True or False (skip calculation of weights)
export smy=False # True or False (skip make yaml)
export scv=False # True or False (skip cut variation)
export spm=False # True or False (skip projection for MC)
export seff=False # True or False (skip efficiency)
export svn=False # True or False (skip vn extraction)
export sfcv=False # True or False (skip fraction by cut variation)
export sddf=False # True or False (skip fraction by data-driven method)
export sv2vf=False # True or False (skip v2 vs fraction)

if [ $spw = False ]; then
export skip_calc_weights=""
else
export skip_calc_weights="--skip_calc_weights"
fi

if [ $smy = False ]; then
export skip_make_yaml=""
Expand Down Expand Up @@ -85,6 +93,7 @@ else
fi

python3 run_cutvar.py $config_flow $anres_dir -c $cent -r $res_file -o $output_dir -s $suffix -vn $vn_method \
$skip_calc_weights \
$skip_make_yaml \
$skip_cut_variation \
$skip_proj_mc \
Expand Down
Loading