diff --git a/machine_learning_hep/analysis/analyzerdhadrons.py b/machine_learning_hep/analysis/analyzerdhadrons.py index 9b77f0819a..f541dfe0db 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons.py +++ b/machine_learning_hep/analysis/analyzerdhadrons.py @@ -29,18 +29,9 @@ TF1, TH1, TH1F, - TH2F, - TArrow, TCanvas, - TDirectory, TFile, TLegend, - TLine, - TPad, - TPaveLabel, - TPaveText, - TText, - gInterpreter, gPad, gROOT, gStyle, @@ -131,7 +122,7 @@ def __init__(self, datap, case, typean, period): self.fit_func_bkg = {} self.fit_range = {} - self.path_fig = Path(f'{os.path.expandvars(self.d_resultsallpdata)}/fig') + self.path_fig = Path(f"{os.path.expandvars(self.d_resultsallpdata)}/fig") for folder in ["qa", "fit", "roofit", "sideband", "signalextr", "fd", "uf"]: (self.path_fig / folder).mkdir(parents=True, exist_ok=True) @@ -147,6 +138,7 @@ def __init__(self, datap, case, typean, period): self.p_anahpt = datap["analysis"]["anahptspectrum"] self.p_fd_method = datap["analysis"]["fd_method"] + self.p_crosssec_prompt = datap["analysis"]["crosssec_prompt"] self.p_cctype = datap["analysis"]["cctype"] self.p_inputfonllpred = datap["analysis"]["inputfonllpred"] self.p_triggereff = datap["analysis"][self.typean].get("triggereff", [1]) @@ -160,7 +152,7 @@ def __init__(self, datap, case, typean, period): # region helpers def _save_canvas(self, canvas, filename): # folder = self.d_resultsallpmc if mcordata == 'mc' else self.d_resultsallpdata - canvas.SaveAs(f'{self.path_fig}/{filename}') + canvas.SaveAs(f"{self.path_fig}/{filename}") def _save_hist(self, hist, filename, option=""): if not hist: @@ -333,12 +325,10 @@ def fit(self): if self.cfg("mass_roofit"): for entry in self.cfg("mass_roofit", []): - if lvl := entry.get("level"): - if lvl != level: - continue - if ptspec := entry.get("ptrange"): - if ptspec[0] > ptrange[0] or ptspec[1] < ptrange[1]: - continue + if (lvl := entry.get("level")) and lvl != level: + continue + if (ptspec := entry.get("ptrange")) and (ptspec[0] > ptrange[0] or ptspec[1] < ptrange[1]): + continue fitcfg = entry break self.logger.debug("Using fit config for %i: %s", ipt, fitcfg) @@ -555,6 +545,7 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b nameyield, selnorm, self.p_sigmamb, + self.p_crosssec_prompt, output_prompt, fileoutcross, ) diff --git a/machine_learning_hep/analysis/analyzerdhadrons_mult.py b/machine_learning_hep/analysis/analyzerdhadrons_mult.py index 3048078932..cd62e163eb 100644 --- a/machine_learning_hep/analysis/analyzerdhadrons_mult.py +++ b/machine_learning_hep/analysis/analyzerdhadrons_mult.py @@ -28,22 +28,10 @@ from ROOT import ( TF1, TH1, - TH1D, TH1F, - TH2F, - TArrow, TCanvas, - TDirectory, TFile, TLegend, - TLine, - TPad, - TPaveLabel, - TPaveText, - TStyle, - TText, - gInterpreter, - gPad, gROOT, gStyle, kBlue, @@ -62,7 +50,6 @@ ) from machine_learning_hep.hf_pt_spectrum import hf_pt_spectrum from machine_learning_hep.logger import get_logger -from machine_learning_hep.root import save_root_object from machine_learning_hep.utils.hist import get_dim, project_hist @@ -144,7 +131,8 @@ def __init__(self, datap, case, typean, period): self.fit_func_bkg = {} self.fit_range = {} - self.path_fig = Path(f"fig/{self.case}/{self.typean}") + # self.path_fig = Path(f"fig/{self.case}/{self.typean}") + self.path_fig = Path(f"{os.path.expandvars(self.d_resultsallpdata)}/fig") for folder in ["qa", "fit", "roofit", "sideband", "signalextr", "fd", "uf"]: (self.path_fig / folder).mkdir(parents=True, exist_ok=True) @@ -156,6 +144,7 @@ def __init__(self, datap, case, typean, period): self.p_anahpt = datap["analysis"]["anahptspectrum"] self.p_fd_method = datap["analysis"]["fd_method"] + self.p_crosssec_prompt = datap["analysis"]["crosssec_prompt"] self.p_cctype = datap["analysis"]["cctype"] self.p_sigmamb = datap["analysis"]["sigmamb"] self.p_inputfonllpred = datap["analysis"]["inputfonllpred"] @@ -190,7 +179,8 @@ def __init__(self, datap, case, typean, period): # region helpers def _save_canvas(self, canvas, filename): # folder = self.d_resultsallpmc if mcordata == 'mc' else self.d_resultsallpdata - canvas.SaveAs(f"fig/{self.case}/{self.typean}/{filename}") + # canvas.SaveAs(f"fig/{self.case}/{self.typean}/{filename}") + canvas.SaveAs(f"{self.path_fig}/{filename}") def _save_hist(self, hist, filename, option=""): if not hist: @@ -378,12 +368,12 @@ def fit(self): if self.cfg("mass_roofit"): for entry in self.cfg("mass_roofit", []): - if lvl := entry.get("level"): - if lvl != level: - continue - if ptspec := entry.get("ptrange"): - if ptspec[0] > ptrange[0] or ptspec[1] < ptrange[1]: - continue + if (lvl := entry.get("level")) and lvl != level: + continue + if (ptspec := entry.get("ptrange")) and ( + ptspec[0] > ptrange[0] or ptspec[1] < ptrange[1] + ): + continue fitcfg = entry break self.logger.debug("Using fit config for %i: %s", ipt, fitcfg) @@ -491,14 +481,8 @@ def efficiency(self): legsl.SetTextSize(0.035) for imult in range(self.p_nbin2): - stringbin2 = "_{}_{:.2f}_{:.2f}".format( - self.v_var2_binning, self.lvar2_binmin[imult], self.lvar2_binmax[imult] - ) - legeffstring = "{:.1f} #leq {} < {:.1f}".format( - self.lvar2_binmin[imult], - self.p_latexbin2var, - self.lvar2_binmax[imult], - ) + stringbin2 = f"_{self.v_var2_binning}_{self.lvar2_binmin[imult]:.2f}_{self.lvar2_binmax[imult]:.2f}" + legeffstring = f"{self.lvar2_binmin[imult]:.1f} #leq {self.p_latexbin2var} < {self.lvar2_binmax[imult]:.1f}" if self.signal_loss: h_gen_pr_sl = lfileeff.Get("h_signal_loss_gen_pr" + stringbin2) @@ -571,13 +555,9 @@ def efficiency(self): legeffFD.SetTextSize(0.035) for imult in range(self.p_nbin2): - stringbin2 = "_{}_{:.2f}_{:.2f}".format( - self.v_var2_binning, self.lvar2_binmin[imult], self.lvar2_binmax[imult] - ) - legeffFDstring = "{:.1f} #leq {} < {:.1f}".format( - self.lvar2_binmin[imult], - self.p_latexbin2var, - self.lvar2_binmax[imult], + stringbin2 = f"_{self.v_var2_binning}_{self.lvar2_binmin[imult]:.2f}_{self.lvar2_binmax[imult]:.2f}" + legeffFDstring = ( + f"{self.lvar2_binmin[imult]:.1f} #leq {self.p_latexbin2var} < {self.lvar2_binmax[imult]:.1f}" ) if self.signal_loss: @@ -674,10 +654,8 @@ def plotter(self): hcross.GetYaxis().SetTitle(f"d#sigma/d#it{{p}}_{{T}} ({self.p_latexnhadron}) {self.typean}") hcross.SetName("hcross%d" % imult) hcross.GetYaxis().SetRangeUser(1e1, 1e10) - legvsvar1endstring = "{:.1f} < {} < {:.1f}".format( - self.lvar2_binmin[imult], - self.p_latexbin2var, - self.lvar2_binmax[imult], + legvsvar1endstring = ( + f"{self.lvar2_binmin[imult]:.1f} < {self.p_latexbin2var} < {self.lvar2_binmax[imult]:.1f}" ) legvsvar1.AddEntry(hcross, legvsvar1endstring, "LEP") hcross.Draw("same") @@ -781,14 +759,12 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b # pylint: disable=logging-not-lazy self.logger.warning("Number of events %d for mult bin %d" % (norm, imult)) + if self.p_nevents is not None: + norm = self.p_nevents if self.p_fprompt_from_mb: if imult == 0: - fileoutcrossmb = "{}/finalcross{}{}mult0.root".format( - self.d_resultsallpdata, self.case, self.typean - ) + fileoutcrossmb = f"{self.d_resultsallpdata}/finalcross{self.case}{self.typean}mult0.root" output_prompt = [] - if self.p_nevents is not None: - norm = self.p_nevents self.logger.warning("Corrected Number of events %d for mult bin %d" % (norm, imult)) hf_pt_spectrum( self.p_anahpt, @@ -803,6 +779,7 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b nameyield, norm, self.p_sigmamb, + self.p_crosssec_prompt, output_prompt, fileoutcrossmb, ) @@ -824,10 +801,12 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b nameyield, norm, self.p_sigmamb, + self.p_crosssec_prompt, output_prompt, fileoutcrossmult, ) else: + output_prompt = [] hf_pt_spectrum( self.p_anahpt, self.p_br, @@ -841,6 +820,7 @@ def makenormyields(self): # pylint: disable=import-outside-toplevel, too-many-b nameyield, norm, self.p_sigmamb, + self.p_crosssec_prompt, output_prompt, fileoutcrossmult, ) @@ -884,17 +864,13 @@ def plotternormyields(self): hcross.GetYaxis().SetTitleOffset(1.3) hcross.GetYaxis().SetTitle(f"Corrected yield/events ({self.p_latexnhadron}) {self.typean}") hcross.GetYaxis().SetRangeUser(1e-10, 1) - legvsvar1endstring = "{:.1f} #leq {} < {:.1f}".format( - self.lvar2_binmin[imult], - self.p_latexbin2var, - self.lvar2_binmax[imult], + legvsvar1endstring = ( + f"{self.lvar2_binmin[imult]:.1f} #leq {self.p_latexbin2var} < {self.lvar2_binmax[imult]:.1f}" ) legvsvar1.AddEntry(hcross, legvsvar1endstring, "LEP") hcross.Draw("same") legvsvar1.Draw() cCrossvsvar1.SaveAs( - "{}/CorrectedYieldsNorm{}{}Vs{}.eps".format( - self.d_resultsallpdata, self.case, self.typean, self.v_var_binning - ) + f"{self.d_resultsallpdata}/CorrectedYieldsNorm{self.case}{self.typean}Vs{self.v_var_binning}.eps" ) fileoutcrosstot.Close() diff --git a/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_newformat.yml b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi.yml similarity index 99% rename from machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_newformat.yml rename to machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi.yml index bec1d33b89..3cc0fc1b3a 100644 --- a/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_newformat.yml +++ b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi.yml @@ -350,10 +350,11 @@ LcpKpi: analysis: anahptspectrum: "LctopKpi" #D0Kpi, DplusKpipi, DstarD0pi, DsKKpi, LctopKpi, LcK0Sp - fd_method: "Nb" # fc, Nb, ext + fd_method: "Nb" # fc, Nb, ext, dd, dd_N + crosssec_prompt: true # true for prompt, false for non-prompt cctype: "pp" sigmamb: 59.4e+9 # 50.87e+9 pp5TeV, 57.8e+9 pp13TeV, 59.4e+9 pp Run3, pb - inputfonllpred: data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root + inputfonllpred: data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root #data/fonll/CutVarLc_pp13TeV_LHC23_pass4_wide.root dir_general_plots: analysis_plots Run3analysis: diff --git a/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_CrystalBall.yml b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_CrystalBall.yml new file mode 100644 index 0000000000..67be528008 --- /dev/null +++ b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_CrystalBall.yml @@ -0,0 +1,612 @@ +--- +############################################################################# +## © Copyright CERN 2023. All rights not expressly granted are reserved. ## +## Author: Gian.Michele.Innocenti@cern.ch ## +## This program is free software: you can redistribute it and/or modify it ## +## under the terms of the GNU General Public License as published by the ## +## Free Software Foundation, either version 3 of the License, or (at your ## +## option) any later version. This program is distributed in the hope that ## +## it will be useful, but WITHOUT ANY WARRANTY; without even the implied ## +## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. ## +## See the GNU General Public License for more details. ## +## You should have received a copy of the GNU General Public License ## +## along with this program. if not, see . ## +############################################################################# + +LcpKpi: + doml: true + mass: 2.286 + sel_reco_unp: "fPt>0" + sel_gen_unp: "fPt>0" + sel_cen_unp: null + sel_good_evt_unp: null + #sel_reco_skim: ["fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015"] + sel_reco_skim: [null, null, null, null, null, null, null, null, null, null] + sel_gen_skim: [null, null, null, null, null, null, null, null, null, null] + sel_skim_binmin: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12] #list of nbins + sel_skim_binmax: [2, 3, 4, 5, 6, 7, 8, 10, 12, 24] #list of nbins + apply_yptacccut: false + var_binning: fPt + dofullevtmerge: false + var_cand: fCandidateSelFlag + var_swap: fIsCandidateSwapped + bitmap_sel: + var_name: fFlagMcMatchRec + var_name_gen: fFlagMcMatchGen + var_name_origgen: fOriginMcGen + var_name_origrec: fOriginMcRec + var_isstd: isstd + var_ismcsignal: ismcsignal + var_ismcprompt: ismcprompt + var_ismcfd: ismcfd + var_ismcbkg: ismcbkg + var_ismcrefl: ismcref + isstd: [[1], []] + ismcsignal: [[1], []] + ismcprompt: [[0], []] + ismcfd: [[1], []] + ismcbkg: [[], [1]] + ismcrefl: [[1], [1]] + + dfs: + read: + evtorig: + level: all + index: fIndexHFLCCOLLBASES + trees: + O2hflccollbase: [fNumContrib, fMultZeqNTracksPV, fCentFT0M] + extra: + fIsEventReject: 0 + + evtoriggen: + level: gen + index: fIndexHFLCMCCOLLBASES + trees: + O2hflcmccollbase: [fPosX, fPosY, fPosZ, fCentFT0M] + O2hflcmcrcollid: [fIndexArrayHFLCCOLLBASES] + rename: {old: fCentFT0M, new: fCentFT0Mmc} + + reco: + level: all + index: fIndexHfLcBases + trees: + O2hflcbase: [fIndexHFLCCOLLBASES, fPt, fY, fEta, fPhi, fM] + O2hflcpar: [fNProngsContributorsPV, fCpa, fCpaXY, fChi2PCA, fDecayLength, fDecayLengthXY, fPtProng0, fPtProng1, fPtProng2, fImpactParameter0, fImpactParameter1, fImpactParameter2, fNSigTpcPi0, fNSigTpcPr0, fNSigTpcKa1, fNSigTpcPi2, fNSigTpcPr2, fNSigTofPi0, fNSigTofPr0, fNSigTofKa1, fNSigTofPi2, fNSigTofPr2, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + O2hflcmc: + level: mc + vars: [fFlagMcMatchRec, fOriginMcRec, fIsCandidateSwapped] + O2hflcsel: + level: mc + vars: [fCandidateSelFlag] + filter: "fPt > 0." + #extra: + #fY: log((sqrt(2.28646**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(2.28646**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated + tags: + isstd: {var: fFlagMcMatchRec, req: [[1], []], level: mc} + ismcsignal: {var: fFlagMcMatchRec, req: [[1], []], abs: true, level: mc} + ismcbkg: {var: fFlagMcMatchRec, req: [[], [1]], abs: true, level: mc} + ismcprompt: {var: fOriginMcRec, req: [[0], []], level: mc} + ismcfd: {var: fOriginMcRec, req: [[1], []], level: mc} + swap: {cand: fCandidateSelFlag, var_swap: fIsCandidateSwapped, vars: [ismcsignal, ismcprompt, ismcfd], level: mc} + + gen: + level: mc + trees: + O2hflcpbase: [fIndexHFLCMCCOLLBASES, fPt, fY, fEta, fPhi, fFlagMcMatchGen, fOriginMcGen] + tags: + isstd: {var: fFlagMcMatchGen, req: [[1], []], level: mc} + ismcsignal: {var: fFlagMcMatchGen, req: [[1], []], abs: true, level: mc} + ismcbkg: {var: fFlagMcMatchGen, req: [[], [1]], abs: true, level: mc} + ismcprompt: {var: fOriginMcGen, req: [[0], []], level: mc} + ismcfd: {var: fOriginMcGen, req: [[1], []], level: mc} + #extra: + #fY: log((sqrt(2.28646**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(2.28646**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated + + merge: + - {base: reco, ref: evtorig, extra: {fMultZeqNTracksPV_sub: fMultZeqNTracksPV - fNProngsContributorsPV}} + - {base: gen, ref: evtoriggen} + # workaround for yamlfmt issue #110 + #- {base: gen, ref: evtorig, left_on: fIndexArrayHFLCCOLLBASES, out: genrec} + + write: + evtorig: + level: all + file: AnalysisResultsEvtOrig.parquet + evt: + level: all + source: evtorig + file: AnalysisResultsEvt.parquet + evtmc: + level: mc + source: evtoriggen + file: AnalysisResultsEvtGen.parquet + reco: + level: all + file: AnalysisResultsReco.parquet + gen: + level: mc + file: AnalysisResultsGen.parquet + + variables: + var_all: [fIndexCollisions, fFlagMcMatchRec, fCandidateSelFlag, fOriginMcRec, fIsCandidateSwapped, fNProngsContributorsPV, fY, fEta, fPt, fCpa, fCpaXY, fM, fChi2PCA, fDecayLength, fDecayLengthXY, fPtProng0, fPtProng1, fPtProng2, fImpactParameter0, fImpactParameter1, fImpactParameter2, fNSigTpcPi0, fNSigTpcPr0, fNSigTpcKa1, fNSigTpcPi2, fNSigTpcPr2, fNSigTofPi0, fNSigTofPr0, fNSigTofKa1, fNSigTofPi2, fNSigTofPr2, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + var_jet: [fJetPt, fJetEta, fJetPhi] + var_jetsub: [fZg, fRg, fNsd] + var_jet_match: [df, fIndexHfCand2Prong] + var_jetsub_match: [df, fIndexLcChargedJets] + var_evt: + data: [fIndexCollisions, fPosX, fPosY, fPosZ, fNumContrib, fMultZeqNTracksPV, fCentFT0M] + mc: [fIndexCollisions, fPosX, fPosY, fPosZ, fMultZeqNTracksPV, fCentFT0M] + #mc: [fIndexCollisions, fPosX, fPosY, fPosZ, fNumContrib, fMultZeqNTracksPV, fCentFT0A, fCentFT0C, fCentFT0M, fCentFV0A] + var_gen: [fIndexMcCollisions, fPosX, fPosY, fPosZ, fPt, fY, fFlagMcMatchGen, fOriginMcGen] + var_training: [[fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2], [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2]] + var_selected: [fM, fY, fEta, fPt, fCpa, fCpaXY, fChi2PCA, fDecayLength, fDecayLengthXY, fPtProng0, fPtProng1, fPtProng2, fImpactParameter0, fImpactParameter1, fImpactParameter2, fNSigTpcPi0, fNSigTpcPr0, fNSigTpcKa1, fNSigTpcPi2, fNSigTpcPr2, fNSigTofPi0, fNSigTofPr0, fNSigTofKa1, fNSigTofPi2, fNSigTofPr2, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + var_boundaries: [fDecayLength, fPt] + var_correlation: + - [fDecayLength, fChi2PCA, fCpa] + - [fPt, fPt, fPt] + var_class: class + var_inv_mass: fM + var_y: fY + var_cuts: + - [fPtProng0, lt, null] + - [fPtProng1, lt, null] + - [fPtProng2, lt, null] + - [fCpa, lt, null] + - [fDecayLength, lt, null] + - [fChi2PCA, lt, null] + + plot_options: + prob_cut_scan: + pt_prong0: + xlim: + - 0 + - 8 + pt_prong1: + xlim: + - 0 + - 8 + pt_prong2: + xlim: + - 0 + - 8 + fDecayLength: + xlim: + - 0 + - 0.08 + fChi2PCA: + xlim: + - 0 + - 20. + fNSigTofPr0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(p)0" + fNSigTofPi0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(\\pi)0" + fNSigTofKa1: + xlim: [-10, 10] + xlabel: "n\\sigma_\\mathrm{TOF}(K)1" + fNSigTofPr2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(p)2" + fNSigTofPi2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(\\pi)2" + fNSigTpcPr0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(p)0" + fNSigTpcPi0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(\\pi)0" + fNSigTpcKa1: + xlim: [-10, 10] + xlabel: "n\\sigma_\\mathrm{TPC}(K)1" + fNSigTpcPr2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(p)2" + fNSigTpcPi2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(\\pi)2" + eff_cut_scan: + pt_prong0: + xlim: + - 0 + - 8 + pt_prong1: + xlim: + - 0 + - 8 + pt_prong2: + xlim: + - 0 + - 8 + fDecayLength: + xlim: + - 0 + - 0.08 + fChi2PCA: + xlim: + - 0 + - 20. + files_names: + namefile_unmerged_tree: AO2D.root + namefile_reco: AnalysisResultsReco.parquet + namefile_evt: AnalysisResultsEvt.parquet + namefile_collcnt: AnalysisResultsCollCnt.parquet + namefile_bccnt: AnalysisResultsBcCnt.parquet + namefile_evtvalroot: AnalysisResultsROOTEvtVal.root + namefile_evtorig: AnalysisResultsEvtOrig.parquet + namefile_gen: AnalysisResultsGen.parquet + namefile_reco_applieddata: AnalysisResultsRecoAppliedData.parquet + namefile_reco_appliedmc: AnalysisResultsRecoAppliedMC.parquet + namefile_mcweights: mcweights.root + treeoutput: "Lctree" + histofilename: "masshisto.root" + efffilename: "effhisto.root" + respfilename: "resphisto.root" + crossfilename: "cross_section_tot.root" + resultfilename: "finalcrossLcpKpiRun3analysis.root" + + multi: + data: + nprocessesparallel: 80 + maxfiles: [-1] #list of periods + chunksizeunp: [100] #list of periods + chunksizeskim: [100] #list of periods + fracmerge: [0.05] #list of periods + seedmerge: [12] #list of periods + period: [LHC23] #list of periods + select_period: [1] + prefix_dir: /data2/MLhep/ + #unmerged_tree_dir: [real/train_339425/alice/cern.ch/user/a/alihyperloop/jobs/0130] #list of periods small + unmerged_tree_dir: [real/train_342351/alice/cern.ch/user/a/alihyperloop/jobs/0133] #list of periods + pkl: [LHC23pp_mb/period_LHC23/pkldata] #list of periods + pkl_skimmed: [LHC23pp_mb/period_LHC23/pklskdata] #list of periods + pkl_skimmed_merge_for_ml: [LHC23pp_mb/period_LHC23/pklskmldata] #list of periods + pkl_skimmed_merge_for_ml_all: LHC23pp_mb/mltotdata + pkl_evtcounter_all: LHC23pp_mb/evttotdata + #select_jobs: [[hy_189959], [hy_189000]] + mcreweights: [../Analyses] + mc: + nprocessesparallel: 80 + maxfiles: [-1] #list of periods + chunksizeunp: [100] #list of periods + chunksizeskim: [100] #list of periods + fracmerge: [1.0] #list of periods + seedmerge: [1] #list of periods + period: [LHC24h1] #list of periods + select_period: [1] + prefix_dir: /data2/MLhep/ + #unmerged_tree_dir: [sim/train_339426/alice/cern.ch/user/a/alihyperloop/jobs/0130] #list of periods training + unmerged_tree_dir: [sim/train_341075/alice/cern.ch/user/a/alihyperloop/jobs/0132] #list of periods + pkl: [LHC23pp_mc_tuner_mb/prod_LHC24h1/pklmc] #list of periods + pkl_skimmed: [LHC23pp_mc_tuner_mb/prod_LHC24h1/pklskmc] #list of periods + pkl_skimmed_merge_for_ml: [LHC23pp_mc_tuner_mb/prod_LHC24h1/pklskmlmc] #list of periods + pkl_skimmed_merge_for_ml_all: LHC23pp_mc_tuner_mb/prod_LHC24_ana/mltotmc + pkl_evtcounter_all: LHC23pp_mc_tuner_mb/prod_LHC24_ana/evttotmc + #select_jobs: [[hy_396609], [hy_396597]] + mcreweights: [../Analyses] + ml: + evtsel: null + triggersel: + data: null + mc: null + + nclasses: [200000, 200000, 200000] + #nclasses: [16557, 16557, 16557] + equalise_sig_bkg: True + #mult_bkg: [30,2,2,3,3,5] + mult_bkg: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + sampletags: [0, 1, 1] + sel_bkg: fM < 2.22 or fM > 2.35 # for plotting significance; should agree with bkg selection in sel_ml + # best to have non-prompt (the smallest class) last, so the plots won't complain about the middle class missing + sel_ml: [fM < 2.22 or fM > 2.35, ismcsignal == 1 and ismcprompt == 1, ismcsignal == 1 and ismcfd == 1] + class_labels: [bkg, prompt, non-prompt] + nkfolds: 5 + rnd_shuffle: 12 + rnd_splt: 12 + rnd_all: 12 # Set to None for pure randomness + test_frac: 0.2 + binmin: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12] # must be equal to sel_skim_binmin (sel_skim_binmin bins) + binmax: [2, 3, 4, 5, 6, 7, 8, 10, 12, 24] # must be equal to sel_skim_binmax (sel_skim_binmin bins) + mltype: MultiClassification + ncorescrossval: 10 + prefix_dir_ml: /data2/ldellost/MLhep/ + mlplot: mlplot_Lc2023 # to be removed + mlout: mlout_Lc2023 # to be removed + + opt: + isFONLLfromROOT: true + filename_fonll: 'data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root' # file with FONLL predictions + fonll_particle: 'hLcpkpipred' + fonll_pred: 'max' # edge of the FONLL prediction + FF: 0.204 # fragmentation fraction + sigma_MB: 57.8e-3 # Minimum Bias cross section (pp) 50.87e-3 [b], 1 for Pb-Pb + Taa: 1 # 23260 [b^-1] in 0-10% Pb-Pb, 3917 [b^-1] in 30-50% Pb-Pb, 1 for pp + BR: 6.23e-2 # branching ratio of the decay Lc -> p K- pi+ + f_prompt: 0.9 # estimated fraction of prompt candidates + bkg_data_fraction: 0.05 # fraction of real data used in the estimation + num_steps: 111 # number of steps used in efficiency and signif. estimation + bkg_function: pol2 # fit function for bkg (among TH1 predefined fit functions, e.g. expo, pol1, pol2, ...) + save_fit: True # save bkg fits with the various cuts on ML output + raahp: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] # sel_skim_binmin bins + presel_gen_eff: "abs(fY) < 0.8" + #presel_gen_eff: "abs(fY) < 0.8 and abs(fPosZ) < 10" + + mlapplication: + data: + prefix_dir_app: /data2/ldellost/MLhep_newformat_mb_DD/ + pkl_skimmed_dec: [LHC23pp/MLapplication/prod_LHC23/skpkldecdata] #list of periods + pkl_skimmed_decmerged: [LHC23pp/MLapplication/prod_LHC23/skpkldecdatamerged] #list of periods + mc: + prefix_dir_app: /data2/ldellost/MLhep_newformat_mb_DD/ + pkl_skimmed_dec: [LHC23pp_mc/MLapplication_mult/prod_LHC24h1/skpkldecmc] #list of periods + pkl_skimmed_decmerged: [LHC23pp_mc/MLapplication_mult/prod_LHC24h1/skpkldecmcmerged] #list of periods + modelname: xgboost + modelsperptbin: [xgboost_classifierLcpKpi_dfselection_fPt_1.0_2.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_2.0_3.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_3.0_4.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_4.0_5.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_5.0_6.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_6.0_7.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_7.0_8.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_8.0_10.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_10.0_12.0.sav, xgboost_classifierLcpKpi_dfselection_fPt_12.0_24.0.sav] + probcutpresel: + data: [[0.05, 0.0, 0.0], [0.05, 0.0, 0.0], [0.08, 0.0, 0.0], [0.15, 0.0, 0.0], [0.20, 0.0, 0.0], [0.30, 0.0, 0.0], [0.40, 0.0, 0.0], [0.50, 0.0, 0.0], [0.60, 0.0, 0.0], [0.80, 0.0, 0.0]] #list of nbins + mc: [[0.05, 0.0, 0.0], [0.05, 0.0, 0.0], [0.08, 0.0, 0.0], [0.15, 0.0, 0.0], [0.20, 0.0, 0.0], [0.30, 0.0, 0.0], [0.40, 0.0, 0.0], [0.50, 0.0, 0.0], [0.60, 0.0, 0.0], [0.80, 0.0, 0.0]] #list of nbins + probcutoptimal: [[0.02, 0.40, 0.0], [0.03, 0.25, 0.0], [0.04, 0.25, 0.0], [0.07, 0.25, 0.0], [0.10, 0.25, 0.0], [0.11, 0.25, 0.0], [0.15, 0.25, 0.0], [0.18, 0.25, 0.0], [0.25, 0.25, 0.0], [0.35, 0.25, 0.0]] #list of nbins + #probcutoptimal: [[0.02, 0.0, 0.0], [0.03, 0.0, 0.0], [0.04, 0.0, 0.0], [0.07, 0.0, 0.0], [0.10, 0.0, 0.0], [0.11, 0.0, 0.0], [0.15, 0.0, 0.0], [0.18, 0.0, 0.0], [0.25, 0.0, 0.0], [0.35, 0.0, 0.0]] #list of nbins + + analysis: + anahptspectrum: "LctopKpi" #D0Kpi, DplusKpipi, DstarD0pi, DsKKpi, LctopKpi, LcK0Sp + fd_method: "dd" # fc, Nb, ext, dd, dd_N + crosssec_prompt: True # True for prompt, False for non-prompt + cctype: "pp" + sigmamb: 59.4e+9 # 50.87e+9 pp5TeV, 57.8e+9 pp13TeV, 59.4e+9 pp Run3, pb + inputfonllpred: data/fonll/CutVarLc_pp13TeV_LHC23_pass4_wide.root #data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root # data/fonll/CutVarLc_pp13TeV_LHC23_pass4_wide.root # + dir_general_plots: analysis_plots + + Run3analysis: + proc_type: Dhadrons + useperiod: [1] + plotbin: [1] + usesinglebineff: null + triggerbit: '' + + variations_db: database_variations_LcToPKPi_mult.yml + observables: + hptspectrum: + bins_gen_var: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 16, 24] + bins_det_var: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 16, 24] + label: "#it{p}_{T}" + label_y: "#Lambda_{c} Cross section" + + use_cuts: False + cuts: + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + + sel_an_binmin: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 16] + sel_an_binmax: [2, 3, 4, 5, 6, 7, 8, 10, 12, 16, 24] + binning_matching: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 9] + presel_gen_eff: "abs(fY) < 0.5" + evtsel: null + #evtsel: "abs(fPosZ)<10" + triggersel: + data: null + mc: null + weighttrig: false + + data: + runselection: [null] #FIXME + prefix_dir_res: "" + results: [/home/ldellost/MLhep_newformat_mb_DD/LHC23pp_forw/ResultsFit/prod_LHC23/resultsdata/default/default] #list of periods + resultsallp: /home/ldellost/MLhep_newformat_mb_DD/LHC23pp_forw/ResultsFit/resultsdatatot/default/default + mc: + runselection: [null] #FIXME + prefix_dir_res: "" + results: [/home/ldellost/MLhep_newformat_mb_DD/LHC23pp_mc_forw/ResultsFit/prod_LHC24h1/resultsmc/default/default] #list of periods + resultsallp: /home/ldellost/MLhep_newformat_mb_DD/LHC23pp_mc_forw/ResultsFit/prod_LHC23/resultsmctot/default/default + + fitcase: Lc + latexnamehadron: "#Lambda_{c}^{pK#pi}" + #nevents: 258442910841 + nevents: 290860860000 + dobkgfromsideband: false + mass_fit_lim: [2.10, 2.47] # region for the fit of the invariant mass distribution [GeV/c^2] + bin_width: 0.001 # bin width of the invariant mass histogram + n_rebin: [2, 2, 2, 2, 2, 3, 3, 4, 5, 6, 6] # number of mass bins to merge + + pdf_names: + pdf_sig: "sig" + pdf_bkg: "bkg" + param_names: + mass: "m" + gauss_mean: "mean" + gauss_sigma: "sigma_g1" + double_gauss_sigma: "sigma_wide" + alpha_l: "alpha1" + alpha_r: "alpha2" + n_l: "n1" + n_r: "n2" + fraction_refl: "frac_refl" + + # To initialize the individual fits in pT bins + # Decide whether to take the sigma from MC or data for individual fits + mass_roofit: + - level: mc + ptrange: [0., 3.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.2832,2.289], sigma_g1[.007,.007,.012])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.007, 0.007, 0.012], alpha1[1.0, 0.8, 2.0], n1[30., 30., 100.0], alpha2[1.0, 0.8, 2.0], n2[30., 30., 100.0])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [3., 4.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.2832,2.289], sigma_g1[.007,.007,.012])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.007, 0.007, 0.012], alpha1[1.0, 0.8, 2.0], n1[20., 20., 100.0], alpha2[1.0, 0.8, 2.0], n2[20., 20., 100.0])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [4., 7.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.012,.01,.016])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.008, 0.008, 0.018], alpha1[1.0, 0.8, 2.0], n1[6.0, 6.0, 100.0], alpha2[1.0, 0.8, 2.0], n2[6.0, 6.0, 100.0])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [7., 12.] + range: [2.10, 2.45] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.010, 0.010, 0.022], alpha1[1.0, 0.8, 2.0], n1[5.0, 5.0, 100.0], alpha2[1.0, 0.8, 2.0], n2[5.0, 5.0, 100.0])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [12., 24.] + range: [2.10, 2.45] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.020,.020,.029])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.016, 0.016, 0.029], alpha1[1.0, 0.8, 2.0], n1[6.0, 6.0, 100.0], alpha2[1.0, 0.8, 2.0], n2[6.0, 6.0, 100.0])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - ptrange: [0., 2.] + range: [2.216, 2.36] #2.216, 2.36 + #fix_params: ['n', 'f_peak'] + components: + #sig: + #fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[-1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [2., 3.] + range: [2.20, 2.37] + #fix_params: ['n', 'f_peak'] + components: + #sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[-1e8, 1e8], a1[-1e8, 1e8], a2[-1e8, 1e8]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [3., 4.] + range: [2.19, 2.38] + #fix_params: ['n', 'f_peak'] + components: + #sig: + #fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[-1e8, 1e8], a1[-1e8, 1e8], a2[-1e8, 1e8]})' + model: + fn: 'SUM::sum(f_sig[0., 0., 0.8]*sig, bkg)' + - ptrange: [4., 5.] + range: [2.19, 2.38] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[2000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [5., 6.] + range: [2.18, 2.39] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[5000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [6., 7.] + range: [2.16, 2.40] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[2000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [7., 8.] + range: [2.16, 2.40] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[1000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [8., 10.] + range: [2.1, 2.46] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[1000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - range: [2.1, 2.46] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[2000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + + systematics: + probvariation: + useperiod: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] #period from where to define prob cuts + ncutvar: 10 #number of looser and tighter variations + maxperccutvar: 0.25 #max diff in efficiency for loosest/tightest var + cutvarminrange: [[0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.7, 0.9], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3]] #Min starting point for scan + cutvarmaxrange: [[0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3]] #Max starting point for scan + fixedmean: true #Fix mean cutvar histo to central fit + fixedsigma: true #Fix sigma cutvar histo to central fit diff --git a/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_newformat_mult_ana.yml b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_mult.yml similarity index 99% rename from machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_newformat_mult_ana.yml rename to machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_mult.yml index 8863a172e8..1d0415c30a 100644 --- a/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_newformat_mult_ana.yml +++ b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_mult.yml @@ -369,11 +369,11 @@ LcpKpi: analysis: anahptspectrum: "LctopKpi" #D0Kpi, DplusKpipi, DstarD0pi, DsKKpi, LctopKpi, LcK0Sp - fd_method: "Nb" # fc, Nb, ext + fd_method: "Nb" # fc, Nb, ext, dd, dd_N crosssec_prompt: true # true for prompt, false for non-prompt cctype: "pp" sigmamb: 59.4e+9 # 50.87e+9 pp5TeV, 57.8e+9 pp13TeV, 59.4e+9 pp Run3, pb - inputfonllpred: data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root + inputfonllpred: data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root #data/fonll/CutVarLc_pp13TeV_LHC23_pass4_wide.root # dir_general_plots: analysis_plots Run3analysis_forward: diff --git a/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_mult_CrystalBall.yml b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_mult_CrystalBall.yml new file mode 100644 index 0000000000..a1f784c02f --- /dev/null +++ b/machine_learning_hep/data/data_run3/database_ml_parameters_LcToPKPi_mult_CrystalBall.yml @@ -0,0 +1,746 @@ +--- +############################################################################# +## © Copyright CERN 2023. All rights not expressly granted are reserved. ## +## Author: Gian.Michele.Innocenti@cern.ch ## +## This program is free software: you can redistribute it and/or modify it ## +## under the terms of the GNU General Public License as published by the ## +## Free Software Foundation, either version 3 of the License, or (at your ## +## option) any later version. This program is distributed in the hope that ## +## it will be useful, but WITHOUT ANY WARRANTY; without even the implied ## +## warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. ## +## See the GNU General Public License for more details. ## +## You should have received a copy of the GNU General Public License ## +## along with this program. if not, see . ## +############################################################################# + +LcpKpi: + doml: true + mass: 2.286 + sel_reco_unp: "fPt>0" + sel_gen_unp: "fPt>0" + sel_cen_unp: null + sel_good_evt_unp: null + #sel_reco_skim: ["fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015","fDecayLength > 0.015"] + sel_reco_skim: [null, null, null, null, null, null, null, null, null, null] + sel_gen_skim: [null, null, null, null, null, null, null, null, null, null] + sel_skim_binmin: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12] #list of nbins + sel_skim_binmax: [2, 3, 4, 5, 6, 7, 8, 10, 12, 24] #list of nbins + apply_yptacccut: false + var_binning: fPt + dofullevtmerge: false + var_cand: fCandidateSelFlag + var_swap: fIsCandidateSwapped + bitmap_sel: + var_name: fFlagMcMatchRec + var_name_gen: fFlagMcMatchGen + var_name_origgen: fOriginMcGen + var_name_origrec: fOriginMcRec + var_isstd: isstd + var_ismcsignal: ismcsignal + var_ismcprompt: ismcprompt + var_ismcfd: ismcfd + var_ismcbkg: ismcbkg + var_ismcrefl: ismcref + isstd: [[1], []] + ismcsignal: [[1], []] + ismcprompt: [[0], []] + ismcfd: [[1], []] + ismcbkg: [[], [1]] + ismcrefl: [[1], [1]] + + dfs: + read: + evtorig: + level: all + index: fIndexHFLCCOLLBASES + trees: + O2hflccollbase: [fNumContrib, fMultZeqNTracksPV, fCentFT0M] + extra: + fIsEventReject: 0 + + evtoriggen: + level: gen + index: fIndexHFLCMCCOLLBASES + trees: + O2hflcmccollbase: [fPosX, fPosY, fPosZ, fCentFT0M] + O2hflcmcrcollid: [fIndexArrayHFLCCOLLBASES] + + reco: + level: all + index: fIndexHfLcBases + trees: + O2hflcbase: [fIndexHFLCCOLLBASES, fPt, fY, fEta, fPhi, fM] + O2hflcpar: [fNProngsContributorsPV, fCpa, fCpaXY, fChi2PCA, fDecayLength, fDecayLengthXY, fPtProng0, fPtProng1, fPtProng2, fImpactParameter0, fImpactParameter1, fImpactParameter2, fNSigTpcPi0, fNSigTpcPr0, fNSigTpcKa1, fNSigTpcPi2, fNSigTpcPr2, fNSigTofPi0, fNSigTofPr0, fNSigTofKa1, fNSigTofPi2, fNSigTofPr2, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + O2hflcmc: + level: mc + vars: [fFlagMcMatchRec, fOriginMcRec, fIsCandidateSwapped] + O2hflcsel: + level: mc + vars: [fCandidateSelFlag] + filter: "fPt > 0." + #extra: + #fY: log((sqrt(2.28646**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(2.28646**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated + tags: + isstd: {var: fFlagMcMatchRec, req: [[1], []], level: mc} + ismcsignal: {var: fFlagMcMatchRec, req: [[1], []], abs: true, level: mc} + ismcbkg: {var: fFlagMcMatchRec, req: [[], [1]], abs: true, level: mc} + ismcprompt: {var: fOriginMcRec, req: [[0], []], level: mc} + ismcfd: {var: fOriginMcRec, req: [[1], []], level: mc} + swap: {cand: fCandidateSelFlag, var_swap: fIsCandidateSwapped, vars: [ismcsignal, ismcprompt, ismcfd], level: mc} + + gen: + level: mc + trees: + O2hflcpbase: [fIndexHFLCMCCOLLBASES, fPt, fY, fEta, fPhi, fFlagMcMatchGen, fOriginMcGen] + tags: + isstd: {var: fFlagMcMatchGen, req: [[1], []], level: mc} + ismcsignal: {var: fFlagMcMatchGen, req: [[1], []], abs: true, level: mc} + ismcbkg: {var: fFlagMcMatchGen, req: [[], [1]], abs: true, level: mc} + ismcprompt: {var: fOriginMcGen, req: [[0], []], level: mc} + ismcfd: {var: fOriginMcGen, req: [[1], []], level: mc} + #extra: + #fY: log((sqrt(2.28646**2 + (fPt * cosh(fEta))**2) + fPt * sinh(fEta)) / sqrt(2.28646**2 + fPt**2)) #TODO : change mass or make sure Lc mass is updated + + merge: + - {base: reco, ref: evtorig, extra: {fMultZeqNTracksPV_sub: fMultZeqNTracksPV - fNProngsContributorsPV}} + - {base: gen, ref: evtoriggen} + - {base: gen, ref: evtorig, left_on: fIndexArrayHFLCCOLLBASES, out: genrec} + # workaround for yamlfmt issue #110 + + write: + evtorig: + level: all + file: AnalysisResultsEvtOrig.parquet + evt: + level: all + source: evtorig + file: AnalysisResultsEvt.parquet + evtmc: + level: mc + source: evtoriggen + file: AnalysisResultsEvtGen.parquet + reco: + level: all + file: AnalysisResultsReco.parquet + gen: + level: mc + file: AnalysisResultsGenSl.parquet + genrec: + level: mc + file: AnalysisResultsGen.parquet + + variables: + var_all: [fIndexCollisions, fFlagMcMatchRec, fCandidateSelFlag, fOriginMcRec, fIsCandidateSwapped, fNProngsContributorsPV, fY, fEta, fPt, fCpa, fCpaXY, fM, fChi2PCA, fDecayLength, fDecayLengthXY, fPtProng0, fPtProng1, fPtProng2, fImpactParameter0, fImpactParameter1, fImpactParameter2, fNSigTpcPi0, fNSigTpcPr0, fNSigTpcKa1, fNSigTpcPi2, fNSigTpcPr2, fNSigTofPi0, fNSigTofPr0, fNSigTofKa1, fNSigTofPi2, fNSigTofPr2, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + var_jet: [fJetPt, fJetEta, fJetPhi] + var_jetsub: [fZg, fRg, fNsd] + var_jet_match: [df, fIndexHfCand2Prong] + var_jetsub_match: [df, fIndexLcChargedJets] + var_evt: + data: [fIndexCollisions, fPosX, fPosY, fPosZ, fNumContrib, fMultZeqNTracksPV, fCentFT0M] + mc: [fIndexCollisions, fPosX, fPosY, fPosZ, fMultZeqNTracksPV, fCentFT0M] + #mc: [fIndexCollisions, fPosX, fPosY, fPosZ, fNumContrib, fMultZeqNTracksPV, fCentFT0A, fCentFT0C, fCentFT0M, fCentFV0A] + var_gen: [fIndexMcCollisions, fPosX, fPosY, fPosZ, fPt, fY, fFlagMcMatchGen, fOriginMcGen] + var_training: + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + - [fImpactParameter0, fImpactParameter1, fImpactParameter2, fCpa, fChi2PCA, fDecayLength, fDecayLengthXY, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + var_selected: [fM, fY, fEta, fPt, fCpa, fCpaXY, fChi2PCA, fDecayLength, fDecayLengthXY, fPtProng0, fPtProng1, fPtProng2, fImpactParameter0, fImpactParameter1, fImpactParameter2, fNSigTpcPi0, fNSigTpcPr0, fNSigTpcKa1, fNSigTpcPi2, fNSigTpcPr2, fNSigTofPi0, fNSigTofPr0, fNSigTofKa1, fNSigTofPi2, fNSigTofPr2, fNSigTpcTofPi0, fNSigTpcTofPr0, fNSigTpcTofKa1, fNSigTpcTofPi2, fNSigTpcTofPr2] + var_boundaries: [fDecayLength, fPt] + var_correlation: + - [fDecayLength, fChi2PCA, fCpa] + - [fPt, fPt, fPt] + var_class: class + var_inv_mass: fM + var_y: fY + var_cuts: + - [fPtProng0, lt, null] + - [fPtProng1, lt, null] + - [fPtProng2, lt, null] + - [fCpa, lt, null] + - [fDecayLength, lt, null] + - [fChi2PCA, lt, null] + + plot_options: + prob_cut_scan: + pt_prong0: + xlim: + - 0 + - 8 + pt_prong1: + xlim: + - 0 + - 8 + pt_prong2: + xlim: + - 0 + - 8 + fDecayLength: + xlim: + - 0 + - 0.08 + fChi2PCA: + xlim: + - 0 + - 20. + fNSigTofPr0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(p)0" + fNSigTofPi0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(\\pi)0" + fNSigTofKa1: + xlim: [-10, 10] + xlabel: "n\\sigma_\\mathrm{TOF}(K)1" + fNSigTofPr2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(p)2" + fNSigTofPi2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TOF}(\\pi)2" + fNSigTpcPr0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(p)0" + fNSigTpcPi0: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(\\pi)0" + fNSigTpcKa1: + xlim: [-10, 10] + xlabel: "n\\sigma_\\mathrm{TPC}(K)1" + fNSigTpcPr2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(p)2" + fNSigTpcPi2: + xlim: [-50, 50] + xlabel: "n\\sigma_\\mathrm{TPC}(\\pi)2" + eff_cut_scan: + pt_prong0: + xlim: + - 0 + - 8 + pt_prong1: + xlim: + - 0 + - 8 + pt_prong2: + xlim: + - 0 + - 8 + fDecayLength: + xlim: + - 0 + - 0.08 + fChi2PCA: + xlim: + - 0 + - 20. + files_names: + namefile_unmerged_tree: AO2D.root + namefile_reco: AnalysisResultsReco.parquet + namefile_evt: AnalysisResultsEvt.parquet + namefile_collcnt: AnalysisResultsCollCnt.parquet + namefile_bccnt: AnalysisResultsBcCnt.parquet + namefile_evtvalroot: AnalysisResultsROOTEvtVal.root + namefile_evtorig: AnalysisResultsEvtOrig.parquet + namefile_gen: AnalysisResultsGen.parquet + namefile_gen_sl: AnalysisResultsGenSl.parquet + namefile_reco_applieddata: AnalysisResultsRecoAppliedData.parquet + namefile_reco_appliedmc: AnalysisResultsRecoAppliedMC.parquet + namefile_mcweights: mcweights.root + treeoutput: "Lctree" + histofilename: "masshisto.root" + efffilename: "effhisto.root" + respfilename: "resphisto.root" + crossfilename: "cross_section_tot.root" + + multi: + data: + nprocessesparallel: 80 + maxfiles: [-1] #list of periods + chunksizeunp: [100] #list of periods + chunksizeskim: [100] #list of periods + fracmerge: [0.05] #list of periods + seedmerge: [12] #list of periods + period: [LHC23] #list of periods + select_period: [1] + prefix_dir: /data2/MLhep/ + unmerged_tree_dir: [real/train_342351/alice/cern.ch/user/a/alihyperloop/jobs/0133] #list of periods + pkl: [LHC23pp/period_LHC23/pkldata] #list of periods + pkl_skimmed: [LHC23pp/period_LHC23/pklskdata] #list of periods + pkl_skimmed_merge_for_ml: [LHC23pp/period_LHC23/pklskmldata] #list of periods + pkl_skimmed_merge_for_ml_all: LHC23pp/mltotdata + pkl_evtcounter_all: LHC23pp/evttotdata + #select_jobs: [[hy_189959], [hy_189000]] + mcreweights: [../Analyses] + mc: + nprocessesparallel: 80 + maxfiles: [-1] #list of periods + chunksizeunp: [100] #list of periods + chunksizeskim: [100] #list of periods + fracmerge: [1.0] #list of periods + seedmerge: [1] #list of periods + period: [LHC24h1] #list of periods + select_period: [1] + prefix_dir: /data2/MLhep/ + unmerged_tree_dir: [sim/train_341075/alice/cern.ch/user/a/alihyperloop/jobs/0132] #list of periods + pkl: [LHC23pp_mc_tuner_mult/prod_LHC24h1/pklmc] #list of periods + pkl_skimmed: [LHC23pp_mc_tuner_mult/prod_LHC24h1/pklskmc] #list of periods + pkl_skimmed_merge_for_ml: [LHC23pp_mc_tuner_mult/prod_LHC24h1/pklskmlmc] #list of periods + pkl_skimmed_merge_for_ml_all: LHC23pp_mc_tuner_mult/prod_LHC24_ana/mltotmc + pkl_evtcounter_all: LHC23pp_mc_tuner_mult/prod_LHC24_ana/evttotmc + #select_jobs: [[hy_396609], [hy_396597]] + mcreweights: [../Analyses] + ml: + evtsel: null + triggersel: + data: null + mc: null + + nclasses: [200000, 200000, 200000] + equalise_sig_bkg: True + #mult_bkg: [30,2,2,3,3,5] + mult_bkg: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1] + sampletags: [0, 1, 1] + sel_bkg: fM < 2.22 or fM > 2.35 # for plotting significance; should agree with bkg selection in sel_ml + # best to have non-prompt (the smallest class) last, so the plots won't complain about the middle class missing + sel_ml: [fM < 2.22 or fM > 2.35, ismcsignal == 1 and ismcprompt == 1, ismcsignal == 1 and ismcfd == 1] + class_labels: [bkg, prompt, non-prompt] + nkfolds: 5 + rnd_shuffle: 12 + rnd_splt: 12 + rnd_all: 12 # Set to None for pure randomness + test_frac: 0.2 + binmin: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12] # must be equal to sel_skim_binmin (sel_skim_binmin bins) + binmax: [2, 3, 4, 5, 6, 7, 8, 10, 12, 24] # must be equal to sel_skim_binmax (sel_skim_binmin bins) + mltype: MultiClassification + ncorescrossval: 10 + prefix_dir_ml: /data2/ldellost/MLhep/ + mlplot: mlplot_Lc2023 # to be removed + mlout: mlout_Lc2023 # to be removed + + opt: + isFONLLfromROOT: true + filename_fonll: 'data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root' # file with FONLL predictions + fonll_particle: 'hLcpkpipred' + fonll_pred: 'max' # edge of the FONLL prediction + FF: 0.204 # fragmentation fraction + sigma_MB: 57.8e-3 # Minimum Bias cross section (pp) 50.87e-3 [b], 1 for Pb-Pb + Taa: 1 # 23260 [b^-1] in 0-10% Pb-Pb, 3917 [b^-1] in 30-50% Pb-Pb, 1 for pp + BR: 6.23e-2 # branching ratio of the decay Lc -> p K- pi+ + f_prompt: 0.9 # estimated fraction of prompt candidates + bkg_data_fraction: 0.05 # fraction of real data used in the estimation + num_steps: 111 # number of steps used in efficiency and signif. estimation + bkg_function: pol2 # fit function for bkg (among TH1 predefined fit functions, e.g. expo, pol1, pol2, ...) + save_fit: True # save bkg fits with the various cuts on ML output + raahp: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1] # sel_skim_binmin bins + presel_gen_eff: "abs(fY) < 0.8" + #presel_gen_eff: "abs(fY) < 0.8 and abs(fPosZ) < 10" + + mlapplication: + data: + prefix_dir_app: /data2/ldellost/MLhep_newformat/ + pkl_skimmed_dec: [LHC23pp/MLapplication/prod_LHC23/skpkldecdata] #list of periods + pkl_skimmed_decmerged: [LHC23pp/MLapplication/prod_LHC23/skpkldecdatamerged] #list of periods + mc: + prefix_dir_app: /data2/ldellost/MLhep_newformat/ + pkl_skimmed_dec: [LHC23pp_mc/MLapplication_mult/prod_LHC24h1/skpkldecmc] #list of periods + pkl_skimmed_decmerged: [LHC23pp_mc/MLapplication_mult/prod_LHC24h1/skpkldecmcmerged] #list of periods + modelname: xgboost + modelsperptbin: + - xgboost_classifierLcpKpi_dfselection_fPt_1.0_2.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_2.0_3.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_3.0_4.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_4.0_5.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_5.0_6.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_6.0_7.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_7.0_8.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_8.0_10.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_10.0_12.0.sav + - xgboost_classifierLcpKpi_dfselection_fPt_12.0_24.0.sav + probcutpresel: + data: [[0.05, 0.0, 0.0], [0.05, 0.0, 0.0], [0.08, 0.0, 0.0], [0.15, 0.0, 0.0], [0.20, 0.0, 0.0], [0.30, 0.0, 0.0], [0.40, 0.0, 0.0], [0.50, 0.0, 0.0], [0.60, 0.0, 0.0], [0.80, 0.0, 0.0]] #list of nbins + mc: [[0.05, 0.0, 0.0], [0.05, 0.0, 0.0], [0.08, 0.0, 0.0], [0.15, 0.0, 0.0], [0.20, 0.0, 0.0], [0.30, 0.0, 0.0], [0.40, 0.0, 0.0], [0.50, 0.0, 0.0], [0.60, 0.0, 0.0], [0.80, 0.0, 0.0]] #list of nbins + probcutoptimal: [[0.02, 0.40, 0.0], [0.03, 0.25, 0.0], [0.04, 0.25, 0.0], [0.07, 0.25, 0.0], [0.10, 0.25, 0.0], [0.11, 0.25, 0.0], [0.15, 0.25, 0.0], [0.18, 0.25, 0.0], [0.25, 0.25, 0.0], [0.35, 0.25, 0.0]] #list of nbins + + analysis: + anahptspectrum: "LctopKpi" #D0Kpi, DplusKpipi, DstarD0pi, DsKKpi, LctopKpi, LcK0Sp + fd_method: "dd" # fc, Nb, ext, dd, dd_N + crosssec_prompt: True # True for prompt, False for non-prompt + cctype: "pp" + sigmamb: 59.4e+9 # 50.87e+9 pp5TeV, 57.8e+9 pp13TeV, 59.4e+9 pp Run3, pb + inputfonllpred: data/fonll/CutVarLc_pp13TeV_LHC23_pass4_wide.root #data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root # + #inputfonllpred: data/fonll/DmesonLcPredictions_13TeV_y05_FFptDepLHCb_BRpythia8_PDG2020.root + dir_general_plots: analysis_plots + + Run3analysis_forward: + proc_type: Dhadrons_mult + useperiod: [1] + plotbin: [1, 1, 1, 1, 1, 1, 1] + usesinglebineff: null + fprompt_from_mb: false + corrEffMult: [false, true, true, true, true, true, true, true] + event_cand_validation: False + sel_binmin2: [0, 85, 70, 50, 30, 10, 1, 0] #list of var2 splittng nbins + sel_binmax2: [100, 100, 85, 70, 50, 30, 10, 1] + var_binning2: fCentFT0M + var_binning2_gen: fCentFT0M + var_binning2_weights: fMultZeqNTracksPV + mc_cut_on_binning2: false + signal_loss: true + signal_loss_idx: fIndexArrayHFLCCOLLBASES + nbinshisto: 100 + minvaluehisto: -0.0005 + maxvaluehisto: 100.0005 + triggerbit: '' + + variations_db: database_variations_LcToPKPi.yml + observables: + hptspectrum: + bins_gen_var: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 24] + bins_det_var: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12, 24] + label: "#it{p}_{T}" + label_y: "#Lambda_{c} Cross section" + + event_weighting_mc: + LHC24h1: + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_0_1 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_85_100 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_70_85 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_50_70 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_30_50 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_10_30 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_1_10 + according_to: fMultZeqNTracksPV + - filepath: data/event_weighting_mc/MultWeigths.root + histo_name: MultWeigths_0_1 + according_to: fMultZeqNTracksPV + + use_cuts: False + cuts: + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + - "fDecayLength > 0.02" + + sel_an_binmin: [1, 2, 3, 4, 5, 6, 7, 8, 10, 12] + sel_an_binmax: [2, 3, 4, 5, 6, 7, 8, 10, 12, 24] + binning_matching: [0, 1, 2, 3, 4, 5, 6, 7, 8, 9] + presel_gen_eff: "abs(fY) < 0.5" + evtsel: null + #evtsel: "abs(fPosZ)<10" + triggersel: + data: null + mc: null + weighttrig: false + + data: + runselection: [null] #FIXME + prefix_dir_res: /data2/ldellost/MLhep_newformat_fixed/ + results: [LHC23pp_forw/Results/prod_LHC23/resultsdata] #list of periods + resultsallp: LHC23pp_forw/Results/resultsdatatot + mc: + runselection: [null] #FIXME + prefix_dir_res: /data2/ldellost/MLhep_newformat_fixed/ + results: [LHC23pp_mc_forw/Results/prod_LHC24h1/resultsmc] #list of periods + resultsallp: LHC23pp_mc_forw/Results/prod_LHC23/resultsmctot + + fitcase: Lc + latexnamehadron: "#Lambda_{c}^{pK#pi}" + latexbin2var: "FT0M" + #nevents: 258442910841 + nevents: 290860860000 + dobkgfromsideband: false + mass_fit_lim: [2.10, 2.47] # region for the fit of the invariant mass distribution [GeV/c^2] + bin_width: 0.001 # bin width of the invariant mass histogram + n_rebin: [2, 2, 2, 2, 2, 3, 3, 4, 5, 6] + + pdf_names: + pdf_sig: "sig" + pdf_bkg: "bkg" + param_names: + mass: "m" + gauss_mean: "mean" + gauss_sigma: "sigma_g1" + double_gauss_sigma: "sigma_wide" + fraction_refl: "frac_refl" + alpha_l: "alpha1" + n_l: "n1" + + # To initialize the individual fits in pT bins + # Decide whether to take the sigma from MC or data for individual fits + mass_roofit: + - level: mc + ptrange: [0., 2.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.2832,2.289], sigma_g1[.007,.007,.012]) + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0085, 0.0085], alpha1[1.3506, 1.3506], n1[40.0972, 40.0972], alpha1[1.3506, 1.3506], n1[40.0972, 40.0972])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [2., 3.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.2832,2.289], sigma_g1[.007,.007,.012])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0088, 0.0088], alpha1[1.64, 1.64], n1[30., 30.], alpha1[1.64, 1.64], n1[30., 30.])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [3., 4.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.2832,2.289], sigma_g1[.007,.007,.012])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0102, 0.0102], alpha1[1.53, 1.53], n1[67.37, 67.37], alpha1[1.53, 1.53], n1[67.37, 67.37])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [4., 5.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.2832,2.289], sigma_g1[.007,.007,.012])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0112, 0.0112], alpha1[1.21, 1.21], n1[94.58, 94.58], alpha1[1.21, 1.21], n1[94.58, 94.58])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [5., 6.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.012,.01,.016])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0127, 0.0127], alpha1[1.147, 1.147], n1[77.2867, 77.2867], alpha1[1.147, 1.147], n1[77.2867, 77.2867])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [6., 7.] + range: [2.18, 2.39] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.286, 2.283,2.289], sigma_g1[.012,.01,.016]) + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0131, 0.0131], alpha1[1.026, 1.026], n1[90.26, 90.26], alpha1[1.026, 1.026], n1[90.26, 90.26])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [7., 8.] + range: [2.10, 2.45] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0151, 0.0151], alpha1[1.249, 1.249], n1[99.955, 99.955], alpha1[1.249, 1.249], n1[99.955, 99.955])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [8., 10.] + range: [2.10, 2.45] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0174, 0.0174], alpha1[1.315, 1.315], n1[4., 4.], alpha1[1.315, 1.315], n1[4., 4.])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [10., 12.] + range: [2.10, 2.45] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.013,.013,.020])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.0187, 0.0187], alpha1[1.013, 1.013], n1[5.0, 5.0], alpha1[1.013, 1.013], n1[5.0, 5.0])' + #fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.010, 0.010, 0.022], alpha1[1.0, 0.8, 2.0], n1[5.0, 5.0, 100.0], alpha2[1.0, 0.8, 2.0], n2[5.0, 5.0, 100.0])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - level: mc + ptrange: [12., 24.] + range: [2.10, 2.45] + components: + sig: + #fn: 'Gaussian::sig(m[1., 10], mean[2.2865, 2.283,2.289], sigma_g1[.020,.020,.029])' + fn: 'RooCrystalBall::sig(m[1., 10], mean[2.286, 2.283, 2.289], sigma_g1[0.020, 0.020], alpha1[0.8, 0.8], n1[80.01, 80.01], alpha1[0.8, 0.8], n1[80.01, 80.01])' + #wide: + # fn: 'Gaussian::wide(m, mean, expr("n*sigma_g1", n[1.,5.], sigma_g1))' + #model: + # fn: 'SUM::sig(f_peak[0.,1.]*peak, wide)' + bkg: + fn: 'Exponential::mcbkg(m, mcalpha[0.])' + model: + fn: 'SUM::mctot(mcfrac[0., 0., 1.0]*sig, mcbkg)' + - ptrange: [0., 2.] + range: [2.21, 2.36] #2.216, 2.36 + #fix_params: ['n', 'f_peak'] + components: + #sig: + #fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[-1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [2., 3.] + range: [2.20, 2.368] + #fix_params: ['n', 'f_peak'] + components: + #sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[-1e8, 1e8], a1[-1e8, 1e8], a2[-1e8, 1e8]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [3., 4.] + range: [2.19, 2.38] + #fix_params: ['n', 'f_peak'] + components: + #sig: + #fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[-1e8, 1e8], a1[-1e8, 1e8], a2[-1e8, 1e8]})' + model: + fn: 'SUM::sum(f_sig[0., 0., 0.8]*sig, bkg)' + - ptrange: [4., 5.] + range: [2.19, 2.38] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[2000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [5., 6.] + range: [2.18, 2.39] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.005,.015])' + bkg: + fn: 'Polynomial::bkg(m, {a0[5000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [6., 7.] + range: [2.16, 2.40] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[2000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [7., 8.] + range: [2.16, 2.40] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[1000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - ptrange: [8., 10.] + range: [2.1, 2.46] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[1000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + - range: [2.1, 2.46] + #fix_params: ['n', 'f_peak'] + components: + # sig: + # fn: 'Gaussian::sig(m, mean[2.28,2.29], sigma_g1[.005,.03])' + bkg: + fn: 'Polynomial::bkg(m, {a0[2000, -1e10, 1e10], a1[-1e10, 1e10], a2[-1e10, 1e10]})' + model: + fn: 'SUM::sum(f_sig[0.,1.]*sig, bkg)' + + systematics: + probvariation: + useperiod: [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] #period from where to define prob cuts + ncutvar: 10 #number of looser and tighter variations + maxperccutvar: 0.25 #max diff in efficiency for loosest/tightest var + cutvarminrange: [[0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.7, 0.9], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3]] #Min starting point for scan + cutvarmaxrange: [[0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.7, 0.9], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3], [0.1, 0.3]] #Max starting point for scan + fixedmean: true #Fix mean cutvar histo to central fit + fixedsigma: true #Fix sigma cutvar histo to central fit diff --git a/machine_learning_hep/hf_analysis_utils.py b/machine_learning_hep/hf_analysis_utils.py index 4414b841aa..5314488d87 100644 --- a/machine_learning_hep/hf_analysis_utils.py +++ b/machine_learning_hep/hf_analysis_utils.py @@ -16,7 +16,7 @@ file: hf_analysis_utils.py brief: script with miscellanea utils methods for the HF analyses author: Fabrizio Grosa , CERN -Macro committed and manteined in O2Physics: +Macro committed and mantained in O2Physics: https://github.com/AliceO2Group/O2Physics/tree/master/PWGHF/D2H/Macros """ @@ -62,6 +62,11 @@ def compute_crosssection( if rawy <= 0: crosssection = -9999 crosssec_unc = -1 + elif method_frac == "dd_N": + crosssection = frac * sigma_mb / (2 * delta_pt * delta_y * n_events * b_ratio) + # TODO: How to calculate the uncertainty? + # frac_unc / frac * crosssection? + crosssec_unc = 0.0 else: crosssection = rawy * frac * sigma_mb / (2 * delta_pt * delta_y * eff_times_acc * n_events * b_ratio) if method_frac in ("Nb", "ext"): @@ -124,8 +129,8 @@ def compute_fraction_fc( frac_fd = [frac_fd_cent, frac_fd_cent, frac_fd_cent] return frac_prompt, frac_fd - for i_sigma, (sigma_p, sigma_f) in enumerate(zip(cross_sec_prompt, cross_sec_fd)): - for i_raa, (raa_p, raa_f) in enumerate(zip(raa_prompt, raa_fd)): + for i_sigma, (sigma_p, sigma_f) in enumerate(zip(cross_sec_prompt, cross_sec_fd, strict=False)): + for i_raa, (raa_p, raa_f) in enumerate(zip(raa_prompt, raa_fd, strict=False)): if i_sigma == 0 and i_raa == 0: frac_prompt_cent = 1.0 / (1 + acc_eff_fd / acc_eff_prompt * sigma_f / sigma_p * raa_f / raa_p) frac_fd_cent = 1.0 / (1 + acc_eff_prompt / acc_eff_fd * sigma_p / sigma_f * raa_p / raa_f) @@ -218,34 +223,20 @@ def compute_fraction_nb( frac_cent * rawy * sigma_mb / 2 / acc_eff_same / delta_pt / delta_y / b_ratio / n_events ) delta_raa = abs((raa_other - raa_other_old) / raa_other) - else: - if raa_rat == 1.0 and taa == 1.0: # pp - frac.append( - 1 - sigma * delta_pt * delta_y * acc_eff_other * b_ratio * n_events * 2 / rawy / sigma_mb + elif raa_rat == 1.0 and taa == 1.0: # pp + frac.append(1 - sigma * delta_pt * delta_y * acc_eff_other * b_ratio * n_events * 2 / rawy / sigma_mb) + else: # p-Pb or Pb-Pb: iterative evaluation of Raa needed + delta_raa = 1.0 + frac_tmp = 1.0 + while delta_raa > 1.0e-3: + raw_fd = ( + taa * raa_rat * raa_other * sigma * delta_pt * delta_y * acc_eff_other * b_ratio * n_events * 2 ) - else: # p-Pb or Pb-Pb: iterative evaluation of Raa needed - delta_raa = 1.0 - frac_tmp = 1.0 - while delta_raa > 1.0e-3: - raw_fd = ( - taa - * raa_rat - * raa_other - * sigma - * delta_pt - * delta_y - * acc_eff_other - * b_ratio - * n_events - * 2 - ) - frac_tmp = 1 - raw_fd / rawy - raa_other_old = raa_other - raa_other = ( - frac_tmp * rawy * sigma_mb / 2 / acc_eff_same / delta_pt / delta_y / b_ratio / n_events - ) - delta_raa = abs((raa_other - raa_other_old) / raa_other) - frac.append(frac_tmp) + frac_tmp = 1 - raw_fd / rawy + raa_other_old = raa_other + raa_other = frac_tmp * rawy * sigma_mb / 2 / acc_eff_same / delta_pt / delta_y / b_ratio / n_events + delta_raa = abs((raa_other - raa_other_old) / raa_other) + frac.append(frac_tmp) if frac: frac.sort() @@ -256,6 +247,44 @@ def compute_fraction_nb( return frac +def compute_fraction_dd(acc_eff_same, acc_eff_other, corryields_same, corryields_other, cov_same, cov_other, cov_comb): + """ + Method to get fraction of prompt / FD fraction with data-driven method + Parameters + ---------- + - acc_eff_same: efficiency times acceptance of prompt (non-prompt) D + - acc_eff_other: efficiency times acceptance of non-prompt (prompt) D + - corryields_same: corrected yield of prompt (non-prompt) D computed with cut variation + - corryields_other: corrected yield of non-prompt (prompt) D computed with cut variation + - cov_same: covariance of prompt (non-prompt) D computed with cut variation + - cov_other: covariance of non-prompt (prompt) D computed with cut variation + - cov_comb: covariance of prompt and non-prompt D computed with cut variation + Returns + ---------- + - frac: list of fraction of prompt (non-prompt) D (central, min, max) + where min = frac - unc, max = frac + unc, and unc is the statistical + uncertainty propagated from efficiency and corrected yield uncertainties + """ + yield_times_acceff_same = corryields_same * acc_eff_same + yield_times_acceff_other = corryields_other * acc_eff_other + frac_v = yield_times_acceff_same / (yield_times_acceff_same + yield_times_acceff_other) + # print(f"same yield times acceff: {yield_times_acceff_same} " \ + # f"other {yield_times_acceff_other} final frac: {frac_v}") + + denom = (yield_times_acceff_same + yield_times_acceff_other) ** 2 + der_same_same = ( + acc_eff_same * (yield_times_acceff_same + yield_times_acceff_other) - acc_eff_same**2 * corryields_same + ) / denom + der_same_other = -acc_eff_same * acc_eff_other * corryields_same / denom + unc = np.sqrt( + der_same_same**2 * cov_same + der_same_other**2 * cov_other + 2 * der_same_same * der_same_other * cov_comb + ) + # print(f"denom {denom} der_same_same {der_same_same} der_same_other {der_same_other} " \ + # f"cov same {cov_same} cov other {cov_other} cov comb {cov_comb} final unc {unc}") + + return [frac_v, frac_v - unc, frac_v + unc] + + def get_hist_binlimits(histo): """ Method to retrieve bin limits of ROOT.TH1 diff --git a/machine_learning_hep/hf_pt_spectrum.py b/machine_learning_hep/hf_pt_spectrum.py index ac1a6e2f4b..1db777dc7c 100644 --- a/machine_learning_hep/hf_pt_spectrum.py +++ b/machine_learning_hep/hf_pt_spectrum.py @@ -18,7 +18,7 @@ usage: python3 HfPtSpectrum.py CONFIG authors: Fabrizio Grosa , CERN Luigi Dello Stritto , CERN -Macro committed and manteined in O2Physics: +Macro committed and mantained in O2Physics: https://github.com/AliceO2Group/O2Physics/tree/master/PWGHF/D2H/Macros """ @@ -40,6 +40,7 @@ from machine_learning_hep.hf_analysis_utils import ( # pylint: disable=import-error compute_crosssection, + compute_fraction_dd, compute_fraction_fc, compute_fraction_nb, get_hist_binlimits, @@ -49,7 +50,7 @@ def hf_pt_spectrum( channel, # pylint: disable=too-many-locals, too-many-arguments, too-many-statements, too-many-branches b_ratio, - inputfonllpred, + input_fonll_or_fdd_pred, frac_method, prompt_frac, eff_filename, @@ -59,6 +60,7 @@ def hf_pt_spectrum( yield_histoname, norm, sigmamb, + crosssec_prompt, output_prompt, output_file, ): @@ -89,7 +91,7 @@ def hf_pt_spectrum( print(f"\033[91mERROR: channel {channel} not supported. Exit\033[0m") sys.exit(2) - if frac_method not in ["Nb", "fc", "ext"]: + if frac_method not in ["Nb", "fc", "ext", "dd", "dd_N"]: print(f"\033[91mERROR: method to subtract nonprompt {frac_method} not supported. Exit\033[0m") sys.exit(5) @@ -104,16 +106,24 @@ def hf_pt_spectrum( histos = {} - histos["FONLL"] = {"prompt": {}, "nonprompt": {}} - infile_fonll = TFile.Open(inputfonllpred) - for pred in ("central", "min", "max"): - histos["FONLL"]["nonprompt"][pred] = infile_fonll.Get(f"{fonll_hist_name[channel]}fromBpred_{pred}_corr") - histos["FONLL"]["nonprompt"][pred].SetDirectory(0) - if frac_method == "fc": - histos["FONLL"]["prompt"][pred] = infile_fonll.Get(f"{fonll_hist_name[channel]}pred_{pred}") - histos["FONLL"]["prompt"][pred].SetDirectory(0) + infile_pred = TFile.Open(input_fonll_or_fdd_pred) + if frac_method in ("dd", "dd_N"): + histos["corryields_fdd"] = [infile_pred.Get("hCorrYieldsPrompt"), infile_pred.Get("hCorrYieldsNonPrompt")] + histos["covariances"] = [ + infile_pred.Get("hCovPromptPrompt"), + infile_pred.Get("hCovNonPromptNonPrompt"), + infile_pred.Get("hCovPromptNonPrompt"), + ] + else: + histos["FONLL"] = {"prompt": {}, "nonprompt": {}} + for pred in ("central", "min", "max"): + histos["FONLL"]["nonprompt"][pred] = infile_pred.Get(f"{fonll_hist_name[channel]}fromBpred_{pred}_corr") + histos["FONLL"]["nonprompt"][pred].SetDirectory(0) + if frac_method == "fc": + histos["FONLL"]["prompt"][pred] = infile_pred.Get(f"{fonll_hist_name[channel]}pred_{pred}") + histos["FONLL"]["prompt"][pred].SetDirectory(0) - infile_fonll.Close() + infile_pred.Close() infile_rawy = TFile.Open(yield_filename) histos["rawyields"] = infile_rawy.Get(yield_histoname) @@ -151,7 +161,7 @@ def hf_pt_spectrum( axistit_cross = "d#sigma/d#it{p}_{T} (pb GeV^{-1} #it{c})" axistit_cross_times_br = "d#sigma/d#it{p}_{T} #times BR (pb GeV^{-1} #it{c})" axistit_pt = "#it{p}_{T} (GeV/#it{c})" - axistit_fprompt = "#if{f}_{prompt}" + axistit_fprompt = "#it{f}_{prompt}" gfraction = TGraphAsymmErrors() gfraction.SetNameTitle("gfraction", f";{axistit_pt};{axistit_fprompt}") @@ -169,19 +179,21 @@ def hf_pt_spectrum( ) hnorm = TH1F("hnorm", "hnorm", 1, 0, 1) - for i_pt, (ptmin, ptmax) in enumerate(zip(ptlims["rawyields"][:-1], ptlims["rawyields"][1:])): + crosssec_nonprompt_fonll = [] + for i_pt, (ptmin, ptmax) in enumerate(zip(ptlims["rawyields"][:-1], ptlims["rawyields"][1:], strict=False)): pt_cent = (ptmax + ptmin) / 2 pt_delta = ptmax - ptmin rawy = histos["rawyields"].GetBinContent(i_pt + 1) rawy_unc = histos["rawyields"].GetBinError(i_pt + 1) eff_times_acc_prompt = histos["acceffp"].GetBinContent(i_pt + 1) eff_times_acc_nonprompt = histos["acceffnp"].GetBinContent(i_pt + 1) - ptmin_fonll = histos["FONLL"]["nonprompt"]["central"].GetXaxis().FindBin(ptmin * 1.0001) - ptmax_fonll = histos["FONLL"]["nonprompt"]["central"].GetXaxis().FindBin(ptmax * 0.9999) - crosssec_nonprompt_fonll = [ - histos["FONLL"]["nonprompt"][pred].Integral(ptmin_fonll, ptmax_fonll, "width") / (ptmax - ptmin) - for pred in histos["FONLL"]["nonprompt"] - ] + if frac_method not in ("dd", "dd_N"): + ptmin_fonll = histos["FONLL"]["nonprompt"]["central"].GetXaxis().FindBin(ptmin * 1.0001) + ptmax_fonll = histos["FONLL"]["nonprompt"]["central"].GetXaxis().FindBin(ptmax * 0.9999) + crosssec_nonprompt_fonll = [ + histos["FONLL"]["nonprompt"][pred].Integral(ptmin_fonll, ptmax_fonll, "width") / (ptmax - ptmin) + for pred in histos["FONLL"]["nonprompt"] + ] # compute prompt fraction frac = [0, 0, 0] @@ -210,13 +222,34 @@ def hf_pt_spectrum( ) elif frac_method == "ext": frac[0] = prompt_frac[i_pt] + elif frac_method == "dd": + eff_times_acc_own = eff_times_acc_prompt if crosssec_prompt else eff_times_acc_nonprompt + eff_times_acc_other = eff_times_acc_nonprompt if crosssec_prompt else eff_times_acc_prompt + pnp_ind = 0 if crosssec_prompt else 1 + frac = compute_fraction_dd( + eff_times_acc_own, + eff_times_acc_other, + histos["corryields_fdd"][pnp_ind].GetBinContent(i_pt + 1), + histos["corryields_fdd"][1 - pnp_ind].GetBinContent(i_pt + 1), + histos["covariances"][pnp_ind].GetBinContent(i_pt + 1), + histos["covariances"][1 - pnp_ind].GetBinContent(i_pt + 1), + histos["covariances"][2].GetBinContent(i_pt + 1), + ) + print("FRACTION HEREEEEEEEE: ") + print(histos["corryields_fdd"][pnp_ind].GetBinContent(i_pt + 1)) + print(histos["corryields_fdd"][1 - pnp_ind].GetBinContent(i_pt + 1)) + print("FRACTION HEREEEEEEEE: ") + elif frac_method == "dd_N": + pnp_ind = 0 if crosssec_prompt else 1 + frac = [histos["corryields_fdd"][pnp_ind].GetBinContent(i_pt + 1)] * 3 # compute cross section times BR + eff_times_acc = eff_times_acc_prompt if crosssec_prompt else eff_times_acc_nonprompt crosssec, crosssec_unc = compute_crosssection( rawy, rawy_unc, frac[0], - eff_times_acc_prompt, + eff_times_acc, ptmax - ptmin, 1.0, sigmamb,