From 98d636924c601eb43d4015454eaa6e978ea27c82 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Fri, 8 Nov 2024 17:28:28 +0100 Subject: [PATCH 1/6] [WIP] kde producer --- run3/flow/kde_producer.py | 102 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 102 insertions(+) create mode 100644 run3/flow/kde_producer.py diff --git a/run3/flow/kde_producer.py b/run3/flow/kde_producer.py new file mode 100644 index 00000000..14663387 --- /dev/null +++ b/run3/flow/kde_producer.py @@ -0,0 +1,102 @@ +import argparse +import yaml +import pandas as pd +import numpy as np +import uproot +import ROOT +import ctypes +from ROOT import TFile, TKDE, TCanvas, TH1D + +def kde_producer(tree, var, cent, pt_min, pt_max, inv_mass_bins, flag, outfile): + + # convert the tree to a pandas dataframe + dfsData = [] + # print(tree) + with uproot.open(f'{tree}') as f: + # print(f.keys()) + for iKey, key in enumerate(f.keys()): + if 'O2hfcanddplite' in key: + # print(key) + dfData = f[key].arrays(library='pd') + dfsData.append(dfData) + + full_dataset = pd.concat([df for df in dfsData], ignore_index=True) + print(full_dataset['fFlagMcMatchRec'].tolist()) + print(full_dataset['fOriginMcRec'].tolist()) + print(full_dataset['fFlagMcDecayChanRec'].tolist()) + pt_filtered_df = full_dataset.query(f"{pt_min} <= fPt < {pt_max}") + filtered_df = pt_filtered_df.query(f"fFlagMcMatchRec == {flag} or fFlagMcMatchRec == {-flag}") + # filtered_df = pt_filtered_df.query(f"fFlagMcMatchRec == 1") + var_values = filtered_df[f'{var}'].tolist() # Or use `.tolist()` to get a list + # print(var_values) + + kde = TKDE(len(var_values), np.asarray(var_values, 'd'), 1.7, 2.0) + kde_func = kde.GetFunction(1000) + + binned_var_values = TH1D(f'hBinned_{var}_{flag}', f'hBinned_{var}_{flag}', 5000, 1, 3) + for var_value in var_values: + binned_var_values.Fill(var_value) + + cOverlap = TCanvas('cOverlap', 'cOverlap', 600, 600) + cOverlap.cd() + binned_var_values.Draw() + kde_func.Draw('SAME') + + outfile.mkdir(f'pT_{pt_low}_{pt_max}') + outfile.cd(f'pT_{pt_low}_{pt_max}') + kde_func.Write() + binned_var_values.Write() + cOverlap.Write() + + return kde_func, binned_var_values + + +if __name__ == "__main__": + parser = argparse.ArgumentParser(description="Arguments") + parser.add_argument("--var", "-v", metavar="text", + default="fM", help="variable of interest") + parser.add_argument("--config", "-cfg", metavar="text", + default="config.yaml", help="configuration file") + parser.add_argument("--centrality", "-c", metavar="text", + default="k2050", help="centrality class") + parser.add_argument("--ptmin", "-pmin", metavar="text", + default="2.", help="min pt") + parser.add_argument("--ptmax", "-pmax", metavar="text", + default="4.", help="max pt") + parser.add_argument("--invmassbins", "-imbins", metavar="text", + default="[1, 1.5, 2]", help="inv mass bins") + parser.add_argument("--flag", "-f", metavar="chn flag", + default="2", help="channel flag") + parser.add_argument("--input", "-in", metavar="path/input.root", + default="AnalysisResults.root", help="path to tree") + parser.add_argument("--wagon_id", "-w", metavar="text", + default="", help="wagon ID", required=False) + parser.add_argument("--outputdir", "-o", metavar="text", + default=".", help="output directory") + parser.add_argument("--suffix", "-s", metavar="text", + default="", help="suffix for output files") + + args = parser.parse_args() + print(args) + + KDEs = [] + binned_histos = [] + if args.config != parser.get_default("config"): + with open(args.config, 'r') as ymlCfgFile: + config = yaml.load(ymlCfgFile, yaml.FullLoader) + + output_dir = config["outputdir"] if args.outputdir == parser.get_default("outputdir") else args.outputdir + suffix = config["suffix"] if args.suffix == parser.get_default("suffix") else args.suffix + outfile = ROOT.TFile(f'{output_dir}/kde_{suffix}.root', 'RECREATE') + + if args.config != parser.get_default("config"): + for pt_low, pt_max, inv_mass_bins in zip(config['pt_mins'], config['pt_maxs'], config['inv_mass_bins']): + KDE, histo = kde_producer(config['input'], config['variable'], config['centrality'], \ + pt_low, pt_max, inv_mass_bins, config['chn_flag'], + outfile) + else: + KDE, histo = kde_producer(args.input, args.var, args.centrality, \ + args.ptmin, args.ptmax, args.invmassbins, args.flag, + outfile) + outfile.Close() + \ No newline at end of file From 324eea02a285b0a89ac72434ec9bb11c0a407553 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Wed, 13 Nov 2024 09:38:47 +0100 Subject: [PATCH 2/6] adjust spaces --- run3/flow/invmassfitter/InvMassFitter.h | 114 ++++++++-------- run3/flow/invmassfitter/VnVsMassFitter.h | 164 +++++++++++------------ 2 files changed, 139 insertions(+), 139 deletions(-) diff --git a/run3/flow/invmassfitter/InvMassFitter.h b/run3/flow/invmassfitter/InvMassFitter.h index 41e8b11d..58615fd1 100644 --- a/run3/flow/invmassfitter/InvMassFitter.h +++ b/run3/flow/invmassfitter/InvMassFitter.h @@ -187,63 +187,63 @@ class InvMassFitter : public TNamed { Bool_t PrepareHighPolFit(TF1 *fback); Double_t BackFitFuncPolHelper(Double_t *x,Double_t *par); - TH1F* fHistoInvMass; /// histogram to fit - Double_t fMinMass; /// lower mass limit - Double_t fMaxMass; /// upper mass limit - Int_t fTypeOfFit4Bkg; /// background fit func - Int_t fPolDegreeBkg; /// degree of polynomial expansion for back fit (option 6 for back) - Int_t fCurPolDegreeBkg; /// help variable - Double_t fMassParticle; /// pdg value of particle mass - Int_t fTypeOfFit4Sgn; /// signal fit func - Double_t fMass; /// signal gaussian mean value - Double_t fMassErr; /// unc on signal gaussian mean value - Double_t fMassLowerLim; /// lower limit of the allowed mass range - Double_t fMassUpperLim; /// upper limit of the allowed mass range - Bool_t fBoundMean; /// switch for bound mean of gaussian - Double_t fSigmaSgn; /// signal gaussian sigma - Double_t fSigmaSgnErr; /// unc on signal gaussian sigma - Double_t fSigmaSgn2Gaus; /// signal second gaussian sigma in case of k2Gaus - Bool_t fFixedMean; /// switch for fix mean of gaussian - Bool_t fFixedSigma; /// switch for fix Sigma of gaussian - Bool_t fBoundSigma; /// switch for bound Sigma of gaussian - Double_t fSigmaVar; /// value of bound Sigma of gaussian - Double_t fParSig; /// +/- range variation of bound Sigma of gaussian in % - Bool_t fFixedSigma2Gaus; /// switch for fix Sigma of second gaussian in case of k2Gaus - Double_t fFixedRawYield; /// initialization for wa yield - Double_t fFrac2Gaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar - Bool_t fFixedFrac2Gaus; /// switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar - Double_t fRatio2GausSigma; /// initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar - Bool_t fFixedRatio2GausSigma; /// switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar - Int_t fNParsSig; /// fit parameters in signal fit function - Int_t fNParsBkg; /// fit parameters in background fit function - Bool_t fOnlySideBands; /// kTRUE = only side bands considered - Double_t fNSigma4SideBands; /// number of sigmas to veto the signal peak - Bool_t fCheckSignalCountsAfterFirstFit; /// switch for check after first fit - TString fFitOption; /// L, LW or Chi2 - Double_t fRawYield; /// signal gaussian integral - Double_t fRawYieldErr; /// err on signal gaussian integral - TF1* fSigFunc; /// Signal fit function - TF1* fBkgFuncSb; /// background fit function (1st step, side bands only) - TF1* fBkgFunc; /// background fit function (1st step, extended in peak region) - TF1* fBkgFuncRefit; /// background fit function (2nd step) - Bool_t fReflections; /// flag use/not use reflections - Int_t fNParsRfl; /// fit parameters in reflection fit function - Double_t fRflOverSig; /// reflection/signal - Bool_t fFixRflOverSig; /// switch for fix refl/signal - TH1F* fHistoTemplRfl; /// histogram with reflection template - Bool_t fSmoothRfl; /// switch for smoothing of reflection template - Double_t fRawYieldHelp; /// internal variable for fit with reflections - TF1* fRflFunc; /// fit function for reflections - TF1* fBkRFunc; /// fit function for reflections - Bool_t fSecondPeak; /// switch off/on second peak (for D+->KKpi in Ds) - Int_t fNParsSec; /// fit parameters in 2nd peak fit function - Double_t fSecMass; /// position of the 2nd peak - Double_t fSecWidth; /// width of the 2nd peak - Bool_t fFixSecMass; /// flag to fix the position of the 2nd peak - Bool_t fFixSecWidth; /// flag to fix the width of the 2nd peak - TF1* fSecFunc; /// fit function for second peak - TF1* fTotFunc; /// total fit function - Bool_t fAcceptValidFit; /// accept a fit when IsValid() gives true, nevertheless the status code + TH1F* fHistoInvMass; /// histogram to fit + Double_t fMinMass; /// lower mass limit + Double_t fMaxMass; /// upper mass limit + Int_t fTypeOfFit4Bkg; /// background fit func + Int_t fPolDegreeBkg; /// degree of polynomial expansion for back fit (option 6 for back) + Int_t fCurPolDegreeBkg; /// help variable + Double_t fMassParticle; /// pdg value of particle mass + Int_t fTypeOfFit4Sgn; /// signal fit func + Double_t fMass; /// signal gaussian mean value + Double_t fMassErr; /// unc on signal gaussian mean value + Double_t fMassLowerLim; /// lower limit of the allowed mass range + Double_t fMassUpperLim; /// upper limit of the allowed mass range + Bool_t fBoundMean; /// switch for bound mean of gaussian + Double_t fSigmaSgn; /// signal gaussian sigma + Double_t fSigmaSgnErr; /// unc on signal gaussian sigma + Double_t fSigmaSgn2Gaus; /// signal second gaussian sigma in case of k2Gaus + Bool_t fFixedMean; /// switch for fix mean of gaussian + Bool_t fFixedSigma; /// switch for fix Sigma of gaussian + Bool_t fBoundSigma; /// switch for bound Sigma of gaussian + Double_t fSigmaVar; /// value of bound Sigma of gaussian + Double_t fParSig; /// +/- range variation of bound Sigma of gaussian in % + Bool_t fFixedSigma2Gaus; /// switch for fix Sigma of second gaussian in case of k2Gaus + Double_t fFixedRawYield; /// initialization for wa yield + Double_t fFrac2Gaus; /// initialization for fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar + Bool_t fFixedFrac2Gaus; /// switch for fixed fraction of 2nd gaussian in case of k2Gaus or k2GausSigmaRatioPar + Double_t fRatio2GausSigma; /// initialization for ratio between two gaussian sigmas in case of k2GausSigmaRatioPar + Bool_t fFixedRatio2GausSigma; /// switch for fixed ratio between two gaussian sigmas in case of k2GausSigmaRatioPar + Int_t fNParsSig; /// fit parameters in signal fit function + Int_t fNParsBkg; /// fit parameters in background fit function + Bool_t fOnlySideBands; /// kTRUE = only side bands considered + Double_t fNSigma4SideBands; /// number of sigmas to veto the signal peak + Bool_t fCheckSignalCountsAfterFirstFit; /// switch for check after first fit + TString fFitOption; /// L, LW or Chi2 + Double_t fRawYield; /// signal gaussian integral + Double_t fRawYieldErr; /// err on signal gaussian integral + TF1* fSigFunc; /// Signal fit function + TF1* fBkgFuncSb; /// background fit function (1st step, side bands only) + TF1* fBkgFunc; /// background fit function (1st step, extended in peak region) + TF1* fBkgFuncRefit; /// background fit function (2nd step) + Bool_t fReflections; /// flag use/not use reflections + Int_t fNParsRfl; /// fit parameters in reflection fit function + Double_t fRflOverSig; /// reflection/signal + Bool_t fFixRflOverSig; /// switch for fix refl/signal + TH1F* fHistoTemplRfl; /// histogram with reflection template + Bool_t fSmoothRfl; /// switch for smoothing of reflection template + Double_t fRawYieldHelp; /// internal variable for fit with reflections + TF1* fRflFunc; /// fit function for reflections + TF1* fBkRFunc; /// fit function for background and reflections + Bool_t fSecondPeak; /// switch off/on second peak (for D+->KKpi in Ds) + Int_t fNParsSec; /// fit parameters in 2nd peak fit function + Double_t fSecMass; /// position of the 2nd peak + Double_t fSecWidth; /// width of the 2nd peak + Bool_t fFixSecMass; /// flag to fix the position of the 2nd peak + Bool_t fFixSecWidth; /// flag to fix the width of the 2nd peak + TF1* fSecFunc; /// fit function for second peak + TF1* fTotFunc; /// total fit function + Bool_t fAcceptValidFit; /// accept a fit when IsValid() gives true, nevertheless the status code /// \cond CLASSIMP ClassDef(InvMassFitter,9); /// class for invariant mass fit diff --git a/run3/flow/invmassfitter/VnVsMassFitter.h b/run3/flow/invmassfitter/VnVsMassFitter.h index e88009cc..cb4c7653 100755 --- a/run3/flow/invmassfitter/VnVsMassFitter.h +++ b/run3/flow/invmassfitter/VnVsMassFitter.h @@ -178,88 +178,88 @@ class VnVsMassFitter : public TObject { void SetParNames(); ///data members - TH1F* fMassHisto; /// mass histogram to fit - TH1F* fVnVsMassHisto; /// vn vs. mass histogram to fit - Int_t fMassSgnFuncType; /// type of mass signal fit function - Int_t fMassBkgFuncType; /// type of mass bkg fit function - Int_t fVnBkgFuncType; /// type of vn bkg fit function - TF1* fMassFuncFromPrefit; /// mass fit function (1st step, from prefit) - TF1* fMassBkgFunc; /// mass bkg fit function (final, after simultaneus fit) - TF1* fMassSgnFunc; /// mass signal fit function (final, after simultaneus fit) - TF1* fMassTotFunc; /// mass fit function (final, after simultaneus fit) - TF1* fVnBkgFuncSb; /// vn bkg fit function (1st step from SB prefit) - TF1* fVnBkgFunc; /// vn bkg fit function (final, after simultaneus fit) - TF1* fVnTotFunc; /// vn fit function (final, after simultaneus fit) - InvMassFitter* fMassFitter; /// mass fitter for mass prefit - Double_t fMassMin; /// upper mass limit - Double_t fMassMax; /// lower mass limit - Double_t fVn; /// vn of the signal from fit - Double_t fVnUncertainty; /// uncertainty on vn of the signal from simultaneus fit - Double_t fSigma; /// mass peak width from simultaneus fit - Double_t fSigmaUncertainty; /// uncertainty on mass peak width from simultaneus fit - Double_t fMean; /// mass peak position from simultaneus fit - Double_t fMeanUncertainty; /// uncertainty on mass peak position from simultaneus fit - Double_t fRawYield; /// raw yield from simultaneus fit - Double_t fRawYieldUncertainty; /// uncertainty raw yield from simultaneus fit - Double_t fChiSquare; /// simultaneus fit chi square - Int_t fNDF; /// simultaneus fit number of degree of freedom - Double_t fProb; /// simultaneus fit probability - Double_t fSBVnPrefitChiSquare; /// vn SB prefit chi square - Int_t fSBVnPrefitNDF; /// vn SB prefit number of degree of freedom - Double_t fSBVnPrefitProb; /// vn SB prefit probability - Double_t fMassPrefitChiSquare; /// Mass prefit chi square - Int_t fMassPrefitNDF; /// Mass prefit number of degree of freedom - Double_t fMassPrefitProb; /// Mass prefit probability - Int_t fNSigmaForSB; /// number of sigma for sidebands region (vn bkg prefit) - Double_t fSigmaInit; /// initialization for peak width - Double_t fMeanInit; /// initialization for peak position - Double_t fSigma2GausInit; /// initialization for second peak width in case of k2Gaus - Double_t fFrac2GausInit; /// initialization for fraction of second gaussian in case of k2Gaus - Bool_t fMeanFixedFromMassFit; /// flag to fix peak position from mass prefit - Bool_t fSigmaFixedFromMassFit; /// flag to fix peak width from mass prefit - Bool_t fSigma2GausFixedFromMassFit; /// flag to fix second peak width from mass prefit in case of k2Gaus - Bool_t fFrac2GausFixedFromMassFit; /// flag to fix fraction of second gaussian in case of k2Gaus - Double_t fMassParticle; /// mass of selected particle - Int_t fNParsMassSgn; /// number of parameters in mass signal fit function - Int_t fNParsMassBkg; /// number of parameters in mass bkg fit function - Int_t fNParsVnBkg; /// number of parameters in vn bkg fit function - Int_t fNParsVnSgn; /// number of parameters in vn sgn fit function (1) - Int_t fNParsVnSecPeak; /// number of parameters in vn sec peak fit function (1 if included, 0 otherwise) - Int_t fNParsVnRfl; /// number of parameters in vn refl fit function (1 if included, 0 otherwise) - Int_t fSigmaFixed; /// flag to fix peak width - Int_t fMeanFixed; /// flag to fix peak position - Int_t fSigma2GausFixed; /// flag to fix second peak width in case of k2Gaus - Int_t fFrac2GausFixed; /// flag to fix fraction of second gaussian in case of k2Gaus - Int_t fPolDegreeBkg; /// degree of polynomial expansion for back fit (option 6 for back) - Int_t fPolDegreeVnBkg; /// degree of polynomial expansion for vn back fit (option 6 for back) - Bool_t fReflections; /// flag use/not use reflections - Int_t fNParsRfl; /// fit parameters in reflection fit function - Double_t fRflOverSig; /// reflection/signal - Bool_t fFixRflOverSig; /// switch for fix refl/signal - TH1F* fHistoTemplRfl; /// histogram with reflection template - TH1F* fHistoTemplRflInit; /// initial histogram with reflection template - TF1* fMassRflFunc; /// fit function for reflections - TF1* fMassBkgRflFunc; /// mass bkg fit function plus reflections (final, after simultaneus fit) - TString fRflOpt; /// refelction option - Double_t fMinRefl; /// minimum for refelction histo - Double_t fMaxRefl; /// maximum for refelction histo - Bool_t fSmoothRfl; /// switch for smoothing of reflection template - Double_t fRawYieldHelp; /// internal variable for fit with reflections - Int_t fVnRflOpt; /// option for reflection vn type - Bool_t fVnRflLimited; /// flag to limit or not the vn of reflections - Double_t fVnRflMin; /// minimum vn of reflections - Double_t fVnRflMax; /// maximum vn of reflections - Bool_t fSecondPeak; /// switch off/on second peak (for D+->KKpi in Ds) - TF1* fMassSecPeakFunc; /// fit function for second peak - Int_t fNParsSec; /// number of parameters in second peak fit function - Double_t fSecMass; /// position of the 2nd peak - Double_t fSecWidth; /// width of the 2nd peak - Bool_t fFixSecMass; /// flag to fix the position of the 2nd peak - Bool_t fFixSecWidth; /// flag to fix the width of the 2nd peak - Double_t fVnSecPeak; /// vn of second peak from fit - Bool_t fDoSecondPeakVn; /// flag to introduce second peak vn in the vn vs. mass fit - Double_t fVnSecPeakUncertainty; /// vn uncertainty of second peak from fit - Int_t fHarmonic; /// harmonic number for drawing + TH1F* fMassHisto; /// mass histogram to fit + TH1F* fVnVsMassHisto; /// vn vs. mass histogram to fit + Int_t fMassSgnFuncType; /// type of mass signal fit function + Int_t fMassBkgFuncType; /// type of mass bkg fit function + Int_t fVnBkgFuncType; /// type of vn bkg fit function + TF1* fMassFuncFromPrefit; /// mass fit function (1st step, from prefit) + TF1* fMassBkgFunc; /// mass bkg fit function (final, after simultaneus fit) + TF1* fMassSgnFunc; /// mass signal fit function (final, after simultaneus fit) + TF1* fMassTotFunc; /// mass fit function (final, after simultaneus fit) + TF1* fVnBkgFuncSb; /// vn bkg fit function (1st step from SB prefit) + TF1* fVnBkgFunc; /// vn bkg fit function (final, after simultaneus fit) + TF1* fVnTotFunc; /// vn fit function (final, after simultaneus fit) + InvMassFitter* fMassFitter; /// mass fitter for mass prefit + Double_t fMassMin; /// upper mass limit + Double_t fMassMax; /// lower mass limit + Double_t fVn; /// vn of the signal from fit + Double_t fVnUncertainty; /// uncertainty on vn of the signal from simultaneus fit + Double_t fSigma; /// mass peak width from simultaneus fit + Double_t fSigmaUncertainty; /// uncertainty on mass peak width from simultaneus fit + Double_t fMean; /// mass peak position from simultaneus fit + Double_t fMeanUncertainty; /// uncertainty on mass peak position from simultaneus fit + Double_t fRawYield; /// raw yield from simultaneus fit + Double_t fRawYieldUncertainty; /// uncertainty raw yield from simultaneus fit + Double_t fChiSquare; /// simultaneus fit chi square + Int_t fNDF; /// simultaneus fit number of degree of freedom + Double_t fProb; /// simultaneus fit probability + Double_t fSBVnPrefitChiSquare; /// vn SB prefit chi square + Int_t fSBVnPrefitNDF; /// vn SB prefit number of degree of freedom + Double_t fSBVnPrefitProb; /// vn SB prefit probability + Double_t fMassPrefitChiSquare; /// Mass prefit chi square + Int_t fMassPrefitNDF; /// Mass prefit number of degree of freedom + Double_t fMassPrefitProb; /// Mass prefit probability + Int_t fNSigmaForSB; /// number of sigma for sidebands region (vn bkg prefit) + Double_t fSigmaInit; /// initialization for peak width + Double_t fMeanInit; /// initialization for peak position + Double_t fSigma2GausInit; /// initialization for second peak width in case of k2Gaus + Double_t fFrac2GausInit; /// initialization for fraction of second gaussian in case of k2Gaus + Bool_t fMeanFixedFromMassFit; /// flag to fix peak position from mass prefit + Bool_t fSigmaFixedFromMassFit; /// flag to fix peak width from mass prefit + Bool_t fSigma2GausFixedFromMassFit; /// flag to fix second peak width from mass prefit in case of k2Gaus + Bool_t fFrac2GausFixedFromMassFit; /// flag to fix fraction of second gaussian in case of k2Gaus + Double_t fMassParticle; /// mass of selected particle + Int_t fNParsMassSgn; /// number of parameters in mass signal fit function + Int_t fNParsMassBkg; /// number of parameters in mass bkg fit function + Int_t fNParsVnBkg; /// number of parameters in vn bkg fit function + Int_t fNParsVnSgn; /// number of parameters in vn sgn fit function (1) + Int_t fNParsVnSecPeak; /// number of parameters in vn sec peak fit function (1 if included, 0 otherwise) + Int_t fNParsVnRfl; /// number of parameters in vn refl fit function (1 if included, 0 otherwise) + Int_t fSigmaFixed; /// flag to fix peak width + Int_t fMeanFixed; /// flag to fix peak position + Int_t fSigma2GausFixed; /// flag to fix second peak width in case of k2Gaus + Int_t fFrac2GausFixed; /// flag to fix fraction of second gaussian in case of k2Gaus + Int_t fPolDegreeBkg; /// degree of polynomial expansion for back fit (option 6 for back) + Int_t fPolDegreeVnBkg; /// degree of polynomial expansion for vn back fit (option 6 for back) + Bool_t fReflections; /// flag use/not use reflections + Int_t fNParsRfl; /// fit parameters in reflection fit function + Double_t fRflOverSig; /// reflection/signal + Bool_t fFixRflOverSig; /// switch for fix refl/signal + TH1F* fHistoTemplRfl; /// histogram with reflection template + TH1F* fHistoTemplRflInit; /// initial histogram with reflection template + TF1* fMassRflFunc; /// fit function for reflections + TF1* fMassBkgRflFunc; /// mass bkg fit function plus reflections (final, after simultaneus fit) + TString fRflOpt; /// refelction option + Double_t fMinRefl; /// minimum for refelction histo + Double_t fMaxRefl; /// maximum for refelction histo + Bool_t fSmoothRfl; /// switch for smoothing of reflection template + Double_t fRawYieldHelp; /// internal variable for fit with reflections + Int_t fVnRflOpt; /// option for reflection vn type + Bool_t fVnRflLimited; /// flag to limit or not the vn of reflections + Double_t fVnRflMin; /// minimum vn of reflections + Double_t fVnRflMax; /// maximum vn of reflections + Bool_t fSecondPeak; /// switch off/on second peak (for D+->KKpi in Ds) + TF1* fMassSecPeakFunc; /// fit function for second peak + Int_t fNParsSec; /// number of parameters in second peak fit function + Double_t fSecMass; /// position of the 2nd peak + Double_t fSecWidth; /// width of the 2nd peak + Bool_t fFixSecMass; /// flag to fix the position of the 2nd peak + Bool_t fFixSecWidth; /// flag to fix the width of the 2nd peak + Double_t fVnSecPeak; /// vn of second peak from fit + Bool_t fDoSecondPeakVn; /// flag to introduce second peak vn in the vn vs. mass fit + Double_t fVnSecPeakUncertainty; /// vn uncertainty of second peak from fit + Int_t fHarmonic; /// harmonic number for drawing /// \cond CLASSDEF ClassDef(VnVsMassFitter,5); From 3a779173376aa711fff00d68f7f39ec33bd104ea Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Wed, 13 Nov 2024 19:13:25 +0100 Subject: [PATCH 3/6] [WIP] implemented fitting with kde --- run3/flow/flow_analysis_utils.py | 2 + run3/flow/get_vn_vs_mass.py | 56 ++++++- run3/flow/invmassfitter/InvMassFitter.cxx | 83 +++++++++- run3/flow/invmassfitter/InvMassFitter.h | 19 ++- run3/flow/invmassfitter/VnVsMassFitter.cxx | 179 +++++++++++++++++++-- run3/flow/invmassfitter/VnVsMassFitter.h | 32 ++++ run3/flow/kde_producer.py | 70 +++++--- 7 files changed, 398 insertions(+), 43 deletions(-) diff --git a/run3/flow/flow_analysis_utils.py b/run3/flow/flow_analysis_utils.py index b451f6f3..c780a98c 100644 --- a/run3/flow/flow_analysis_utils.py +++ b/run3/flow/flow_analysis_utils.py @@ -417,6 +417,8 @@ def get_vnfitter_results(vnFitter, secPeak, useRefl): vn_results['fBkgFuncMass'] = vnFitter.GetMassBkgFitFunc() vn_results['fBkgFuncVn'] = vnFitter.GetVnVsMassBkgFitFunc() vn_results['fSgnFuncMass'] = vnFitter.GetMassSignalFitFunc() + vn_results['fMassTemplFuncts'] = vnFitter.GetMassTemplFuncts() + vn_results['fVnTemplFuncts'] = vnFitter.GetVnTemplFuncts() bkg, bkgUnc = ctypes.c_double(), ctypes.c_double() vnFitter.Background(3, bkg, bkgUnc) vn_results['bkg'] = bkg.value diff --git a/run3/flow/get_vn_vs_mass.py b/run3/flow/get_vn_vs_mass.py index 25d91fb9..044738e3 100644 --- a/run3/flow/get_vn_vs_mass.py +++ b/run3/flow/get_vn_vs_mass.py @@ -19,6 +19,7 @@ from ROOT import InvMassFitter, VnVsMassFitter from utils.StyleFormatter import SetGlobalStyle, SetObjectStyle, DivideCanvas from utils.FitUtils import SingleGaus, DoubleGaus, DoublePeakSingleGaus, DoublePeakDoubleGaus, RebinHisto +from kde_producer import kde_producer, kde_producer_sim def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, outputdir, suffix, vn_method, batch): @@ -119,7 +120,31 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, else: print('ERROR: only kGaus, k2Gaus and k2GausSigmaRatioPar signal functions supported! Exit!') sys.exit() - + + useKDEtemplates = True if fitConfig.get('UseKDETemplates') else False + KDEtemplates = [[None]*len(fitConfig['MCTemplatesFlags'][iPt]) for iPt in range(len(ptMins))] + KDEtemplatesNames = [[None]*len(fitConfig['MCTemplatesFlags'][iPt]) for iPt in range(len(ptMins))] + if useKDEtemplates: + for iPt in range(len(ptMins)): + for iFlag, flag in enumerate(fitConfig['MCTemplatesFlags'][iPt]): + if fitConfig.get('ProduceNew') and not fitConfig.get('FromSim'): + KDEtemplates[iPt][iFlag], hist_templ = kde_producer(fitConfig['MCTemplatesInputs'][iPt][iFlag], + 'fM', ptMins[iPt], ptMaxs[iPt], + fitConfig['inv_mass_bins'][iPt], flag) + KDEtemplatesNames[iPt][iFlag] = fitConfig['MCTemplatesNames'][iPt][iFlag] + elif fitConfig.get('ProduceNew') and fitConfig.get('FromSim'): + KDEtemplates[iPt][iFlag] = kde_producer_sim(fitConfig['MCTemplatesInputs'][iPt][iFlag], + fitConfig['MCTemplatesTreeNames'][iPt][iFlag], + ptMins[iPt], ptMaxs[iPt]) + KDEtemplatesNames[iPt][iFlag] = fitConfig['MCTemplatesTreeNames'][iPt][iFlag] + else: + templFile = TFile.Open(f'{fitConfig["MCTemplatesFiles"][iPt][iFlag]}', 'r') + KDEtemplates[iPt][iFlag] = templFile.Get(f'KDE_pt_{ptMins[iPt]}_{ptMaxs[iPt]}_flag{flag}/kde_') + print(f"Number of KDE templates loaded: {len(KDEtemplates)}") + print(KDEtemplates) + KDEtemplatesFuncts = [[KDE.GetFunction() for KDE in KDEtemplatesPt] for KDEtemplatesPt in KDEtemplates] + print(KDEtemplatesFuncts) + # set particle configuration if particleName == 'Ds': _, massAxisTit, decay, massForFit = get_particle_info(particleName) @@ -137,6 +162,8 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, hMassIns, hMassOuts, hMassInsForFit, hMassOutsForFit = [], [], [], [] fTotFuncMass, fTotFuncVn, fSgnFuncMass, fBkgFuncMass, fMassBkgRflFunc,fBkgFuncVn = [], [], [], [], [], [] hMCSgn, hMCRefl = [], [] + fMassTemplFuncts = [[None]*len(fitConfig['MCTemplatesFlags'][iPt]) for iPt in range(len(ptMins))] + fVnTemplFuncts = [[None]*len(fitConfig['MCTemplatesFlags'][iPt]) for iPt in range(len(ptMins))] hist_reso = infile.Get('hist_reso') hist_reso.SetDirectory(0) reso = hist_reso.GetBinContent(1) @@ -366,6 +393,11 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, vnFitter[iPt].SetTemplateReflections(hMCRefl[iPt],reflFuncStr,massMin,massMax) vnFitter[iPt].SetFixReflOverS(SoverR) vnFitter[iPt].SetReflVnOption(0) # kSameVnSignal + if useKDEtemplates: + vnFitter[iPt].SetKDETemplates(KDEtemplatesFuncts[iPt], fitConfig['MCTemplatesTreeNames'][iPt], + fitConfig['InitWeights'][iPt], fitConfig['MinWeights'][iPt], + fitConfig['MaxWeights'][iPt]) + # collect fit results vnFitter[iPt].SimultaneousFit(False) vnResults = get_vnfitter_results(vnFitter[iPt], secPeak, useRefl) @@ -374,6 +406,8 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, fSgnFuncMass.append(vnResults['fSgnFuncMass']) fBkgFuncMass.append(vnResults['fBkgFuncMass']) fBkgFuncVn.append(vnResults['fBkgFuncVn']) + fMassTemplFuncts[iPt] = vnResults['fMassTemplFuncts'] + fVnTemplFuncts[iPt] = vnResults['fVnTemplFuncts'] if useRefl: fMassBkgRflFunc.append(vnResults['fMassBkgRflFunc']) @@ -442,6 +476,14 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, latex.DrawLatex(0.18, 0.60, f'Signif. (3#sigma) = {round(vnResults["signif"], 2)}') if useRefl: latex.DrawLatex(0.18, 0.20, f'RoverS = {SoverR:.2f}') + print(f'Drawing fMassTemplFuncts') + print(fMassTemplFuncts) + for iMassTemplFunct in fMassTemplFuncts[iPt]: + print(f'Eval: {iMassTemplFunct.Eval(1.8)}') + iMassTemplFunct.SetLineColor(1) + iMassTemplFunct.Draw('same') + cSimFit[iCanv].Modified() + cSimFit[iCanv].Update() cSimFit[iPt].cd(2) hVnForFit[iPt].GetYaxis().SetRangeUser(-1, 1) hVnForFit[iPt].GetYaxis().SetTitle(f'#it{{v}}_{{{harmonic}}} ({vn_method})') @@ -458,6 +500,13 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, if secPeak: latex.DrawLatex(0.18, 0.75, f'#it{{v}}{harmonic}(D^{{+}}) = {vnResults["vnSecPeak"]:.3f} #pm {vnResults["vnSecPeakUnc"]:.3f}') + print(f'Drawing fVnTemplFuncts') + print(fVnTemplFuncts) + for iVnTemplFunct in fVnTemplFuncts[iPt]: + print(f'Eval: {iVnTemplFunct.Eval(1.8)}') + iVnTemplFunct.Draw('same') + cSimFit[iCanv].Modified() + cSimFit[iCanv].Update() cSimFit[iCanv].Modified() cSimFit[iCanv].Update() @@ -717,6 +766,11 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, fSgnFuncMass[ipt].Write(f'fSgnFuncMass_pt{ptmin*10:.0f}_{ptmax*10:.0f}') fBkgFuncMass[ipt].Write(f'fBkgFuncMass_pt{ptmin*10:.0f}_{ptmax*10:.0f}') fBkgFuncVn[ipt].Write(f'fBkgFuncVn_pt{ptmin*10:.0f}_{ptmax*10:.0f}') + print(f'Writing templates for pt{iPt}') + print(KDEtemplatesFuncts) + for iFlag in range(len(KDEtemplatesFuncts[ipt])): + KDEtemplatesFuncts[ipt][iFlag].Write(f'Templ_pt{ptmin*10:.0f}_{ptmax*10:.0f}_flag{fitConfig["MCTemplatesFlags"][iPt][iFlag]}') + hSigmaSimFit.Write() hMeanSimFit.Write() hMeanSecPeakFitMass.Write() diff --git a/run3/flow/invmassfitter/InvMassFitter.cxx b/run3/flow/invmassfitter/InvMassFitter.cxx index 18e607f1..cf63bfed 100644 --- a/run3/flow/invmassfitter/InvMassFitter.cxx +++ b/run3/flow/invmassfitter/InvMassFitter.cxx @@ -93,6 +93,8 @@ InvMassFitter::InvMassFitter() : fFixRflOverSig(kFALSE), fHistoTemplRfl(0x0), fSmoothRfl(kFALSE), + fTemplates(kFALSE), + fNParsTempls(0), fRawYieldHelp(0), fRflFunc(0x0), fBkRFunc(0x0), @@ -157,6 +159,8 @@ InvMassFitter::InvMassFitter(const TH1F *histoToFit, Double_t minvalue, Double_t fFixRflOverSig(kFALSE), fHistoTemplRfl(0x0), fSmoothRfl(kFALSE), + fTemplates(kFALSE), + fNParsTempls(0), fRawYieldHelp(0), fRflFunc(0x0), fBkRFunc(0x0), @@ -242,6 +246,9 @@ void InvMassFitter::SetNumberOfParams(){ if(fReflections) fNParsRfl=1; else fNParsRfl=0; + if(fTemplates) fNParsTempls=this->fTemplatesFuncts.size(); + else fNParsTempls=0; + if(fSecondPeak) fNParsSec=3; else fNParsSec=0; @@ -252,9 +259,16 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){ /// returns 0 if the fit fails /// returns 1 if the fit succeeds /// returns 2 if there is no signal and the fit is performed with only background - + cout << "Entering InvMassFitter::MassFitter" << endl; TVirtualFitter::SetDefaultFitter("Minuit"); + // for(int iFunc=0; iFuncIntegral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width"); fOnlySideBands = kTRUE; @@ -314,8 +328,21 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){ fRflFunc = CreateReflectionFunction("freflect"); fBkRFunc = CreateBackgroundPlusReflectionFunction("fbkgrfl"); } + // cout << "sfTemplates: " << fTemplates << endl; + if(fTemplates){ + printf(" ---> Final fit includes templates\n"); + fTemplFunc = CreateTemplatesFunction("ftempl"); + } + // for(int iFunc=0; iFunccd(); + fTotFunc->Write(); + funcFile->Close(); + if(doFinalFit){ printf("\n--- Final fit with signal+background on the full range ---\n"); TFitResultPtr resultptr=fHistoInvMass->Fit("funcmass",Form("R,S,%s,+,0",fFitOption.Data())); @@ -360,6 +387,13 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){ fBkRFunc->SetLineColor(kRed+1); fBkRFunc->SetLineStyle(7); } + if(fTemplates){ + for(Int_t ipar=0; iparSetParameter(ipar,fTotFunc->GetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl)); + fTemplFunc->SetParError(ipar,fTotFunc->GetParError(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl)); + } + fTemplFunc->SetLineColor(kGreen+1); + } fMass=fSigFunc->GetParameter(1); fMassErr=fSigFunc->GetParError(1); fSigmaSgn=fSigFunc->GetParameter(2); @@ -473,6 +507,21 @@ TF1* InvMassFitter::CreateReflectionFunction(TString fname){ return funcrfl; } //______________________________________________________________________________ +TF1* InvMassFitter::CreateTemplatesFunction(TString fname){ + /// Creates a function for templates in the D+ inv. mass distribution + TF1* functempl = new TF1(fname.Data(),this,&InvMassFitter::FitFunction4Templ,fMinMass,fMaxMass,this->fTemplatesFuncts.size(),"InvMassFitter","FitFunction4Templ"); + for(int iPar=0; iParGetNpar(); iPar++) { + functempl->SetParName(iPar, Form("w_%s",this->fTemplatesFuncts[iPar].GetName())); + if(this->fWeightsLowerLims[iPar] >= this->fWeightsUpperLims[iPar]) { + functempl->FixParameter(iPar,this->fInitWeights[iPar]); + } else { + functempl->SetParameter(iPar,this->fInitWeights[iPar]); + functempl->SetParLimits(iPar,this->fWeightsLowerLims[iPar],this->fWeightsUpperLims[iPar]); + } + } + return functempl; +} +//______________________________________________________________________________ TF1* InvMassFitter::CreateBackgroundPlusReflectionFunction(TString fname){ /// Creates the function with sum of background and reflections /// @@ -549,11 +598,11 @@ TF1* InvMassFitter::CreateSignalFitFunction(TString fname, Double_t integsig){ //______________________________________________________________________________ TF1* InvMassFitter::CreateTotalFitFunction(TString fname){ - /// Creates the total fit fucntion (signal+background+possible second peak) + /// Creates the total fit function (signal+background+possible second peak) /// SetNumberOfParams(); - Int_t totParams=fNParsBkg+fNParsRfl+fNParsSec+fNParsSig; + Int_t totParams=fNParsBkg+fNParsRfl+fNParsTempls+fNParsSec+fNParsSig; TF1* ftot=new TF1(fname.Data(),this,&InvMassFitter::FitFunction4Mass,fMinMass,fMaxMass,totParams,"InvMassFitter","FitFunction4Mass"); for(Int_t ipar=0; iparSetParameter(ipar,fBkgFunc->GetParameter(ipar)); @@ -587,6 +636,22 @@ TF1* InvMassFitter::CreateTotalFitFunction(TString fname){ ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec,parmin,parmax); } } + if(fTemplates && fTemplFunc){ + cout << "Setting template parameters" << endl; + cout << "Number of template parameters " << fNParsTempls << endl; + cout << "Number of parameters of fTemplFunc " << fTemplFunc->GetNpar() << endl; + for(Int_t ipar=0; iparGetParLimits(ipar,parmin,parmax); + ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl,parmin,parmax); + cout << "Setting iPar" << ipar << " to " << fTemplFunc->GetParameter(ipar) + << ", parname: " << fTemplFunc->GetParName(ipar) << ", parlimits: [" + << parmin << "," << parmax << "]" << endl; + ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl,fTemplFunc->GetParameter(ipar)); + ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl,fTemplFunc->GetParName(ipar)); + } + } + cout << "Returning fTot" << endl; return ftot; } //__________________________________________________________________________ @@ -773,6 +838,14 @@ Double_t InvMassFitter::FitFunction4SecPeak (Double_t *x, Double_t *par){ return secgaval; } //_________________________________________________________________________ +Double_t InvMassFitter::FitFunction4Templ(Double_t *x, Double_t *par){ + Double_t totalTempl=0; + for(int iFunc=0; iFuncfTemplatesFuncts.size(); iFunc++) { + totalTempl += par[iFunc] * this->fTemplatesFuncts[iFunc].Eval(x[0]); + } + return totalTempl; +} +//_________________________________________________________________________ Double_t InvMassFitter::FitFunction4Mass(Double_t *x, Double_t *par){ /// Total fit function (signal+background+possible second peak) /// @@ -783,7 +856,9 @@ Double_t InvMassFitter::FitFunction4Mass(Double_t *x, Double_t *par){ if(fSecondPeak) sec=FitFunction4SecPeak(x,&par[fNParsBkg+fNParsSig]); Double_t refl=0; if(fReflections) refl=FitFunction4Refl(x,&par[fNParsBkg+fNParsSig+fNParsSec]); - return bkg+sig+sec+refl; + Double_t templ=0; + if(fTemplates) templ=FitFunction4Templ(x,&par[fNParsBkg+fNParsSig+fNParsSec+fNParsRfl]); + return bkg+sig+sec+refl+templ; } //_________________________________________________________________________ diff --git a/run3/flow/invmassfitter/InvMassFitter.h b/run3/flow/invmassfitter/InvMassFitter.h index 58615fd1..84b952fd 100644 --- a/run3/flow/invmassfitter/InvMassFitter.h +++ b/run3/flow/invmassfitter/InvMassFitter.h @@ -99,6 +99,14 @@ class InvMassFitter : public TNamed { fFixRflOverSig=kTRUE; } void SetSmoothReflectionTemplate(Bool_t opt){fSmoothRfl=opt;} + void SetTemplates(std::vector templates, std::vector initweights, + std::vector minweights, std::vector maxweights) { + fTemplatesFuncts=templates; + fInitWeights=initweights; + fWeightsLowerLims=minweights; + fWeightsUpperLims=maxweights; + fTemplates=kTRUE; + } void IncludeSecondGausPeak(Double_t mass, Bool_t fixm, Double_t width, Bool_t fixw){ fSecondPeak=kTRUE; fSecMass=mass; fSecWidth=width; @@ -151,6 +159,7 @@ class InvMassFitter : public TNamed { Double_t FitFunction4Sgn (Double_t* x, Double_t* par); Double_t FitFunction4Bkg (Double_t* x, Double_t* par); Double_t FitFunction4Refl(Double_t *x,Double_t *par); + Double_t FitFunction4Templ(Double_t *x,Double_t *par); Double_t FitFunction4BkgAndRefl(Double_t *x,Double_t *par); Double_t FitFunction4SecPeak (Double_t* x, Double_t* par); Double_t FitFunction4Mass (Double_t* x, Double_t* par); @@ -182,6 +191,7 @@ class InvMassFitter : public TNamed { TF1* CreateSignalFitFunction(TString fname, Double_t integral); TF1* CreateSecondPeakFunction(TString fname, Double_t integral); TF1* CreateReflectionFunction(TString fname); + TF1* CreateTemplatesFunction(TString fname); TF1* CreateBackgroundPlusReflectionFunction(TString fname); TF1* CreateTotalFitFunction(TString fname); Bool_t PrepareHighPolFit(TF1 *fback); @@ -244,7 +254,14 @@ class InvMassFitter : public TNamed { TF1* fSecFunc; /// fit function for second peak TF1* fTotFunc; /// total fit function Bool_t fAcceptValidFit; /// accept a fit when IsValid() gives true, nevertheless the status code - + std::vector fTemplatesFuncts; /// vector storing TF1 to be added as templates + TF1* fTemplFunc; /// fit function for templates + Bool_t fTemplates; /// flag use/not use templates fit functions + Int_t fNParsTempls; /// fit parameters in templates fit function + std::vector fWeightsUpperLims; /// upper limit of the templates' weights + std::vector fWeightsLowerLims; /// lower limit of the templates' weights + std::vector fInitWeights; /// init value of the templates' weights + /// \cond CLASSIMP ClassDef(InvMassFitter,9); /// class for invariant mass fit /// \endcond diff --git a/run3/flow/invmassfitter/VnVsMassFitter.cxx b/run3/flow/invmassfitter/VnVsMassFitter.cxx index 1fe15cbd..ab269fc9 100755 --- a/run3/flow/invmassfitter/VnVsMassFitter.cxx +++ b/run3/flow/invmassfitter/VnVsMassFitter.cxx @@ -79,6 +79,8 @@ VnVsMassFitter::VnVsMassFitter() ,fPolDegreeVnBkg(3) ,fReflections(kFALSE) ,fNParsRfl(0) + ,fTemplates(kFALSE) + ,fNParsTempls(0) ,fRflOverSig(0.) ,fFixRflOverSig(kFALSE) ,fHistoTemplRfl(0x0) @@ -164,6 +166,8 @@ VnVsMassFitter::VnVsMassFitter(TH1F* hMass, TH1F* hvn, Double_t min, Double_t ma ,fPolDegreeVnBkg(3) ,fReflections(kFALSE) ,fNParsRfl(0) + ,fTemplates(kFALSE) + ,fNParsTempls(0) ,fRflOverSig(0.) ,fFixRflOverSig(kFALSE) ,fHistoTemplRfl(0x0) @@ -221,33 +225,55 @@ VnVsMassFitter::~VnVsMassFitter() { //________________________________________________________________ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { + cout << "PERFORMING SIMULTANEOUS FIT" << endl; if(!fMassHisto || !fVnVsMassHisto) {printf("Histograms not set! Exit."); return kFALSE;} DefineNumberOfParameters(); - const Int_t nparsmass = fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl; + const Int_t nparsmass = fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls; + cout << "nParsMass " << nparsmass << endl; Int_t NvnParsSgn = 1; if(fSecondPeak && fDoSecondPeakVn) {NvnParsSgn+=1;} if(fReflections && fVnRflOpt==kFreePar) {NvnParsSgn+=1;} + if(fTemplates) {NvnParsSgn+=fNParsTempls;} + cout << "Added " << fNParsTempls << " to vn sgn pars" << endl; const Int_t nparsvn = nparsmass+fNParsVnBkg+NvnParsSgn; + cout << "nparsvn " << nparsvn << endl; + cout << "Starting PREFITS" << endl; Bool_t massprefit=MassPrefit(); if(!massprefit) {printf("Impossible to perform the mass prefit"); return kFALSE;} + cout << "Mass PREFIT" << endl; Bool_t vnprefit=VnSBPrefit(); + cout << "vnprefit: " << vnprefit << endl; if(!vnprefit) {printf("Impossible to perform the bkg vn prefit"); return kFALSE;} + cout << "Vn PREFIT" << endl; + cout << "PREFITS DONE" << endl; + cout << "Mass function obtained from the prefit" << endl; + cout << "Npar: " << fMassFuncFromPrefit->GetNpar() << endl; std::vector initpars; for(Int_t iBkgPar=0; iBkgParGetParName(iBkgPar) << ", value: " << fMassFuncFromPrefit->GetParameter(iBkgPar) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iBkgPar)); } for(Int_t iSgnPar=0; iSgnParGetParName(iSgnPar+fNParsMassBkg) << ", value: " << fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg)); } for(Int_t iSecPeakPar=0; iSecPeakParGetParName(iSecPeakPar+fNParsMassBkg+fNParsMassSgn) << ", value: " << fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn)); } for(Int_t iReflPar=0; iReflParGetParName(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec) << ", value: " << fMassFuncFromPrefit->GetParameter(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec)); } + for(Int_t iTemplPar=0; iTemplParGetParName(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl) << ", value: " << fMassFuncFromPrefit->GetParameter(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl) << endl; + initpars.push_back(fMassFuncFromPrefit->GetParameter(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl)); + } + cout << "------------" << endl; + cout << "Background parameters for the VnBkgFunc: " << fNParsVnBkg << endl; for(Int_t iVnBkgPar=0; iVnBkgParGetParameter(iVnBkgPar));} else {initpars.push_back(0.05);} @@ -256,9 +282,16 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { if(fSecondPeak && fDoSecondPeakVn) {initpars.push_back(0.10);} //initial parameter for second peak vn fMassTotFunc = new TF1("fMassTotFunc",this,&VnVsMassFitter::MassFunc,fMassMin,fMassMax,nparsmass,"VnVsMassFitter","MassFunc"); + cout << "Number of parameters for vn function: " << nparsvn << endl; fVnTotFunc = new TF1("fVnTotFunc",this,&VnVsMassFitter::vnFunc,fMassMin,fMassMax,nparsvn,"VnVsMassFitter","vnFunc"); SetParNames(); + for(int iMassPar=0; iMassParGetNpar(); iMassPar++) { + cout << "[Mass] iPar" << iMassPar << ", name: " << fMassTotFunc->GetParName(iMassPar) << endl; + } + for(int iVnPar=0; iVnParGetNpar(); iVnPar++) { + cout << "[Vn] iPar" << iVnPar << ", name: " << fVnTotFunc->GetParName(iVnPar) << endl; + } ROOT::Math::WrappedMultiTF1 wfTotMass(*fMassTotFunc,1); ROOT::Math::WrappedMultiTF1 wfTotVn(*fVnTotFunc,1); @@ -297,13 +330,49 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { if(fFixRflOverSig) fitter.Config().ParSettings(fNParsMassBkg+fNParsMassSgn+fNParsSec).Fix(); if(fVnRflLimited) fitter.Config().ParSettings(nparsmass+fNParsVnBkg+NvnParsSgn-1).SetLimits(fVnRflMin,fVnRflMax); } + if(fTemplates) { + for(int iTemplPar=0; iTemplParfInitWeights.size()/2; iTemplPar++) { + cout << "Setting template parameter " << iTemplPar << endl; + cout << "Mass func parameter " << iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl << endl; + cout << fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Name() << endl; + cout << "[" << fInitWeights[iTemplPar] << ", " << fWeightsLowerLims[iTemplPar] + << ", " << fWeightsUpperLims[iTemplPar] << "]" << endl; + cout << "Vn func parameter " << iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn << endl; + cout << "[" << fInitWeights[iTemplPar+this->fInitWeights.size()/2] << ", " << fWeightsLowerLims[iTemplPar+this->fInitWeights.size()/2] + << ", " << fWeightsUpperLims[iTemplPar+this->fInitWeights.size()/2] << "]" << endl; + cout << fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).Name() << endl; + + cout << "---" << endl; + if(this->fWeightsLowerLims[iTemplPar] > this->fWeightsUpperLims[iTemplPar]) { + cout << "Fix mass" << endl; + fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); + fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Fix(); + } else { + cout << "Vary mass" << endl; + fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); + fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetLimits(fWeightsLowerLims[iTemplPar],fWeightsUpperLims[iTemplPar]); + } + if(this->fWeightsLowerLims[iTemplPar] > this->fWeightsUpperLims[iTemplPar]) { + cout << "Fix vn" << endl; + fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).SetValue(fInitWeights[iTemplPar+this->fInitWeights.size()/2]); + fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).Fix(); + } else { + cout << "Vary vn" << endl; + fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).SetValue(fInitWeights[iTemplPar+this->fInitWeights.size()/2]); + fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).SetLimits(fWeightsLowerLims[iTemplPar+this->fInitWeights.size()/2],fWeightsUpperLims[iTemplPar+this->fInitWeights.size()/2]); + } + cout << "---" << endl; + } + } fitter.Config().MinimizerOptions().SetPrintLevel(0); fitter.Config().SetMinimizer("Minuit2","Migrad"); for(Int_t iPar=0; iParGetParName(iPar));} + cout << "nparsvn: " << nparsvn << endl; // fit FCN function directly // (specify optionally data size and flag to indicate that is a chi2 fit Bool_t isFitOk = fitter.FitFCN(nparsvn,globalChi2,0,dataMass.Size()+dataVn.Size(),kFALSE); + cout << "isFitOk: " << isFitOk << endl; if(!isFitOk) return kFALSE; ROOT::Fit::FitResult result = fitter.Result(); @@ -349,6 +418,32 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { fMassSecPeakFunc->SetParError(iPar-(fNParsMassBkg+fNParsMassSgn),result.ParError(iPar)); } } + if(fTemplates) { + int nParBeforeTemplsMass = fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl; + cout << "Coeff first func mass: " << result.Parameter(nParBeforeTemplsMass) << endl; + cout << "Eval func: " << this->fKDETemplates[0].Eval(1.8) << endl; + for(int iTempl=0; iTemplfKDETemplates[iTempl].Eval(x[0]); + return paramValue * kdeEval; + }, fMassMin, fMassMax, 0)); + } + int nParBeforeTemplsVn = nParBeforeTemplsMass+fNParsTempls+fNParsVnBkg+fNParsVnSgn; + cout << "Coeff first func vn: " << result.Parameter(nParBeforeTemplsVn) << endl; + cout << "Eval func: " << this->fKDETemplates[fKDETemplates.size() / 2].Eval(1.8) << endl; + for(int iTempl=0; iTemplfKDETemplates[iTempl + fKDETemplates.size() / 2].Eval(x[0]); + return paramValue * kdeEval; + }, fMassMin, fMassMax, 0)); + } + } + + cout << "DrawFit: " << drawFit << endl; if(drawFit) {DrawFit();} fVn = fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-NvnParsSgn); @@ -366,7 +461,7 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { fChiSquare = result.MinFcnValue(); fNDF = result.Ndf(); fProb = result.Prob(); - + cout << "Combined fit completed" << endl; return kTRUE; } @@ -374,6 +469,7 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { void VnVsMassFitter::DrawHere(TVirtualPad* c){ /// Core method to draw the fit output + cout << "[VnVsMassFitter] DrawHere" << endl; gStyle->SetOptStat(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); @@ -421,6 +517,15 @@ void VnVsMassFitter::DrawHere(TVirtualPad* c){ fMassTotFunc->SetRange(fMassMin,fMassMax); fMassTotFunc->Draw("same"); } + if(fTemplates) { + cout << "Drawing templates for mass spectrum" << endl; + for(int iMassTempl=0; iMassTemplfKDETemplates.size()/2; iMassTempl++) { + fKDEMassTemplatesDraw[iMassTempl]->SetLineColor(kMagenta); + fKDEMassTemplatesDraw[iMassTempl]->SetRange(fMassMin,fMassMax); + fKDEMassTemplatesDraw[iMassTempl]->Draw("same"); + } + } + TPaveText* massinfo = new TPaveText(0.45,0.7,1.,0.87,"NDC"); massinfo->SetTextFont(42); massinfo->SetTextSize(0.05); @@ -459,6 +564,14 @@ void VnVsMassFitter::DrawHere(TVirtualPad* c){ fVnTotFunc->SetRange(fMassMin,fMassMax); fVnTotFunc->Draw("same"); } + if(fTemplates) { + cout << "Drawing templates for vn" << endl; + for(int iVnTempl=this->fKDETemplates.size()/2; iVnTemplfKDETemplates.size(); iVnTempl++) { + fKDEVnTemplatesDraw[iVnTempl]->SetLineColor(kMagenta); + fKDEVnTemplatesDraw[iVnTempl]->SetRange(fMassMin,fMassMax); + fKDEVnTemplatesDraw[iVnTempl]->Draw("same"); + } + } TPaveText* vninfo = new TPaveText(-0.45,0.7,1.,0.87,"NDC"); vninfo->SetTextFont(42); @@ -504,6 +617,7 @@ Bool_t VnVsMassFitter::MassPrefit() { if(fRflOverSig>0) {fMassFitter->SetInitialReflOverS(fRflOverSig);} if(fFixRflOverSig) {fMassFitter->SetFixReflOverS(fRflOverSig);} } + if(fTemplates) {fMassFitter->SetTemplates(fKDETemplates, fInitWeights, fWeightsLowerLims, fWeightsUpperLims);} Bool_t status = fMassFitter->MassFitter(kFALSE); if(status) { @@ -542,6 +656,7 @@ Bool_t VnVsMassFitter::VnSBPrefit() { TGraphErrors* gVnVsMassSB = new TGraphErrors(nSBbins); for(Int_t iBin=0; iBinGetBinContent(iBin+1) << endl; gVnVsMassSB->SetPoint(iBin,fVnVsMassHisto->GetBinCenter(iBin+1),fVnVsMassHisto->GetBinContent(iBin+1)); gVnVsMassSB->SetPointError(iBin,fVnVsMassHisto->GetBinWidth(iBin+1)/2,fVnVsMassHisto->GetBinError(iBin+1)); } @@ -644,6 +759,13 @@ void VnVsMassFitter::DefineNumberOfParameters() { fNParsRfl=0; fNParsVnRfl=0; } + + if(fTemplates) { + fNParsTempls=fKDETemplates.size(); + } + else { + fNParsTempls=0; + } if(fSecondPeak) { fNParsSec=3; @@ -710,19 +832,27 @@ void VnVsMassFitter::SetParNames() { break; } - for(Int_t iPar=0; iParSetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl+iPar,fVnBkgFuncSb->GetParName(iPar));} + for(Int_t iPar=0; iParSetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl+fNParsTempls+iPar,fVnBkgFuncSb->GetParName(iPar));} if(fReflections) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec,"ReflOverS");} + if(fTemplates) { + cout << "Setting names of template parameters" << endl; + for(int iTempl=0; iTemplfNParsTempls; iTempl++) { + cout << "iPar" << iTempl+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg+fNParsVnSgn << ", name:" << Form("w_%s", fKDETemplates[iTempl].GetName()) << endl; + fVnTotFunc->SetParName(iTempl+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl,Form("wm_%s", fKDETemplates[iTempl].GetName())); + fVnTotFunc->SetParName(iTempl+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn,Form("wf_%s", fKDETemplates[iTempl].GetName())); + } + } if(fSecondPeak) { fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn,"SecPeakInt"); fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+1,"SecPeakMean"); fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+2,"SecPeakSigma"); } - fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg,Form("v%dSgn",fHarmonic)); + fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsTempls+fNParsVnBkg,Form("v%dSgn",fHarmonic)); - if(fSecondPeak && fDoSecondPeakVn) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg+1,Form("v%dSecPeak",fHarmonic));} - if(fReflections && fVnRflOpt==kFreePar) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg+1+fNParsVnSecPeak,Form("v%dRefl",fHarmonic));} + if(fSecondPeak && fDoSecondPeakVn) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsTempls+fNParsVnBkg+1,Form("v%dSecPeak",fHarmonic));} + if(fReflections && fVnRflOpt==kFreePar) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsRfl+fNParsSec+fNParsVnBkg+fNParsTempls+1+fNParsVnSecPeak,Form("v%dRefl",fHarmonic));} } //_________________________________________________________________________ @@ -943,6 +1073,16 @@ Double_t VnVsMassFitter::MassRfl(Double_t *m,Double_t *pars){ return pars[0]*value/norm*fRawYieldHelp*fMassHisto->GetBinWidth(1); } +Double_t VnVsMassFitter::MassTemplates(Double_t *m,Double_t *pars){ + // Add the contributions of the templates loaded in fKDETemplates, each + // scaled by a multiplicative constant, left as free fit parameter + Double_t totalTemplates = 0.; + for(int iTempl=0; iTempl #include #include +#include #include "Fit/Fitter.h" #include "Fit/Chi2FCN.h" #include "Math/WrappedMultiTF1.h" @@ -61,6 +62,21 @@ class VnVsMassFitter : public TObject { fMaxRefl=maxRange; fReflections=kTRUE; } + void SetKDETemplates(std::vector templs, std::vector templsnames, + std::vector initweights, std::vector minweights, + std::vector maxweights) { + fKDETemplates=templs; + fInitWeights=initweights; + fWeightsLowerLims=minweights; + fWeightsUpperLims=maxweights; + for(int iFunc=0; iFunc GetMassTemplFuncts() const { + if(fTemplates) return fKDEMassTemplatesDraw; + else return {}; + } + std::vector GetVnTemplFuncts() const { + if(fTemplates) return fKDEVnTemplatesDraw; + else return {}; } //struct for global chi2 (for simultaneus fit) @@ -164,6 +187,7 @@ class VnVsMassFitter : public TObject { Double_t MassSignal(Double_t *m, Double_t *pars); Double_t MassBkg(Double_t *m, Double_t *pars); Double_t MassRfl(Double_t *m,Double_t *par); + Double_t MassTemplates(Double_t *m,Double_t *pars); Double_t MassBkgRfl(Double_t *m,Double_t *par); Double_t MassSecondPeak(Double_t *m,Double_t *par); Double_t vnBkgFunc(Double_t *m, Double_t *pars); @@ -260,6 +284,14 @@ class VnVsMassFitter : public TObject { Bool_t fDoSecondPeakVn; /// flag to introduce second peak vn in the vn vs. mass fit Double_t fVnSecPeakUncertainty; /// vn uncertainty of second peak from fit Int_t fHarmonic; /// harmonic number for drawing + Bool_t fTemplates; /// flag use/not use templates + Int_t fNParsTempls; /// fit parameters to include templates + std::vector fKDETemplates; /// vector to store TKDE to be added as templates to the fit function + std::vector fKDEVnTemplatesDraw; /// vector to store TKDE to be added as templates to the fit function + std::vector fKDEMassTemplatesDraw; /// vector to store TKDE to be added as templates to the fit function + std::vector fWeightsUpperLims; /// upper limit of the templates' weights + std::vector fWeightsLowerLims; /// lower limit of the templates' weights + std::vector fInitWeights; /// init values of the templates' weights /// \cond CLASSDEF ClassDef(VnVsMassFitter,5); diff --git a/run3/flow/kde_producer.py b/run3/flow/kde_producer.py index 14663387..6db62347 100644 --- a/run3/flow/kde_producer.py +++ b/run3/flow/kde_producer.py @@ -7,49 +7,69 @@ import ctypes from ROOT import TFile, TKDE, TCanvas, TH1D -def kde_producer(tree, var, cent, pt_min, pt_max, inv_mass_bins, flag, outfile): +def kde_producer(tree_file, var, pt_min, pt_max, inv_mass_bins, flag, outfile=''): - # convert the tree to a pandas dataframe + print(f"Producing KDE from {tree_file} for var {var}, {pt_min} <= pt < {pt_max}, flag {flag}") + # convert the tree_file to a pandas dataframe dfsData = [] - # print(tree) - with uproot.open(f'{tree}') as f: + with uproot.open(f'{tree_file}') as f: # print(f.keys()) for iKey, key in enumerate(f.keys()): if 'O2hfcanddplite' in key: # print(key) dfData = f[key].arrays(library='pd') dfsData.append(dfData) - + full_dataset = pd.concat([df for df in dfsData], ignore_index=True) - print(full_dataset['fFlagMcMatchRec'].tolist()) - print(full_dataset['fOriginMcRec'].tolist()) - print(full_dataset['fFlagMcDecayChanRec'].tolist()) pt_filtered_df = full_dataset.query(f"{pt_min} <= fPt < {pt_max}") filtered_df = pt_filtered_df.query(f"fFlagMcMatchRec == {flag} or fFlagMcMatchRec == {-flag}") - # filtered_df = pt_filtered_df.query(f"fFlagMcMatchRec == 1") var_values = filtered_df[f'{var}'].tolist() # Or use `.tolist()` to get a list - # print(var_values) + print(var_values) kde = TKDE(len(var_values), np.asarray(var_values, 'd'), 1.7, 2.0) kde_func = kde.GetFunction(1000) - binned_var_values = TH1D(f'hBinned_{var}_{flag}', f'hBinned_{var}_{flag}', 5000, 1, 3) + binned_var_values = TH1D(f'hBinned', f'hBinned', 5000, 1, 3) for var_value in var_values: binned_var_values.Fill(var_value) - cOverlap = TCanvas('cOverlap', 'cOverlap', 600, 600) - cOverlap.cd() - binned_var_values.Draw() - kde_func.Draw('SAME') + if outfile != '': + cOverlap = TCanvas('cOverlap', 'cOverlap', 600, 600) + cOverlap.cd() + binned_var_values.Draw() + kde_func.Draw('SAME') + outfile.mkdir(f'KDE_pT_{pt_low}_{pt_max}_flag{flag}') + outfile.cd(f'KDE_pT_{pt_low}_{pt_max}_flag{flag}') + binned_var_values.Write() + kde_func.Write() + cOverlap.Write() - outfile.mkdir(f'pT_{pt_low}_{pt_max}') - outfile.cd(f'pT_{pt_low}_{pt_max}') - kde_func.Write() - binned_var_values.Write() - cOverlap.Write() + return kde, binned_var_values + +def kde_producer_sim(tree_file, tree_name, pt_min, pt_max, outfile=''): - return kde_func, binned_var_values + print(f"Producing KDE from sim {tree_file} for process {tree_name}, {pt_min} <= pt < {pt_max}") + # convert the tree_file to a pandas dataframe + dfsData = [] + with uproot.open(f'{tree_file}') as f: + # print(f.keys()) + for iKey, key in enumerate(f.keys()): + if f'{tree_name}' in key: + # print(key) + dfData = f[key].arrays(library='pd') + dfsData.append(dfData) + full_dataset = pd.concat([df for df in dfsData], ignore_index=True) + pt_filtered_df = full_dataset.query(f"{pt_min} <= pt < {pt_max}") + var_values = pt_filtered_df["mass"].tolist() # Or use `.tolist()` to get a list + + kde = TKDE(len(var_values), np.asarray(var_values, 'd'), 1.7, 2.0) + + binned_var_values = TH1D(f'hBinned', f'hBinned', 5000, 1, 3) + for var_value in var_values: + binned_var_values.Fill(var_value) + + return kde if __name__ == "__main__": parser = argparse.ArgumentParser(description="Arguments") @@ -57,8 +77,6 @@ def kde_producer(tree, var, cent, pt_min, pt_max, inv_mass_bins, flag, outfile): default="fM", help="variable of interest") parser.add_argument("--config", "-cfg", metavar="text", default="config.yaml", help="configuration file") - parser.add_argument("--centrality", "-c", metavar="text", - default="k2050", help="centrality class") parser.add_argument("--ptmin", "-pmin", metavar="text", default="2.", help="min pt") parser.add_argument("--ptmax", "-pmax", metavar="text", @@ -68,7 +86,7 @@ def kde_producer(tree, var, cent, pt_min, pt_max, inv_mass_bins, flag, outfile): parser.add_argument("--flag", "-f", metavar="chn flag", default="2", help="channel flag") parser.add_argument("--input", "-in", metavar="path/input.root", - default="AnalysisResults.root", help="path to tree") + default="AnalysisResults.root", help="path to file containing tree") parser.add_argument("--wagon_id", "-w", metavar="text", default="", help="wagon ID", required=False) parser.add_argument("--outputdir", "-o", metavar="text", @@ -91,11 +109,11 @@ def kde_producer(tree, var, cent, pt_min, pt_max, inv_mass_bins, flag, outfile): if args.config != parser.get_default("config"): for pt_low, pt_max, inv_mass_bins in zip(config['pt_mins'], config['pt_maxs'], config['inv_mass_bins']): - KDE, histo = kde_producer(config['input'], config['variable'], config['centrality'], \ + KDE, histo = kde_producer(config['input'], config['variable'], \ pt_low, pt_max, inv_mass_bins, config['chn_flag'], outfile) else: - KDE, histo = kde_producer(args.input, args.var, args.centrality, \ + KDE, histo = kde_producer(args.input, args.var, \ args.ptmin, args.ptmax, args.invmassbins, args.flag, outfile) outfile.Close() From 22dcf1c66138fde3a17002455f2913429cc9c308 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Sun, 17 Nov 2024 19:24:53 +0100 Subject: [PATCH 4/6] Eliminated reduntant fit parameters --- run3/flow/flow_analysis_utils.py | 59 +++++++++++++++++++- run3/flow/get_vn_vs_mass.py | 42 ++++++++------- run3/flow/invmassfitter/InvMassFitter.cxx | 20 +++---- run3/flow/invmassfitter/VnVsMassFitter.cxx | 63 ++++++++++++++-------- run3/flow/invmassfitter/VnVsMassFitter.h | 3 +- run3/flow/kde_producer.py | 18 +++++-- run3/flow/project_thnsparse.py | 15 +++++- run3/flow/run_full_flow_analysis.py | 1 + 8 files changed, 157 insertions(+), 64 deletions(-) diff --git a/run3/flow/flow_analysis_utils.py b/run3/flow/flow_analysis_utils.py index c780a98c..4c63b836 100644 --- a/run3/flow/flow_analysis_utils.py +++ b/run3/flow/flow_analysis_utils.py @@ -1,11 +1,12 @@ ''' Analysis utilities for flow analysis ''' -import ROOT + import os import sys import ctypes from itertools import combinations +import ROOT import numpy as np def get_vn_versus_mass(thnSparse, inv_mass_bins, mass_axis, vn_axis, debug=False): @@ -28,6 +29,9 @@ def get_vn_versus_mass(thnSparse, inv_mass_bins, mass_axis, vn_axis, debug=False - hist_mass_proj: TH1D, histogram with vn as a function of mass ''' + print('GETTING VNVSMASS') + print(f'AXIS MASS: {mass_axis}') + print(f'AXIS VN: {vn_axis}') hist_vn_proj = thnSparse.Projection(vn_axis, mass_axis) hist_mass_proj = thnSparse.Projection(mass_axis) hist_mass_proj.Reset() @@ -53,6 +57,59 @@ def get_vn_versus_mass(thnSparse, inv_mass_bins, mass_axis, vn_axis, debug=False return hist_mass_proj +def get_occupancy(thnSparse, occupancy_axis, debug=False): + ''' + Project occupancy versus mass + + Input: + - thnSparse: + THnSparse, input THnSparse obeject (already projected in centrality and pt) + - occupancy_axis: + int, axis number for occupancy + - debug: + bool, if True, create a debug file with the projections (default: False) + + Output: + - hist_occupancy: + TH1D, histogram with vn as a function of mass + ''' + print('GETTING OCCUPANCY') + print(f'AXIS: {occupancy_axis}') + hist_occupancy = thnSparse.Projection(occupancy_axis) + + if debug: + outfile = ROOT.TFile('debug.root', 'RECREATE') + hist_occupancy.Write() + outfile.Close() + + return hist_occupancy + +def get_evselbits(thnSparse, evselbits_axis, debug=False): + ''' + Project evselbits versus mass + + Input: + - thnSparse: + THnSparse, input THnSparse obeject (already projected in centrality and pt) + - evselbits_axis: + int, axis number for evselbits + - debug: + bool, if True, create a debug file with the projections (default: False) + + Output: + - hist_evselbits: + TH1D, histogram with vn as a function of mass + ''' + print('GETTING EVSELBITS') + print(f'AXIS: {evselbits_axis}') + hist_evselbits = thnSparse.Projection(evselbits_axis) + + if debug: + outfile = ROOT.TFile('debug.root', 'RECREATE') + hist_evselbits.Write() + outfile.Close() + + return hist_evselbits def get_resolution(dets, det_lables, cent_min_max): ''' diff --git a/run3/flow/get_vn_vs_mass.py b/run3/flow/get_vn_vs_mass.py index 044738e3..e0cb323b 100644 --- a/run3/flow/get_vn_vs_mass.py +++ b/run3/flow/get_vn_vs_mass.py @@ -135,7 +135,7 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, elif fitConfig.get('ProduceNew') and fitConfig.get('FromSim'): KDEtemplates[iPt][iFlag] = kde_producer_sim(fitConfig['MCTemplatesInputs'][iPt][iFlag], fitConfig['MCTemplatesTreeNames'][iPt][iFlag], - ptMins[iPt], ptMaxs[iPt]) + ptMins[iPt], ptMaxs[iPt], fitConfig['MCTemplatesFlags'][iPt][iFlag]) KDEtemplatesNames[iPt][iFlag] = fitConfig['MCTemplatesTreeNames'][iPt][iFlag] else: templFile = TFile.Open(f'{fitConfig["MCTemplatesFiles"][iPt][iFlag]}', 'r') @@ -407,7 +407,7 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, fBkgFuncMass.append(vnResults['fBkgFuncMass']) fBkgFuncVn.append(vnResults['fBkgFuncVn']) fMassTemplFuncts[iPt] = vnResults['fMassTemplFuncts'] - fVnTemplFuncts[iPt] = vnResults['fVnTemplFuncts'] + # fVnTemplFuncts[iPt] = vnResults['fVnTemplFuncts'] if useRefl: fMassBkgRflFunc.append(vnResults['fMassBkgRflFunc']) @@ -478,10 +478,12 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, latex.DrawLatex(0.18, 0.20, f'RoverS = {SoverR:.2f}') print(f'Drawing fMassTemplFuncts') print(fMassTemplFuncts) - for iMassTemplFunct in fMassTemplFuncts[iPt]: - print(f'Eval: {iMassTemplFunct.Eval(1.8)}') - iMassTemplFunct.SetLineColor(1) - iMassTemplFunct.Draw('same') + for iMassTemplFunct, massTemplFunct in enumerate(fMassTemplFuncts[iPt]): + print("Setting object styleeee") + SetObjectStyle(massTemplFunct, color=kMagenta+iMassTemplFunct*2, linewidth=3) + print(f'Eval: {massTemplFunct.Eval(1.8)}') + massTemplFunct.SetLineColor(1) + massTemplFunct.Draw('same') cSimFit[iCanv].Modified() cSimFit[iCanv].Update() cSimFit[iPt].cd(2) @@ -500,13 +502,13 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, if secPeak: latex.DrawLatex(0.18, 0.75, f'#it{{v}}{harmonic}(D^{{+}}) = {vnResults["vnSecPeak"]:.3f} #pm {vnResults["vnSecPeakUnc"]:.3f}') - print(f'Drawing fVnTemplFuncts') - print(fVnTemplFuncts) - for iVnTemplFunct in fVnTemplFuncts[iPt]: - print(f'Eval: {iVnTemplFunct.Eval(1.8)}') - iVnTemplFunct.Draw('same') - cSimFit[iCanv].Modified() - cSimFit[iCanv].Update() + # print(f'Drawing fVnTemplFuncts') + # print(fVnTemplFuncts) + # for iVnTemplFunct in fVnTemplFuncts[iPt]: + # print(f'Eval: {iVnTemplFunct.Eval(1.8)}') + # iVnTemplFunct.Draw('same') + # cSimFit[iCanv].Modified() + # cSimFit[iCanv].Update() cSimFit[iCanv].Modified() cSimFit[iCanv].Update() @@ -737,8 +739,8 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, gvnUncSecPeak.Draw('pez same') canvVnUnc.Modified() canvVnUnc.Update() - if not batch: - input('Press Enter to continue...') + # if not batch: + # input('Press Enter to continue...') #save output histos print(f'Saving output to {outputdir}') @@ -750,7 +752,10 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, suffix_pdf = ')' else: suffix_pdf = '' - cSimFit[iPt].SaveAs(f'{outputdir}/SimFit{suffix}_{particleName}.pdf{suffix_pdf}') + if len(ptMins)==1: + cSimFit[iPt].SaveAs(f'{outputdir}/SimFit{suffix}_{particleName}.pdf') + else: + cSimFit[iPt].SaveAs(f'{outputdir}/SimFit{suffix}_{particleName}.pdf{suffix_pdf}') outfile_name = f'{outputdir}/raw_yields{suffix}.root' outFile = TFile(outfile_name, 'recreate') if vn_method == 'sp' or vn_method == 'ep': @@ -825,9 +830,8 @@ def get_vn_vs_mass(fitConfigFileName, centClass, inFileName, outFile.Close() - if not batch: - input('Press enter to exit') - + # if not batch: + # input('Press enter to exit') if __name__ == "__main__": parser = argparse.ArgumentParser(description='Arguments') diff --git a/run3/flow/invmassfitter/InvMassFitter.cxx b/run3/flow/invmassfitter/InvMassFitter.cxx index cf63bfed..6158eb9f 100644 --- a/run3/flow/invmassfitter/InvMassFitter.cxx +++ b/run3/flow/invmassfitter/InvMassFitter.cxx @@ -259,16 +259,15 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){ /// returns 0 if the fit fails /// returns 1 if the fit succeeds /// returns 2 if there is no signal and the fit is performed with only background - cout << "Entering InvMassFitter::MassFitter" << endl; TVirtualFitter::SetDefaultFitter("Minuit"); - // for(int iFunc=0; iFuncIntegral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width"); fOnlySideBands = kTRUE; @@ -338,11 +337,6 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){ // } fTotFunc = CreateTotalFitFunction("funcmass"); - TFile *funcFile = new TFile("~/flowDplus/FitFuncInvMassFitter.root", "recreate"); - funcFile->cd(); - fTotFunc->Write(); - funcFile->Close(); - if(doFinalFit){ printf("\n--- Final fit with signal+background on the full range ---\n"); TFitResultPtr resultptr=fHistoInvMass->Fit("funcmass",Form("R,S,%s,+,0",fFitOption.Data())); diff --git a/run3/flow/invmassfitter/VnVsMassFitter.cxx b/run3/flow/invmassfitter/VnVsMassFitter.cxx index ab269fc9..1c630374 100755 --- a/run3/flow/invmassfitter/VnVsMassFitter.cxx +++ b/run3/flow/invmassfitter/VnVsMassFitter.cxx @@ -234,7 +234,7 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { Int_t NvnParsSgn = 1; if(fSecondPeak && fDoSecondPeakVn) {NvnParsSgn+=1;} if(fReflections && fVnRflOpt==kFreePar) {NvnParsSgn+=1;} - if(fTemplates) {NvnParsSgn+=fNParsTempls;} + // if(fTemplates) {NvnParsSgn+=fNParsTempls;} cout << "Added " << fNParsTempls << " to vn sgn pars" << endl; const Int_t nparsvn = nparsmass+fNParsVnBkg+NvnParsSgn; cout << "nparsvn " << nparsvn << endl; @@ -332,17 +332,17 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { } if(fTemplates) { for(int iTemplPar=0; iTemplParfInitWeights.size()/2; iTemplPar++) { - cout << "Setting template parameter " << iTemplPar << endl; - cout << "Mass func parameter " << iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl << endl; - cout << fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Name() << endl; - cout << "[" << fInitWeights[iTemplPar] << ", " << fWeightsLowerLims[iTemplPar] - << ", " << fWeightsUpperLims[iTemplPar] << "]" << endl; - cout << "Vn func parameter " << iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn << endl; - cout << "[" << fInitWeights[iTemplPar+this->fInitWeights.size()/2] << ", " << fWeightsLowerLims[iTemplPar+this->fInitWeights.size()/2] - << ", " << fWeightsUpperLims[iTemplPar+this->fInitWeights.size()/2] << "]" << endl; - cout << fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).Name() << endl; - - cout << "---" << endl; + // cout << "Setting template parameter " << iTemplPar << endl; + // cout << "Mass func parameter " << iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl << endl; + // cout << fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Name() << endl; + // cout << "[" << fInitWeights[iTemplPar] << ", " << fWeightsLowerLims[iTemplPar] + // << ", " << fWeightsUpperLims[iTemplPar] << "]" << endl; + // cout << "Vn func parameter " << iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn << endl; + // cout << "[" << fInitWeights[iTemplPar+this->fInitWeights.size()/2] << ", " << fWeightsLowerLims[iTemplPar+this->fInitWeights.size()/2] + // << ", " << fWeightsUpperLims[iTemplPar+this->fInitWeights.size()/2] << "]" << endl; + // cout << fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).Name() << endl; + + // cout << "---" << endl; if(this->fWeightsLowerLims[iTemplPar] > this->fWeightsUpperLims[iTemplPar]) { cout << "Fix mass" << endl; fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); @@ -420,9 +420,12 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { } if(fTemplates) { int nParBeforeTemplsMass = fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl; + cout << "Size of fKDETemplates: " << this->fKDETemplates.size() << endl; + cout << "Eval func 1: " << this->fKDETemplates[0].Eval(1.8) << endl; + // cout << "Eval func 2: " << this->fKDETemplates[1].Eval(1.8) << endl; cout << "Coeff first func mass: " << result.Parameter(nParBeforeTemplsMass) << endl; - cout << "Eval func: " << this->fKDETemplates[0].Eval(1.8) << endl; - for(int iTempl=0; iTemplfKDETemplates[fKDETemplates.size() / 2].Eval(1.8) << endl; - for(int iTempl=0; iTemplfKDETemplates[iTempl + fKDETemplates.size() / 2].Eval(x[0]); + double paramValue = result.Parameter(iTempl + fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg); + double kdeEval = this->fKDETemplates[iTempl].Eval(x[0]); return paramValue * kdeEval; }, fMassMin, fMassMax, 0)); } @@ -618,7 +621,9 @@ Bool_t VnVsMassFitter::MassPrefit() { if(fFixRflOverSig) {fMassFitter->SetFixReflOverS(fRflOverSig);} } if(fTemplates) {fMassFitter->SetTemplates(fKDETemplates, fInitWeights, fWeightsLowerLims, fWeightsUpperLims);} + cout << "Fitting" << endl; Bool_t status = fMassFitter->MassFitter(kFALSE); + cout << "Status " << status << endl; if(status) { fMassFuncFromPrefit = (TF1*)fMassFitter->GetMassFunc(); @@ -629,7 +634,6 @@ Bool_t VnVsMassFitter::MassPrefit() { } if(fReflections) fRawYieldHelp=fMassFitter->GetRawYield(); - cout << "Status " << status << endl; return status; } @@ -1073,6 +1077,7 @@ Double_t VnVsMassFitter::MassRfl(Double_t *m,Double_t *pars){ return pars[0]*value/norm*fRawYieldHelp*fMassHisto->GetBinWidth(1); } +//_________________________________________________________________________ Double_t VnVsMassFitter::MassTemplates(Double_t *m,Double_t *pars){ // Add the contributions of the templates loaded in fKDETemplates, each // scaled by a multiplicative constant, left as free fit parameter @@ -1083,6 +1088,17 @@ Double_t VnVsMassFitter::MassTemplates(Double_t *m,Double_t *pars){ return totalTemplates; } +//_________________________________________________________________________ +Double_t VnVsMassFitter::VnTemplates(Double_t *m,Double_t *pars){ + // Add the contributions of the templates loaded in fKDETemplates, each + // scaled by a multiplicative constant, left as free fit parameter + Double_t totalTemplates = 0.; + for(int iTempl=0; iTempl Date: Sat, 4 Jan 2025 16:52:59 +0100 Subject: [PATCH 5/6] fixed representation of v2 components --- run3/flow/WriteTF1.root | Bin 0 -> 3929 bytes run3/flow/get_vn_vs_mass.py | 26 ++++---- run3/flow/invmassfitter/VnVsMassFitter.cxx | 73 +++++++++------------ run3/flow/invmassfitter/VnVsMassFitter.h | 1 + 4 files changed, 46 insertions(+), 54 deletions(-) create mode 100644 run3/flow/WriteTF1.root diff --git a/run3/flow/WriteTF1.root b/run3/flow/WriteTF1.root new file mode 100644 index 0000000000000000000000000000000000000000..8725abbe2ae867a619155d58b51df8535591ff4f GIT binary patch literal 3929 zcma)9cQ72@{$0IBCy3syx_T!;bV7*U%WARFMU5^*jowQ_MAT)~)uM-psEbvivxv5d z_~LE+{`k$jc{A_bna?eA?>V3Q*O_yCAP_$Q0KE+W0Gt2-a&t_vVfM0^!orjiE@pBD z0B|u6h>!rR5uG%PF(~D68m-u)iN5Y*%}}c$pdK`96>i~!T8frTAa5TRBL94rhF7eFH!g+ndnry9Vs}nY9QB0c2sd<| z?ESoFjWp{7!zld^Ltu!ox@#dSvmwXH7nSopD~D~023t*dTXh+In|8I+ zM=QeeA=mVs>hr1Njje4h*E-`2jv+e>i3Dh@ZFILlz`HRU??!TxwvSLY?p38Sm`VBL z+)1VDWdA>1^-Ad2PxsM^^ZZRBtI7{wm-$A8MbrNBAN}KcH7A6a1f$x)t8Y6Ry}hrw zuC(&wLb8WX5C%*C)T6eo#IFz`IoA41wH`fo-y@@Un(G!s^iii*i;{!kdDWAlTZ_JE zfp9D27Qc~wPq0gmPT&E~l|NHWfRcr^ThBu@+r0Z*nych;ff^G1+$Bip*68(R$eZg* zM4^Ush~*+=E)Dqfgc7!@p!G8%Y1r~=LORDJ3zW0~xv}*=9L%^xTjnYgbQG_EdoJ76 zZ_gQ?qNMAMyT)Z6PK86OeyxKe3Q^L`((d+Gp#$M>jjZ<05nbp96?A9C9n_9JWo$i1 z$BcAW{#E70$lEeVPvtV_C#%w-@U8fE@6;C6!yI{?eVgrk*PeGavFt$itm}T}aTmb| zX~d)>Vq3`F-AQ>vi&VJ;9D7akxhXeq-sh0Ffo_sk?vux=x3rjh`7!BLCBdxum&#CB z{rYs9znb(=H`kJesmM51662f-3HTCA^sH?j%JM7!VeGB_JjYGv1 z-W6gw-+`6BP3crvO&)x7TPisr`XxF);twN>aCWolnyA#GjBh#oMY`XgT6iTIg-?uM z5^}qAJttE8wlH|%D9r5S`9X!qAq$H|)0IFE6uC7m5WWONEHWc!%1<~Ubx5>R(kgoa z2kKd2FncC{&}Xl0y+aL?t;7K9dg zm*r=XL<9?jtmkd%0l4#Wb9o`*=p~GV$`q~u;_QQYAkhowXoscmF-q)%$zCub=@x0J zasHMIqo222-TScGlcn)+G*-fl` zw6;{juPebYFlQgJCU-@CUa{fi75oN1e+)lral6G~EG1?bZqUnGRTwa|*5%lfaF)2= z#=d?CZl2&-Og=&ow0Tv3tkRR<LumT&d~l)wZk(cJYs-3Oo$hS+>U*xb#w*=1F#M ze3vykYn`yn(oZX>%*Vt0YX>V&M;{U744=$mCnUhyO|g&lfjeK1W>vK=+jWrv=cwWK zbx;O=V=K)_+|Q^19(UNN#gKJI?1A2pqSFtuHgFf$%F{C{th3+uD0%fWh8(!xFse1< zfRA4#j<*Cn-e%@GyTm3DRxm5l&LFX2+ncTx_f7LGyJh+#AXU|Y7$6wuv3eD6*U@Xn zP8Swu#qBG1{@GQtMa^3Bgxr}NMs#$*f^*_S&z~PWV5U14nrA=ok*I8vX(NALGh6%( z?JSB}wK1Z!sPT62-a^SmXbH=W+Yf5)IiBEeF6g~_BuX5z%eKUD1xb&s@<6E&d_f81 zjo*;XQjs~VQRkLwut`N?HEYT9rL1^>nULTVG+lJ$=X?C*@Ljc5$`3=f|EwJk2%6BH z#GX!_AqBUT7QYGRtvArT)>zgBo=If1l)QRZ?44W51b)keP?Jy!*GiZiMvWAC$ZGBHt{r@Jx^y|u0W32U9+ z_L?z!dZMwCp>=k79PR@W={GT)*ja1-#kBD=)AoHY-I!SQHUybFDELl&UwfKekiBS{Ii10+f0p{ zHgvIWvt+CrdNI%7wHzbP^W#H6y+*|zqC+^p9^U>iTkl+|f*CZ0zkHx4vJl8_m7e|~ z*6rJfXI^hpI<+8wE6WV8*8!^~qrJu1d%YRzVra*)cFl$E&;8&~=@IPUX9U>}5{Zr3 z)Q80e1(JZcx?bvKvGjBw8rlJtyp5|Si~cE5NQN`DkXSM|1XuU7$5FmsbnuVZ$lho* zSt@+MF70MK<0gx_)RZ;8EO_xFtE2gK2Pl(bagZY!v@`S1Qz^^D!nr2Xz);H3XIzdn zB`Mz%g6B8#=Nuz1avc@B?e;-yK6>*>Bla=H^R^{C2No z>V8tz7Fl~(t|79v7*vh_Tx^#ry4xixY90zDkWaYRPNT&j0P&u3jnT`@Gyrtj;04d% zviMZxz<=u}cC26a`-{pLPTD%=7UYJNvFuT$&a9`ih+`2TlWf)`r?-{D?gug}XQ32@ zly#|EU@y1RzPjV*AxjufbKXl1Ib1}IRUN<8HmOnAoZe20U)Ir=pjNl1Vcj&T6sy?; z&l>s0sjBgd8sA{iv9v!4ITEJJvG)^83a$Fa{5l)I?DO#R;GEQp^JvVEsE<^<9+Li8Afscg1kTZqcg8McbBsWh>Mpm=_mUq zDg@<4e*B)!J?JhL=QkI%7|||W87xtXLtcuy%PDDh4fS=6f=ptB$k6dt!>NXZwXfs+ zie?OqTc9LhcFr8n(VK?S5E-q!+RPcYtg6t4Z~N$OT|zoCzmJCzo| zT~pVAl16U2a{CPYezc`mo+rn+xUUJYfBj0(>KK{<$<-$7fU4+*?X(>UWjryykf&8T zSZJga_#*b?5@@7beF2ns#!4v?C_7Wuegn$Yfg0A=R%{tJxbXcE_*l9sVm6gu%{Qr3 z#b7IOF6aQMlHp4L#-en zdE-Xdyv-%OoTl^f@?Hai*wLUJ@g=VR;;VTr|RjAKN_jnG;Q-?yQ5V zRQKNPhY6B|8`i+ulJx<+3klg971@h`^WmXL3-xo`V;gl*BJ}Mb5c_R`qqlb06F!~m z_oQ9ZMQD(s97`?YWnGy%5E^?_| zr~kQDz4QH6@p;abXO18dvvi|IP{rhH6^llDKSZ_8eYmPp z`hAfq5Ni0(K|?`gjpy(E7H6V;&PHroE3&7aWgM_1-LC!8awD!BiZ#yh$`RiNP2=YC z?RrA(EN%HEo=aYW_c(amgpVqOhLdzz<<(2N=JMl%XfGq;$tjC=&yIGT%A$$W7Z!dT zUYt3tyTwf)TgFbx&1%J^T0A1OLZ)P!JTFvLmFtRZ9}-&fHV%C9DxEC|I%1&J zJPJBfb{wI&Co7BAuf8>GIWi!NjXbzK z{eZC5VpNWFa6#TCXi>iFPF3O^xCT8uhmhvtv}0_X^fInitWeights.size()/2; iTemplPar++) { - // cout << "Setting template parameter " << iTemplPar << endl; - // cout << "Mass func parameter " << iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl << endl; - // cout << fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Name() << endl; - // cout << "[" << fInitWeights[iTemplPar] << ", " << fWeightsLowerLims[iTemplPar] - // << ", " << fWeightsUpperLims[iTemplPar] << "]" << endl; - // cout << "Vn func parameter " << iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn << endl; - // cout << "[" << fInitWeights[iTemplPar+this->fInitWeights.size()/2] << ", " << fWeightsLowerLims[iTemplPar+this->fInitWeights.size()/2] - // << ", " << fWeightsUpperLims[iTemplPar+this->fInitWeights.size()/2] << "]" << endl; - // cout << fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).Name() << endl; - - // cout << "---" << endl; + for(int iTemplPar=0; iTemplParfInitWeights.size(); iTemplPar++) { + cout << "Setting template parameter " << iTemplPar << endl; + cout << "Mass func parameter " << iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl << endl; + cout << fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Name() << endl; + cout << "[" << fInitWeights[iTemplPar] << ", " << fWeightsLowerLims[iTemplPar] + << ", " << fWeightsUpperLims[iTemplPar] << "]" << endl; + if(this->fWeightsLowerLims[iTemplPar] > this->fWeightsUpperLims[iTemplPar]) { cout << "Fix mass" << endl; fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); @@ -352,15 +347,6 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetLimits(fWeightsLowerLims[iTemplPar],fWeightsUpperLims[iTemplPar]); } - if(this->fWeightsLowerLims[iTemplPar] > this->fWeightsUpperLims[iTemplPar]) { - cout << "Fix vn" << endl; - fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).SetValue(fInitWeights[iTemplPar+this->fInitWeights.size()/2]); - fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).Fix(); - } else { - cout << "Vary vn" << endl; - fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).SetValue(fInitWeights[iTemplPar+this->fInitWeights.size()/2]); - fitter.Config().ParSettings(iTemplPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn).SetLimits(fWeightsLowerLims[iTemplPar+this->fInitWeights.size()/2],fWeightsUpperLims[iTemplPar+this->fInitWeights.size()/2]); - } cout << "---" << endl; } } @@ -393,6 +379,7 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { fMassTotFunc->SetParError(iPar,result.ParError(iPar)); } if(iPar>=nparsmass && iParSetParameter(iPar-nparsmass,result.Parameter(iPar)); fVnBkgFunc->SetParError(iPar-nparsmass,result.ParError(iPar)); } @@ -419,34 +406,34 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { } } if(fTemplates) { - int nParBeforeTemplsMass = fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl; - cout << "Size of fKDETemplates: " << this->fKDETemplates.size() << endl; - cout << "Eval func 1: " << this->fKDETemplates[0].Eval(1.8) << endl; - // cout << "Eval func 2: " << this->fKDETemplates[1].Eval(1.8) << endl; - cout << "Coeff first func mass: " << result.Parameter(nParBeforeTemplsMass) << endl; - cout << "Coeff second func mass: " << result.Parameter(nParBeforeTemplsMass+1) << endl; + int idxParMassTemplsScaling = fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl; + int idxParVnSgn = idxParMassTemplsScaling+fNParsTempls+fNParsVnBkg; + double vnSgn = result.Parameter(idxParVnSgn); for(int iTempl=0; iTemplfKDETemplates[iTempl].Eval(x[0]); - return paramValue * kdeEval; + [&, this, iTempl, idxParMassTemplsScaling, result] (double *x, double *par) { + double templScalingPar = result.Parameter(iTempl + idxParMassTemplsScaling); + double kdeTemplEval = this->fKDETemplates[iTempl].Eval(x[0]); + return templScalingPar * kdeTemplEval; }, fMassMin, fMassMax, 0)); + + fKDEVnTemplatesDraw.push_back(new TF1(Form("vnTempl_%s", fKDETemplates[iTempl].GetName()), + [&, this, iTempl, result, vnSgn] (double *x, double *par) { + return (vnSgn * this->fKDEMassTemplatesDraw[iTempl]->Eval(x[0])) / (this->fMassTotFunc->Eval(x[0])); + }, fMassMin, fMassMax, 0)); } - int nParBeforeTemplsVn = nParBeforeTemplsMass+fNParsTempls+fNParsVnBkg+fNParsVnSgn; - cout << "Coeff first func vn: " << result.Parameter(nParBeforeTemplsVn) << endl; - cout << "Coeff second func vn: " << result.Parameter(nParBeforeTemplsMass+1) << endl; - for(int iTempl=0; iTemplfKDETemplates[iTempl].Eval(x[0]); - return paramValue * kdeEval; + + fKDEVnTemplatesDraw.push_back(new TF1("vnBkg", + [&, this, result, vnSgn] (double *x, double *par) { + return (this->fVnBkgFunc->Eval(x[0]) * this->fMassBkgFunc->Eval(x[0])) / (this->fMassTotFunc->Eval(x[0])); + }, fMassMin, fMassMax, 0)); + + fKDEVnTemplatesDraw.push_back(new TF1("vnSgn", + [&, this, result, vnSgn] (double *x, double *par) { + return (vnSgn * this->fMassSgnFunc->Eval(x[0])) / (this->fMassTotFunc->Eval(x[0])); }, fMassMin, fMassMax, 0)); - } } - cout << "DrawFit: " << drawFit << endl; if(drawFit) {DrawFit();} fVn = fVnTotFunc->GetParameter(fVnTotFunc->GetNpar()-NvnParsSgn); @@ -1245,7 +1232,7 @@ Double_t VnVsMassFitter::vnFunc(Double_t *m, Double_t *pars) { Double_t Templates = 0; // if(fTemplates) {Templates += VnTemplates(m,vntemplpars);} if(fTemplates) {Templates += VnTemplates(m,templpars);} - + return (vnSgn*Sgn+vnBkg*Bkg+vnSecPeak*SecPeak+vnRefl*Refl+vnSgn*Templates)/(Sgn+Bkg+SecPeak+Refl+Templates); } diff --git a/run3/flow/invmassfitter/VnVsMassFitter.h b/run3/flow/invmassfitter/VnVsMassFitter.h index d40b5eea..68137ba7 100755 --- a/run3/flow/invmassfitter/VnVsMassFitter.h +++ b/run3/flow/invmassfitter/VnVsMassFitter.h @@ -153,6 +153,7 @@ class VnVsMassFitter : public TObject { TF1* GetMassBkgRflFunc() const { if(fReflections) return fMassBkgRflFunc; else return nullptr; + } std::vector GetMassTemplFuncts() const { if(fTemplates) return fKDEMassTemplatesDraw; else return {}; From 1d9f35b87707299bb65b872462bab98f066fb246 Mon Sep 17 00:00:00 2001 From: marcellocosti Date: Sat, 4 Jan 2025 17:06:16 +0100 Subject: [PATCH 6/6] Removed leftovers --- run3/flow/WriteTF1.root | Bin 3929 -> 0 bytes run3/flow/flow_analysis_utils.py | 9 +-- run3/flow/get_vn_vs_mass.py | 82 +++++++++------------ run3/flow/invmassfitter/InvMassFitter.cxx | 18 ----- run3/flow/invmassfitter/VnVsMassFitter.cxx | 79 ++------------------ run3/flow/invmassfitter/VnVsMassFitter.h | 5 +- run3/flow/kde_producer.py | 44 ++++------- run3/flow/project_thnsparse.py | 20 ++--- run3/flow/run_full_flow_analysis.py | 1 - 9 files changed, 66 insertions(+), 192 deletions(-) delete mode 100644 run3/flow/WriteTF1.root diff --git a/run3/flow/WriteTF1.root b/run3/flow/WriteTF1.root deleted file mode 100644 index 8725abbe2ae867a619155d58b51df8535591ff4f..0000000000000000000000000000000000000000 GIT binary patch literal 0 HcmV?d00001 literal 3929 zcma)9cQ72@{$0IBCy3syx_T!;bV7*U%WARFMU5^*jowQ_MAT)~)uM-psEbvivxv5d z_~LE+{`k$jc{A_bna?eA?>V3Q*O_yCAP_$Q0KE+W0Gt2-a&t_vVfM0^!orjiE@pBD z0B|u6h>!rR5uG%PF(~D68m-u)iN5Y*%}}c$pdK`96>i~!T8frTAa5TRBL94rhF7eFH!g+ndnry9Vs}nY9QB0c2sd<| z?ESoFjWp{7!zld^Ltu!ox@#dSvmwXH7nSopD~D~023t*dTXh+In|8I+ zM=QeeA=mVs>hr1Njje4h*E-`2jv+e>i3Dh@ZFILlz`HRU??!TxwvSLY?p38Sm`VBL z+)1VDWdA>1^-Ad2PxsM^^ZZRBtI7{wm-$A8MbrNBAN}KcH7A6a1f$x)t8Y6Ry}hrw zuC(&wLb8WX5C%*C)T6eo#IFz`IoA41wH`fo-y@@Un(G!s^iii*i;{!kdDWAlTZ_JE zfp9D27Qc~wPq0gmPT&E~l|NHWfRcr^ThBu@+r0Z*nych;ff^G1+$Bip*68(R$eZg* zM4^Ush~*+=E)Dqfgc7!@p!G8%Y1r~=LORDJ3zW0~xv}*=9L%^xTjnYgbQG_EdoJ76 zZ_gQ?qNMAMyT)Z6PK86OeyxKe3Q^L`((d+Gp#$M>jjZ<05nbp96?A9C9n_9JWo$i1 z$BcAW{#E70$lEeVPvtV_C#%w-@U8fE@6;C6!yI{?eVgrk*PeGavFt$itm}T}aTmb| zX~d)>Vq3`F-AQ>vi&VJ;9D7akxhXeq-sh0Ffo_sk?vux=x3rjh`7!BLCBdxum&#CB z{rYs9znb(=H`kJesmM51662f-3HTCA^sH?j%JM7!VeGB_JjYGv1 z-W6gw-+`6BP3crvO&)x7TPisr`XxF);twN>aCWolnyA#GjBh#oMY`XgT6iTIg-?uM z5^}qAJttE8wlH|%D9r5S`9X!qAq$H|)0IFE6uC7m5WWONEHWc!%1<~Ubx5>R(kgoa z2kKd2FncC{&}Xl0y+aL?t;7K9dg zm*r=XL<9?jtmkd%0l4#Wb9o`*=p~GV$`q~u;_QQYAkhowXoscmF-q)%$zCub=@x0J zasHMIqo222-TScGlcn)+G*-fl` zw6;{juPebYFlQgJCU-@CUa{fi75oN1e+)lral6G~EG1?bZqUnGRTwa|*5%lfaF)2= z#=d?CZl2&-Og=&ow0Tv3tkRR<LumT&d~l)wZk(cJYs-3Oo$hS+>U*xb#w*=1F#M ze3vykYn`yn(oZX>%*Vt0YX>V&M;{U744=$mCnUhyO|g&lfjeK1W>vK=+jWrv=cwWK zbx;O=V=K)_+|Q^19(UNN#gKJI?1A2pqSFtuHgFf$%F{C{th3+uD0%fWh8(!xFse1< zfRA4#j<*Cn-e%@GyTm3DRxm5l&LFX2+ncTx_f7LGyJh+#AXU|Y7$6wuv3eD6*U@Xn zP8Swu#qBG1{@GQtMa^3Bgxr}NMs#$*f^*_S&z~PWV5U14nrA=ok*I8vX(NALGh6%( z?JSB}wK1Z!sPT62-a^SmXbH=W+Yf5)IiBEeF6g~_BuX5z%eKUD1xb&s@<6E&d_f81 zjo*;XQjs~VQRkLwut`N?HEYT9rL1^>nULTVG+lJ$=X?C*@Ljc5$`3=f|EwJk2%6BH z#GX!_AqBUT7QYGRtvArT)>zgBo=If1l)QRZ?44W51b)keP?Jy!*GiZiMvWAC$ZGBHt{r@Jx^y|u0W32U9+ z_L?z!dZMwCp>=k79PR@W={GT)*ja1-#kBD=)AoHY-I!SQHUybFDELl&UwfKekiBS{Ii10+f0p{ zHgvIWvt+CrdNI%7wHzbP^W#H6y+*|zqC+^p9^U>iTkl+|f*CZ0zkHx4vJl8_m7e|~ z*6rJfXI^hpI<+8wE6WV8*8!^~qrJu1d%YRzVra*)cFl$E&;8&~=@IPUX9U>}5{Zr3 z)Q80e1(JZcx?bvKvGjBw8rlJtyp5|Si~cE5NQN`DkXSM|1XuU7$5FmsbnuVZ$lho* zSt@+MF70MK<0gx_)RZ;8EO_xFtE2gK2Pl(bagZY!v@`S1Qz^^D!nr2Xz);H3XIzdn zB`Mz%g6B8#=Nuz1avc@B?e;-yK6>*>Bla=H^R^{C2No z>V8tz7Fl~(t|79v7*vh_Tx^#ry4xixY90zDkWaYRPNT&j0P&u3jnT`@Gyrtj;04d% zviMZxz<=u}cC26a`-{pLPTD%=7UYJNvFuT$&a9`ih+`2TlWf)`r?-{D?gug}XQ32@ zly#|EU@y1RzPjV*AxjufbKXl1Ib1}IRUN<8HmOnAoZe20U)Ir=pjNl1Vcj&T6sy?; z&l>s0sjBgd8sA{iv9v!4ITEJJvG)^83a$Fa{5l)I?DO#R;GEQp^JvVEsE<^<9+Li8Afscg1kTZqcg8McbBsWh>Mpm=_mUq zDg@<4e*B)!J?JhL=QkI%7|||W87xtXLtcuy%PDDh4fS=6f=ptB$k6dt!>NXZwXfs+ zie?OqTc9LhcFr8n(VK?S5E-q!+RPcYtg6t4Z~N$OT|zoCzmJCzo| zT~pVAl16U2a{CPYezc`mo+rn+xUUJYfBj0(>KK{<$<-$7fU4+*?X(>UWjryykf&8T zSZJga_#*b?5@@7beF2ns#!4v?C_7Wuegn$Yfg0A=R%{tJxbXcE_*l9sVm6gu%{Qr3 z#b7IOF6aQMlHp4L#-en zdE-Xdyv-%OoTl^f@?Hai*wLUJ@g=VR;;VTr|RjAKN_jnG;Q-?yQ5V zRQKNPhY6B|8`i+ulJx<+3klg971@h`^WmXL3-xo`V;gl*BJ}Mb5c_R`qqlb06F!~m z_oQ9ZMQD(s97`?YWnGy%5E^?_| zr~kQDz4QH6@p;abXO18dvvi|IP{rhH6^llDKSZ_8eYmPp z`hAfq5Ni0(K|?`gjpy(E7H6V;&PHroE3&7aWgM_1-LC!8awD!BiZ#yh$`RiNP2=YC z?RrA(EN%HEo=aYW_c(amgpVqOhLdzz<<(2N=JMl%XfGq;$tjC=&yIGT%A$$W7Z!dT zUYt3tyTwf)TgFbx&1%J^T0A1OLZ)P!JTFvLmFtRZ9}-&fHV%C9DxEC|I%1&J zJPJBfb{wI&Co7BAuf8>GIWi!NjXbzK z{eZC5VpNWFa6#TCXi>iFPF3O^xCT8uhmhvtv}0_X^Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width"); fOnlySideBands = kTRUE; @@ -327,14 +320,10 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){ fRflFunc = CreateReflectionFunction("freflect"); fBkRFunc = CreateBackgroundPlusReflectionFunction("fbkgrfl"); } - // cout << "sfTemplates: " << fTemplates << endl; if(fTemplates){ printf(" ---> Final fit includes templates\n"); fTemplFunc = CreateTemplatesFunction("ftempl"); } - // for(int iFunc=0; iFuncGetNpar() << endl; for(Int_t ipar=0; iparGetParLimits(ipar,parmin,parmax); ftot->SetParLimits(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl,parmin,parmax); - cout << "Setting iPar" << ipar << " to " << fTemplFunc->GetParameter(ipar) - << ", parname: " << fTemplFunc->GetParName(ipar) << ", parlimits: [" - << parmin << "," << parmax << "]" << endl; ftot->SetParameter(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl,fTemplFunc->GetParameter(ipar)); ftot->SetParName(ipar+fNParsBkg+fNParsSig+fNParsSec+fNParsRfl,fTemplFunc->GetParName(ipar)); } } - cout << "Returning fTot" << endl; return ftot; } //__________________________________________________________________________ diff --git a/run3/flow/invmassfitter/VnVsMassFitter.cxx b/run3/flow/invmassfitter/VnVsMassFitter.cxx index 67261b08..cd87108b 100755 --- a/run3/flow/invmassfitter/VnVsMassFitter.cxx +++ b/run3/flow/invmassfitter/VnVsMassFitter.cxx @@ -225,55 +225,35 @@ VnVsMassFitter::~VnVsMassFitter() { //________________________________________________________________ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { - cout << "PERFORMING SIMULTANEOUS FIT" << endl; if(!fMassHisto || !fVnVsMassHisto) {printf("Histograms not set! Exit."); return kFALSE;} DefineNumberOfParameters(); const Int_t nparsmass = fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls; - cout << "nParsMass " << nparsmass << endl; Int_t NvnParsSgn = 1; if(fSecondPeak && fDoSecondPeakVn) {NvnParsSgn+=1;} if(fReflections && fVnRflOpt==kFreePar) {NvnParsSgn+=1;} - // if(fTemplates) {NvnParsSgn+=fNParsTempls;} - cout << "Added " << fNParsTempls << " to vn sgn pars" << endl; const Int_t nparsvn = nparsmass+fNParsVnBkg+NvnParsSgn; - cout << "nparsvn " << nparsvn << endl; - cout << "Starting PREFITS" << endl; Bool_t massprefit=MassPrefit(); if(!massprefit) {printf("Impossible to perform the mass prefit"); return kFALSE;} - cout << "Mass PREFIT" << endl; Bool_t vnprefit=VnSBPrefit(); - cout << "vnprefit: " << vnprefit << endl; if(!vnprefit) {printf("Impossible to perform the bkg vn prefit"); return kFALSE;} - cout << "Vn PREFIT" << endl; - cout << "PREFITS DONE" << endl; - - cout << "Mass function obtained from the prefit" << endl; - cout << "Npar: " << fMassFuncFromPrefit->GetNpar() << endl; std::vector initpars; for(Int_t iBkgPar=0; iBkgParGetParName(iBkgPar) << ", value: " << fMassFuncFromPrefit->GetParameter(iBkgPar) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iBkgPar)); } for(Int_t iSgnPar=0; iSgnParGetParName(iSgnPar+fNParsMassBkg) << ", value: " << fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg)); } for(Int_t iSecPeakPar=0; iSecPeakParGetParName(iSecPeakPar+fNParsMassBkg+fNParsMassSgn) << ", value: " << fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn)); } for(Int_t iReflPar=0; iReflParGetParName(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec) << ", value: " << fMassFuncFromPrefit->GetParameter(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iReflPar+fNParsMassBkg+fNParsMassSgn+fNParsSec)); } for(Int_t iTemplPar=0; iTemplParGetParName(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl) << ", value: " << fMassFuncFromPrefit->GetParameter(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl) << endl; initpars.push_back(fMassFuncFromPrefit->GetParameter(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl)); } - cout << "------------" << endl; - cout << "Background parameters for the VnBkgFunc: " << fNParsVnBkg << endl; for(Int_t iVnBkgPar=0; iVnBkgParGetParameter(iVnBkgPar));} else {initpars.push_back(0.05);} @@ -282,16 +262,9 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { if(fSecondPeak && fDoSecondPeakVn) {initpars.push_back(0.10);} //initial parameter for second peak vn fMassTotFunc = new TF1("fMassTotFunc",this,&VnVsMassFitter::MassFunc,fMassMin,fMassMax,nparsmass,"VnVsMassFitter","MassFunc"); - cout << "Number of parameters for vn function: " << nparsvn << endl; fVnTotFunc = new TF1("fVnTotFunc",this,&VnVsMassFitter::vnFunc,fMassMin,fMassMax,nparsvn,"VnVsMassFitter","vnFunc"); SetParNames(); - for(int iMassPar=0; iMassParGetNpar(); iMassPar++) { - cout << "[Mass] iPar" << iMassPar << ", name: " << fMassTotFunc->GetParName(iMassPar) << endl; - } - for(int iVnPar=0; iVnParGetNpar(); iVnPar++) { - cout << "[Vn] iPar" << iVnPar << ", name: " << fVnTotFunc->GetParName(iVnPar) << endl; - } ROOT::Math::WrappedMultiTF1 wfTotMass(*fMassTotFunc,1); ROOT::Math::WrappedMultiTF1 wfTotVn(*fVnTotFunc,1); @@ -332,33 +305,22 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { } if(fTemplates) { for(int iTemplPar=0; iTemplParfInitWeights.size(); iTemplPar++) { - cout << "Setting template parameter " << iTemplPar << endl; - cout << "Mass func parameter " << iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl << endl; - cout << fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Name() << endl; - cout << "[" << fInitWeights[iTemplPar] << ", " << fWeightsLowerLims[iTemplPar] - << ", " << fWeightsUpperLims[iTemplPar] << "]" << endl; - if(this->fWeightsLowerLims[iTemplPar] > this->fWeightsUpperLims[iTemplPar]) { - cout << "Fix mass" << endl; fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).Fix(); } else { - cout << "Vary mass" << endl; fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetValue(fInitWeights[iTemplPar]); fitter.Config().ParSettings(iTemplPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl).SetLimits(fWeightsLowerLims[iTemplPar],fWeightsUpperLims[iTemplPar]); } - cout << "---" << endl; } } fitter.Config().MinimizerOptions().SetPrintLevel(0); fitter.Config().SetMinimizer("Minuit2","Migrad"); for(Int_t iPar=0; iParGetParName(iPar));} - cout << "nparsvn: " << nparsvn << endl; // fit FCN function directly // (specify optionally data size and flag to indicate that is a chi2 fit Bool_t isFitOk = fitter.FitFCN(nparsvn,globalChi2,0,dataMass.Size()+dataVn.Size(),kFALSE); - cout << "isFitOk: " << isFitOk << endl; if(!isFitOk) return kFALSE; ROOT::Fit::FitResult result = fitter.Result(); @@ -379,7 +341,6 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { fMassTotFunc->SetParError(iPar,result.ParError(iPar)); } if(iPar>=nparsmass && iParSetParameter(iPar-nparsmass,result.Parameter(iPar)); fVnBkgFunc->SetParError(iPar-nparsmass,result.ParError(iPar)); } @@ -418,18 +379,18 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { }, fMassMin, fMassMax, 0)); fKDEVnTemplatesDraw.push_back(new TF1(Form("vnTempl_%s", fKDETemplates[iTempl].GetName()), - [&, this, iTempl, result, vnSgn] (double *x, double *par) { + [&, this, iTempl, vnSgn] (double *x, double *par) { return (vnSgn * this->fKDEMassTemplatesDraw[iTempl]->Eval(x[0])) / (this->fMassTotFunc->Eval(x[0])); }, fMassMin, fMassMax, 0)); } fKDEVnTemplatesDraw.push_back(new TF1("vnBkg", - [&, this, result, vnSgn] (double *x, double *par) { + [&, this, vnSgn] (double *x, double *par) { return (this->fVnBkgFunc->Eval(x[0]) * this->fMassBkgFunc->Eval(x[0])) / (this->fMassTotFunc->Eval(x[0])); }, fMassMin, fMassMax, 0)); fKDEVnTemplatesDraw.push_back(new TF1("vnSgn", - [&, this, result, vnSgn] (double *x, double *par) { + [&, this, vnSgn] (double *x, double *par) { return (vnSgn * this->fMassSgnFunc->Eval(x[0])) / (this->fMassTotFunc->Eval(x[0])); }, fMassMin, fMassMax, 0)); } @@ -451,7 +412,6 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { fChiSquare = result.MinFcnValue(); fNDF = result.Ndf(); fProb = result.Prob(); - cout << "Combined fit completed" << endl; return kTRUE; } @@ -459,7 +419,6 @@ Bool_t VnVsMassFitter::SimultaneousFit(Bool_t drawFit) { void VnVsMassFitter::DrawHere(TVirtualPad* c){ /// Core method to draw the fit output - cout << "[VnVsMassFitter] DrawHere" << endl; gStyle->SetOptStat(0); gStyle->SetCanvasColor(0); gStyle->SetFrameFillColor(0); @@ -508,8 +467,7 @@ void VnVsMassFitter::DrawHere(TVirtualPad* c){ fMassTotFunc->Draw("same"); } if(fTemplates) { - cout << "Drawing templates for mass spectrum" << endl; - for(int iMassTempl=0; iMassTemplfKDETemplates.size()/2; iMassTempl++) { + for(int iMassTempl=0; iMassTemplfKDEMassTemplatesDraw.size(); iMassTempl++) { fKDEMassTemplatesDraw[iMassTempl]->SetLineColor(kMagenta); fKDEMassTemplatesDraw[iMassTempl]->SetRange(fMassMin,fMassMax); fKDEMassTemplatesDraw[iMassTempl]->Draw("same"); @@ -555,8 +513,7 @@ void VnVsMassFitter::DrawHere(TVirtualPad* c){ fVnTotFunc->Draw("same"); } if(fTemplates) { - cout << "Drawing templates for vn" << endl; - for(int iVnTempl=this->fKDETemplates.size()/2; iVnTemplfKDETemplates.size(); iVnTempl++) { + for(int iVnTempl=0; iVnTemplfKDEVnTemplatesDraw.size(); iVnTempl++) { fKDEVnTemplatesDraw[iVnTempl]->SetLineColor(kMagenta); fKDEVnTemplatesDraw[iVnTempl]->SetRange(fMassMin,fMassMax); fKDEVnTemplatesDraw[iVnTempl]->Draw("same"); @@ -608,9 +565,7 @@ Bool_t VnVsMassFitter::MassPrefit() { if(fFixRflOverSig) {fMassFitter->SetFixReflOverS(fRflOverSig);} } if(fTemplates) {fMassFitter->SetTemplates(fKDETemplates, fInitWeights, fWeightsLowerLims, fWeightsUpperLims);} - cout << "Fitting" << endl; Bool_t status = fMassFitter->MassFitter(kFALSE); - cout << "Status " << status << endl; if(status) { fMassFuncFromPrefit = (TF1*)fMassFitter->GetMassFunc(); @@ -647,7 +602,6 @@ Bool_t VnVsMassFitter::VnSBPrefit() { TGraphErrors* gVnVsMassSB = new TGraphErrors(nSBbins); for(Int_t iBin=0; iBinGetBinContent(iBin+1) << endl; gVnVsMassSB->SetPoint(iBin,fVnVsMassHisto->GetBinCenter(iBin+1),fVnVsMassHisto->GetBinContent(iBin+1)); gVnVsMassSB->SetPointError(iBin,fVnVsMassHisto->GetBinWidth(iBin+1)/2,fVnVsMassHisto->GetBinError(iBin+1)); } @@ -827,11 +781,8 @@ void VnVsMassFitter::SetParNames() { if(fReflections) {fVnTotFunc->SetParName(fNParsMassBkg+fNParsMassSgn+fNParsSec,"ReflOverS");} if(fTemplates) { - cout << "Setting names of template parameters" << endl; for(int iTempl=0; iTemplfNParsTempls; iTempl++) { - cout << "iPar" << iTempl+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsVnBkg+fNParsVnSgn << ", name:" << Form("w_%s", fKDETemplates[iTempl].GetName()) << endl; fVnTotFunc->SetParName(iTempl+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl,Form("wm_%s", fKDETemplates[iTempl].GetName())); - fVnTotFunc->SetParName(iTempl+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn,Form("wf_%s", fKDETemplates[iTempl].GetName())); } } @@ -1075,17 +1026,6 @@ Double_t VnVsMassFitter::MassTemplates(Double_t *m,Double_t *pars){ return totalTemplates; } -//_________________________________________________________________________ -Double_t VnVsMassFitter::VnTemplates(Double_t *m,Double_t *pars){ - // Add the contributions of the templates loaded in fKDETemplates, each - // scaled by a multiplicative constant, left as free fit parameter - Double_t totalTemplates = 0.; - for(int iTempl=0; iTempl