diff --git a/run3/flow/WriteTF1.root b/run3/flow/WriteTF1.root
deleted file mode 100644
index 8725abbe..00000000
Binary files a/run3/flow/WriteTF1.root and /dev/null differ
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"):