From 1d9f35b87707299bb65b872462bab98f066fb246 Mon Sep 17 00:00:00 2001
From: marcellocosti <marcellodicostanzo00@gmail.com>
Date: Sat, 4 Jan 2025 17:06:16 +0100
Subject: [PATCH] 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*%<ca)f&qVtbb_=<w*ew|%$RNf0G{z<H&5q#
zrao?d&c=^H!aje0^KVB0-e1%m5&*EqfRr$GhXw%1KL58b;S&DsjLUn62i|pF`HRQF
zw8IPVdJ6XXzu7{Ua25x^iwg+C^0l)wbXoJZvrDtH!+peOC-j9xe=<5vMV1Voht0{!
z4!a5QO3)r3pDl?*aL9GFx7Y4=hm^3e<rYGUPf7t06TwR`FAbIPw{sxPy;s4kvh5QU
z(*Vsq>}}c$pdK`96>i~!T8frTAa5TRBL94rhF7eFH!g+ndnry9Vs}nY9QB0c2sd<|
z?ESoFjWp{7!zld^Ltu!ox@#dSvmwXH7nSopD~D~023t*dTXh<!!N8J*BT>+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`kJesm<VwlBwb0N{QPQ>M51662f-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;7K9<T8so!}I`?&gtR{xi07zcLc0@W0SUXgn;-lOO7u}
zp^ino2$K06Y?C?aLZn;&aoFgFMp4d=_n&tX)@hd+%p1|UyuPSki0+cw%fMr(%F>dg
zm*r=XL<9?jtmkd%0l4#Wb9o`*=p~GV$`q~u;_QQYAkhowXoscmF-q)%$zCub=@x0J
zasHMIqo222-TScGlcn)+G<TNFDJNfaBsBMH1RJ5ZIH{(f_PCyqmJkr%SpkIqnC*2s
z<-=eAy-c~~op@n%;KA9Q4ygRqfxnUlaKbazb@TNjed_1qZ2#QZN6X6vLZt@rckq14
z=IQ3;;p~K2eEl%$kVynMdw@}g3G)eGnqk@q`3x1!el~2fxVGJljaSAhFNmwg@cRha
z@fgzhnnPs<c%Yq%-+NKxVyv2p?B=|UXX}|?0@{Hlqug5&9I_ogW9+}UvQ$$>*-fl`
zw6;{juPebYFlQgJCU-@CUa{fi75oN1e+)lral6G~EG1?bZqUnGRTwa|*5%lfaF)2=
z#=d?CZl2&-Og=&ow0Tv3t<Z?RIX#fpfAkm+est1G(qj44aPw33v6rDOL%BiULt~35
z*^L1sI9w~N&895es>kRR<<uPT-B_?#VzG_2m5tel#I)7567Y-clvH#?_(14Xgp0d3
zzDOx2H@DrJX_9Q3^VX663!rKQB8KI~tc-@Z#JZ$YT2LJ?Cr5sN`~97`;b-bx0MBwr
zYqbyduLGy`;5_e`*U-V*;MXD=5>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?J<FhGfK60WIM)EAbaW
zbnxR5d&NTf=NwMM!!d8L=E;k2vr9>y!*GiZiMvWAC$ZGBHt{r@J<IPe<w1}3q;O8`
zr$@BTRk&DXzhY~UOmk%XG6CC`?Y@tXP+O(+V^B9;@N=?k71h_1n)>x^y|u0W32U9+
z_L?z!dZMwCp>=k79PR@W={GT)*ja1-#kBD=)AoHY-I!SQHUybFDEL<L%;%qjfn5_o
zCzl|EaaqC#d6}QGM5iu7dQMge;*4^;bPyk|FFYN{uW9>l&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^{<CaiW(0-Ie^)#uI6YbtkZ>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;G<t{w=t5z@se_m?#(7{
zxHqXZ`@fmO!s@FPZX>EQp^JvVEsE<^<9+Li8Afscg1kTZqcg8McbBsWh>Mpm=_mUq
zDg@<4e*B)!J?JhL=QkI%7|||W87xtXLtcuy%PDDh4fS=6f=ptB$k6dt!>NXZwXfs+
zie?OqTc9LhcFr8<l`9YX-6b{US>n(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#b7IOF6aQMlHp4<Yu_Jqe};d<7(dFUo8nqaqVXoH5B%bz8GOBwxt&7X)|6F>L#-en
zdE-Xdyv-%OoTl^f@?Hai*wLUJ@g=VR;;V<ZYxv2SdgQq>Tr|RjAKN_jnG;Q-?yQ5V
zRQKNPhY6B|8`i+ulJx<+3klg971@h`^WmXL3-xo`V;gl*BJ}Mb5c_R`qqlb06F!~m
z_oQ9ZMQD(<oD0HS0;+J4IkknO*go-69MZAj6*%^bWIV`_=d~wO@kJ50WaPemnp&BN
zf~m}ctt;(H1&mQ`yvRPK*m0sNKeYtJy1Ioh1X7$L8gEw=y|Gg2qjU+f=L2OqhXGyd
zCq~n5MKd9J@sx^WJ8yI1J18Wb-uYC1g!4Wx50GCvu;L~$pgb@k?pF`W8$oNx&(;zL
z6b?W&IW|86L;Uv?)l;=LE-fB>s97`?YWnGy<S*lv53D&Ep4@&9m*K!Rj_>%5E^?_|
zr~kQDz4QH6@p;abXO18dvvi|IP{rhH6<ey~&asyGrhA2^a)!JV>^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|<fH4Ja-joO
z<;7m}kZiqs6tN3cig6Nsf{0D+RdIIZ#9#Z##AJ!`lYG|FObz-(Y*JJgK-ajFH=zD$
zl)iIx0kIwGHvgia|5;DPfV|v@yDC4qvZBcj=r%V){`b2dK2bUZv`aJ7IdVO&2np22
zVN{6hM3-^<R8O1qSn1ber3istf~g4ohK8a6m`J!H;q(XLsxaOZ-r+R+VYMG>I%1&J
z<v|#Z(dFhsvS|ujvQl|{Wh+#HeSACunx)p>JPJBfb{wI&Co7BAuf8>GIWi!NjXbzK
z{eZC5VpNWFa6#TCXi>iFPF3O^xCT8uhmhvtv}0_X^<BIazB6}s@e&gZ@66qQh?kg)
jzwYmU%!QcOvUlXf|4)9GkTK5jKapdo$;|^LGXehs$*@p@

diff --git a/run3/flow/flow_analysis_utils.py b/run3/flow/flow_analysis_utils.py
index 4c63b836..a5c96655 100644
--- a/run3/flow/flow_analysis_utils.py
+++ b/run3/flow/flow_analysis_utils.py
@@ -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):
@@ -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()
@@ -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:
@@ -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:
diff --git a/run3/flow/get_vn_vs_mass.py b/run3/flow/get_vn_vs_mass.py
index e2d65e49..5a390dfc 100644
--- a/run3/flow/get_vn_vs_mass.py
+++ b/run3/flow/get_vn_vs_mass.py
@@ -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):
@@ -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':
@@ -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)
@@ -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])
 
@@ -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()
@@ -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()
 
@@ -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}')
@@ -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()
@@ -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')
diff --git a/run3/flow/invmassfitter/InvMassFitter.cxx b/run3/flow/invmassfitter/InvMassFitter.cxx
index 6158eb9f..72fedf10 100644
--- a/run3/flow/invmassfitter/InvMassFitter.cxx
+++ b/run3/flow/invmassfitter/InvMassFitter.cxx
@@ -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;
@@ -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){
@@ -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;
 }
 //__________________________________________________________________________
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<Double_t> initpars;
   for(Int_t iBkgPar=0; iBkgPar<fNParsMassBkg; iBkgPar++) {
-    cout << "Bkg par: " << fMassFuncFromPrefit->GetParName(iBkgPar) << ", value: " << fMassFuncFromPrefit->GetParameter(iBkgPar) << endl;
     initpars.push_back(fMassFuncFromPrefit->GetParameter(iBkgPar));
   }
   for(Int_t iSgnPar=0; iSgnPar<fNParsMassSgn; iSgnPar++) {
-    cout << "Sgn par: " << fMassFuncFromPrefit->GetParName(iSgnPar+fNParsMassBkg) << ", value: " << fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg) << endl;
     initpars.push_back(fMassFuncFromPrefit->GetParameter(iSgnPar+fNParsMassBkg));
   }
   for(Int_t iSecPeakPar=0; iSecPeakPar<fNParsSec; iSecPeakPar++) {
-    cout << "SecPeak par: " << fMassFuncFromPrefit->GetParName(iSecPeakPar+fNParsMassBkg+fNParsMassSgn) << ", value: " << fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn) << endl;
     initpars.push_back(fMassFuncFromPrefit->GetParameter(iSecPeakPar+fNParsMassBkg+fNParsMassSgn));
   }
   for(Int_t iReflPar=0; iReflPar<fNParsRfl; iReflPar++) {
-    cout << "Refl par: " << fMassFuncFromPrefit->GetParName(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; iTemplPar<fNParsTempls; iTemplPar++) {
-    cout << "Templ par: " << fMassFuncFromPrefit->GetParName(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; iVnBkgPar<fNParsVnBkg; iVnBkgPar++) {
     if(vnprefit) {initpars.push_back(fVnBkgFuncSb->GetParameter(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; iMassPar<fMassTotFunc->GetNpar(); iMassPar++) {
-    cout << "[Mass] iPar" << iMassPar << ", name: " << fMassTotFunc->GetParName(iMassPar) << endl;
-  }
-  for(int iVnPar=0; iVnPar<fVnTotFunc->GetNpar(); 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; iTemplPar<this->fInitWeights.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; iPar<nparsvn; iPar++) {fitter.Config().ParSettings(iPar).SetName(fVnTotFunc->GetParName(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 && iPar<nparsvn-NvnParsSgn) {
-      cout << "Setting parameter " << iPar-nparsmass << " of fVnBkgFunc to " << result.Parameter(iPar) << endl; 
       fVnBkgFunc->SetParameter(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; iMassTempl<this->fKDETemplates.size()/2; iMassTempl++) {
+    for(int iMassTempl=0; iMassTempl<this->fKDEMassTemplatesDraw.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; iVnTempl<this->fKDETemplates.size(); iVnTempl++) {
+    for(int iVnTempl=0; iVnTempl<this->fKDEVnTemplatesDraw.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; iBin<nMassBins; iBin++) {
     if(SBbins[iBin]==1) {
-      cout << "Setting point " << iBin << " to " << fVnVsMassHisto->GetBinContent(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; iTempl<this->fNParsTempls; 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<fNParsTempls; iTempl++) {
-    totalTemplates += pars[iTempl]*fKDETemplates[iTempl].Eval(m[0]);
-  }
-  return totalTemplates;
-}
-
 //_________________________________________________________________________
 Double_t VnVsMassFitter::MassBkgRfl(Double_t *m,Double_t *pars){
 
@@ -1147,7 +1087,6 @@ Double_t VnVsMassFitter::MassFunc(Double_t *m, Double_t *pars) {
   Double_t rflpars[1]; //maximum number of parameters for rfl = 1 for the implemented functions
   for(Int_t iPar=0; iPar<fNParsRfl; iPar++) {rflpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn+fNParsSec];}
   //templates parameters
-  // cout << "Number of parameters of templates" << fNParsTempls << endl;
   Double_t templpars[fNParsTempls]; //one parameter for each template (i.e. its scaling parameter)
   for(Int_t iPar=0; iPar<fNParsTempls; iPar++) { templpars[iPar] = pars[iPar+fNParsMassBkg+fNParsMassSgn+fNParsSec+fNParsRfl];}
 
@@ -1163,7 +1102,6 @@ Double_t VnVsMassFitter::MassFunc(Double_t *m, Double_t *pars) {
 Double_t VnVsMassFitter::vnFunc(Double_t *m, Double_t *pars) {
 
   //bkg mass parameters
-  // cout << "Background parameters: " << fNParsMassBkg << endl;
   const Int_t nBkgPars = fNParsMassBkg;
   Double_t massbkgpars[nBkgPars];
   for(Int_t iPar=0; iPar<fNParsMassBkg; iPar++) {massbkgpars[iPar] = pars[iPar];}
@@ -1185,7 +1123,6 @@ Double_t VnVsMassFitter::vnFunc(Double_t *m, Double_t *pars) {
   Double_t vnbkgpars[nVnBkgPars];
   for(Int_t iPar=0; iPar<fNParsVnBkg; iPar++) {
     vnbkgpars[iPar] = pars[iPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls];
-    // cout << "vnbkgpars iPar: " << iPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls << endl;
   }
   //signal vn parameter
   Double_t vnSgn = pars[fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg];
@@ -1213,9 +1150,6 @@ Double_t VnVsMassFitter::vnFunc(Double_t *m, Double_t *pars) {
         break;
     }
   }
-  // Double_t vntemplpars[fNParsTempls];
-  // for(Int_t iPar=0; iPar<fNParsTempls; iPar++) {vntemplpars[iPar] = pars[iPar+fNParsMassSgn+fNParsMassBkg+fNParsSec+fNParsRfl+fNParsTempls+fNParsVnBkg+fNParsVnSgn];}
-  
   Double_t vnBkg = vnBkgFunc(m,vnbkgpars);
   Double_t Sgn = MassSignal(m,masssgnpars);
   Double_t Bkg = MassBkg(m,massbkgpars);
@@ -1230,8 +1164,7 @@ Double_t VnVsMassFitter::vnFunc(Double_t *m, Double_t *pars) {
     else Refl += MassRfl(m,rflpars);
   }
   Double_t Templates = 0;
-  // if(fTemplates) {Templates += VnTemplates(m,vntemplpars);}
-  if(fTemplates) {Templates += VnTemplates(m,templpars);}
+  if(fTemplates) {Templates += MassTemplates(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 68137ba7..a19b8186 100755
--- a/run3/flow/invmassfitter/VnVsMassFitter.h
+++ b/run3/flow/invmassfitter/VnVsMassFitter.h
@@ -186,12 +186,11 @@ 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 VnTemplates(Double_t *m,Double_t *pars);
   Double_t MassBkgRfl(Double_t *m,Double_t *par);
+  Double_t MassTemplates(Double_t *m,Double_t *pars);
   Double_t MassSecondPeak(Double_t *m,Double_t *par);
-  Double_t vnBkgFunc(Double_t *m, Double_t *pars);
   Double_t MassFunc(Double_t *m, Double_t *pars);
+  Double_t vnBkgFunc(Double_t *m, Double_t *pars);
   Double_t vnFunc(Double_t *m, Double_t *pars);
 
     ///private methods
diff --git a/run3/flow/kde_producer.py b/run3/flow/kde_producer.py
index 2f560683..8dd85388 100644
--- a/run3/flow/kde_producer.py
+++ b/run3/flow/kde_producer.py
@@ -7,23 +7,20 @@
 import ctypes
 from ROOT import TFile, TKDE, TCanvas, TH1D
 
-def kde_producer(tree_file, var, pt_min, pt_max, inv_mass_bins, flag, outfile=''):
+def kde_producer_grid(tree_file, var, pt_min, pt_max, flag, outfile=''):
     
     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 = []
     with uproot.open(f'{tree_file}') as f:
-        # print(f.keys())
-        for iKey, key in enumerate(f.keys()):
+        for key in 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)
     pt_filtered_df = full_dataset.query(f"{pt_min} <= fPt < {pt_max}")
     filtered_df = pt_filtered_df.query(f"fFlagMcMatchRec == {flag} or fFlagMcMatchRec == {-flag}")
     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)
@@ -36,7 +33,7 @@ def kde_producer(tree_file, var, pt_min, pt_max, inv_mass_bins, flag, outfile=''
         cOverlap = TCanvas('cOverlap', 'cOverlap', 600, 600)
         cOverlap.cd()
         binned_var_values.Draw()
-        kde_func.Draw('SAME')
+        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()
@@ -45,32 +42,25 @@ def kde_producer(tree_file, var, pt_min, pt_max, inv_mass_bins, flag, outfile=''
     
     return kde, binned_var_values 
 
-def kde_producer_sim(tree_file, tree_name, pt_min, pt_max, flag=0, outfile=''):
+def kde_producer_sim(tree_file, tree_name, pt_min, pt_max, flag=0):
     
     print(f"Producing KDE from sim {tree_file} for process {tree_name}, {pt_min} <= pt < {pt_max}, flag {flag}")
     # 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()):
+        for key in 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)
     if flag!=0:
-        # print('Non trivial flag')
-        # print(flag)
         flagged_dataset = full_dataset.query(f"decayflag == {flag}")
     else:
         flagged_dataset = full_dataset
-    # print(flagged_dataset)
     pt_filtered_df = flagged_dataset.query(f"{pt_min} <= pt < {pt_max}")
     var_values = pt_filtered_df["mass"].tolist()  # Or use `.tolist()` to get a list
-    # print(full_dataset)
-    # print(var_values)
-    
+
     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)
@@ -81,29 +71,24 @@ def kde_producer_sim(tree_file, tree_name, pt_min, pt_max, flag=0, outfile=''):
 
 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("--var", "-v", metavar="text",
+                        default="fM", help="variable of interest")
     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 file containing 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 = []
@@ -116,13 +101,10 @@ def kde_producer_sim(tree_file, tree_name, pt_min, pt_max, flag=0, outfile=''):
     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'], \
-                                      pt_low, pt_max, inv_mass_bins, config['chn_flag'],
-                                      outfile)
+        for pt_low, pt_max in zip(config['pt_mins'], config['pt_maxs']):
+            KDE, histo = kde_producer_grid(config['input'], config['variable'], pt_low, 
+                                           pt_max, config['chn_flag'], outfile)
     else:
-        KDE, histo = kde_producer(args.input, args.var, \
-                                  args.ptmin, args.ptmax, args.invmassbins, args.flag,
-                                  outfile)
+        KDE, histo = kde_producer_grid(args.input, args.var, args.ptmin,
+                                       args.ptmax, args.flag, outfile)
     outfile.Close()    
-    
\ No newline at end of file
diff --git a/run3/flow/project_thnsparse.py b/run3/flow/project_thnsparse.py
index bba65b51..5f2876cf 100644
--- a/run3/flow/project_thnsparse.py
+++ b/run3/flow/project_thnsparse.py
@@ -22,8 +22,6 @@ def check_anres(config, an_res_file, centrality, resolution,
     axis_mass = config['axes']['mass']
     axis_sp = config['axes']['sp']
     axis_deltaphi = config['axes']['deltaphi']
-    # axis_occupancy = config['axes']['occupancy']
-    # axis_evselbits = config['axes']['evselbits']
     inv_mass_bins = config['inv_mass_bins']
     use_inv_mass_bins = config['use_inv_mass_bins']
     if use_inv_mass_bins:
@@ -137,15 +135,17 @@ def check_anres(config, an_res_file, centrality, resolution,
                 hist_mass.SetName(f'hist_mass_cent{cent_min}_{cent_max}_pt{pt_min}_{pt_max}')
                 hist_mass.Write()
             
-            # hist_occ = get_occupancy(thnsparse_selcent, axis_occupancy, False)
-            # hist_occ.SetName(f'hist_occ_pt{pt_min}_{pt_max}')
-            # hist_occ.SetDirectory(0)
-            # hist_occ.Write()
+            if config['axes'].get('occupancy'):
+                hist_occ = get_occupancy(thnsparse_selcent, config['axes']['occupancy'], False)
+                hist_occ.SetName(f'hist_occ_pt{pt_min}_{pt_max}')
+                hist_occ.SetDirectory(0)
+                hist_occ.Write()
 
-            # hist_evselbits = get_evselbits(thnsparse_selcent, axis_evselbits, False)
-            # hist_evselbits.SetName(f'hist_evselbits_pt{pt_min}_{pt_max}')
-            # hist_evselbits.SetDirectory(0)
-            # hist_evselbits.Write()
+            if config['axes'].get('evselbits'):
+                hist_evselbits = get_evselbits(thnsparse_selcent, config['axes']['evselbits'], False)
+                hist_evselbits.SetName(f'hist_evselbits_pt{pt_min}_{pt_max}')
+                hist_evselbits.SetDirectory(0)
+                hist_evselbits.Write()
             
             bar()
             outfile.cd('..')
diff --git a/run3/flow/run_full_flow_analysis.py b/run3/flow/run_full_flow_analysis.py
index 7d0d6268..97344813 100644
--- a/run3/flow/run_full_flow_analysis.py
+++ b/run3/flow/run_full_flow_analysis.py
@@ -115,7 +115,6 @@ def run_full_analysis(config,
         print("\n\033[92m Starting vn extraction\033[0m")
         print(f"\033[92m {command_vn}\033[0m")
         os.system(command_vn)
-    exit()
 
     # copy config file
     if not os.path.exists(f"{outputdir}/config"):