Skip to content

Commit

Permalink
Removed leftovers
Browse files Browse the repository at this point in the history
  • Loading branch information
@Marcellocosti
Marcellocosti committed Jan 5, 2025
1 parent 5d26dbd commit 1d9f35b
Show file tree
Hide file tree
Showing 9 changed files with 66 additions and 192 deletions.
Binary file removed run3/flow/WriteTF1.root
View file
Binary file not shown.
9 changes: 1 addition & 8 deletions run3/flow/flow_analysis_utils.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -2,11 +2,11 @@
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):
Expand All @@ -29,9 +29,6 @@ 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()
Expand Down Expand Up @@ -73,8 +70,6 @@ def get_occupancy(thnSparse, occupancy_axis, debug=False):
- 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:
Expand All @@ -100,8 +95,6 @@ def get_evselbits(thnSparse, evselbits_axis, debug=False):
- 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:
Expand Down
82 changes: 34 additions & 48 deletions run3/flow/get_vn_vs_mass.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -19,7 +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
from kde_producer import kde_producer_grid, kde_producer_sim

def get_vn_vs_mass(fitConfigFileName, centClass, inFileName,
outputdir, suffix, vn_method, batch):
Expand Down Expand Up @@ -120,30 +120,26 @@ 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:

KDEtemplates = [[None]*len(fitConfig['TemplsFlags']) for _ in range(len(ptMins))]
if fitConfig.get('IncludeKDETempls'):
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], fitConfig['MCTemplatesFlags'][iPt][iFlag])
KDEtemplatesNames[iPt][iFlag] = fitConfig['MCTemplatesTreeNames'][iPt][iFlag]
for iFlag, flag in enumerate(fitConfig['TemplsFlags']):
if fitConfig.get('FromGrid'):
KDEtemplates[iPt][iFlag], hist_templ = kde_producer_grid(fitConfig['TemplsInputs'][iFlag],
'fM', ptMins[iPt], ptMaxs[iPt], flag)
elif fitConfig.get('FromSim'):
KDEtemplates[iPt][iFlag] = kde_producer_sim(fitConfig['TemplsInputs'][iFlag],
fitConfig['TemplsTreeNames'][iFlag],
ptMins[iPt], ptMaxs[iPt], fitConfig['TemplsFlags'][iFlag])
elif fitConfig.get('FromFile'):
templFile = TFile.Open(f'{fitConfig["FromFile"]}', 'r')
KDEtemplates[iPt][iFlag] = templFile.Get(f'KDE_pt_{ptMins[iPt]}_{ptMaxs[iPt]}_flag{flag}')
templFile.Close()
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)
print(f'ERROR: incorrect setting for including KDEs in fit! Exit!')
sys.exit()
KDEtemplatesFuncts = [[KDE.GetFunction() for KDE in KDEtemplatesPt] for KDEtemplatesPt in KDEtemplates]
print(KDEtemplatesFuncts)

# set particle configuration
if particleName == 'Ds':
Expand All @@ -162,8 +158,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))]
fMassTemplFuncts = [[None]*len(fitConfig['TemplsFlags']) for _ in range(len(ptMins))]
fVnTemplFuncts = [[None]*len(fitConfig['TemplsFlags']) for _ in range(len(ptMins))]
hist_reso = infile.Get('hist_reso')
hist_reso.SetDirectory(0)
reso = hist_reso.GetBinContent(1)
Expand Down Expand Up @@ -393,8 +389,8 @@ 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],
if fitConfig.get('IncludeKDETempls'):
vnFitter[iPt].SetKDETemplates(KDEtemplatesFuncts[iPt], fitConfig['TemplsTreeNames'],
fitConfig['InitWeights'][iPt], fitConfig['MinWeights'][iPt],
fitConfig['MaxWeights'][iPt])

Expand Down Expand Up @@ -476,10 +472,8 @@ 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')
for iMassTemplFunct, massTemplFunct in enumerate(fMassTemplFuncts[iPt]):
SetObjectStyle(massTemplFunct, color=kMagenta+iMassTemplFunct*2, linewidth=3)
print(f'Eval: {massTemplFunct.Eval(1.8)}')
massTemplFunct.SetLineColor(1)
massTemplFunct.Draw('same')
cSimFit[iCanv].Modified()
Expand All @@ -500,19 +494,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, vnTemplFunct in enumerate(fVnTemplFuncts[iPt]):
print(f'Eval: {vnTemplFunct.Eval(1.8)}')
file = TFile.Open("WriteTF1.root", "recreate")
file.cd()
vnTemplFunct.Write()
file.Close()
cSimFit[iPt].cd(2)
vnTemplFunct.SetLineColor(iVnTemplFunct+1)
vnTemplFunct.Draw('same')
cSimFit[iCanv].Modified()
cSimFit[iCanv].Update()

if fitConfig.get('drawvncomps'):
for iVnTemplFunct, vnTemplFunct in enumerate(fVnTemplFuncts[iPt]):
vnTemplFunct.SetLineColor(iVnTemplFunct+1)
vnTemplFunct.Draw('same')
cSimFit[iCanv].Modified()
cSimFit[iCanv].Update()
cSimFit[iCanv].Modified()
cSimFit[iCanv].Update()

Expand Down Expand Up @@ -743,8 +731,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}')
Expand Down Expand Up @@ -775,10 +763,8 @@ 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]}')
KDEtemplatesFuncts[ipt][iFlag].Write(f'{fitConfig["TemplsTreeNames"][iFlag]}_pt{ptmin*10:.0f}_{ptmax*10:.0f}_flag{fitConfig["TemplsFlags"][iFlag]}')

hSigmaSimFit.Write()
hMeanSimFit.Write()
Expand Down Expand Up @@ -834,8 +820,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')
Expand Down
18 changes: 0 additions & 18 deletions run3/flow/invmassfitter/InvMassFitter.cxx
View file
Original file line number Diff line number Diff line change
Expand Up @@ -261,13 +261,6 @@ Int_t InvMassFitter::MassFitter(Bool_t draw){
/// returns 2 if there is no signal and the fit is performed with only background
TVirtualFitter::SetDefaultFitter("Minuit");

for(int iFunc=0; iFunc<fTemplatesFuncts.size(); iFunc++) {
cout << "[InvMassFitter] Evaluating template " << iFunc << ": " << fTemplatesFuncts[iFunc].Eval(2) << endl;
cout << "[InvMassFitter] InitWeight " << fInitWeights[iFunc] << endl;
cout << "[InvMassFitter] MinWeight " << fWeightsLowerLims[iFunc] << endl;
cout << "[InvMassFitter] MaxWeight " << fWeightsUpperLims[iFunc] << endl;
}

Double_t integralHisto=fHistoInvMass->Integral(fHistoInvMass->FindBin(fMinMass),fHistoInvMass->FindBin(fMaxMass),"width");

fOnlySideBands = kTRUE;
Expand Down Expand Up @@ -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; iFunc<fTemplatesFuncts.size(); iFunc++) {
// cout << "[InvMassFitter] Inserting template " << fTemplatesFuncts[iFunc].Eval(2.0) << endl;
// }
fTotFunc = CreateTotalFitFunction("funcmass");

if(doFinalFit){
Expand Down Expand Up @@ -631,21 +620,14 @@ TF1* InvMassFitter::CreateTotalFitFunction(TString fname){
}
}
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; ipar<fNParsTempls; ipar++){
Double_t parmin,parmax;
fTemplFunc->GetParLimits(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;
}
//__________________________________________________________________________
Expand Down
Loading

0 comments on commit 1d9f35b

Please sign in to comment.