diff --git a/MC/CustomGenerators/PWGHF/HF_Muon_GenParam.C b/MC/CustomGenerators/PWGHF/HF_Muon_GenParam.C
new file mode 100644
index 00000000..d748418c
--- /dev/null
+++ b/MC/CustomGenerators/PWGHF/HF_Muon_GenParam.C
@@ -0,0 +1,85 @@
+AliGenerator* GeneratorCustom(TString opt = "")
+{
+  // return one of the parameterized MUON generators according to the parameter opt
+  // if opt = "custom" it uses the external parameterization macro MuonGeneratorCustom.C
+  // otherwise opt must match one of the options listed below
+  // the parameterization macro must implement the function AliGenerator* CreateGenerator();
+
+  // lists of available options and corresponding macros (to be updated each time a new macro is committed)
+  // the third list is to activate whatever special setup this option might needed. Currently implemented:
+  // - EvtGen --> load the libraries needed to use EvtGen (needed for J/psi radiative decay, ...)
+  const Int_t nParametrizations = 4;
+  TString optNames[nParametrizations] = {
+    "HF_single_pPbcmsl16r",
+    "HF_single_Pbpcmsl16s",
+    "HF_single_pPbcmsh16r",
+    "HF_single_Pbpcmsh16s",
+//     "jpsi_PbPb5TeV_1",
+//     "Charmonia_pp13TeV_1"
+  };
+  TString optMacros[nParametrizations] = {
+    "HF_Muon_GenParamSingle_pPb_CMSL_16r.C",
+    "HF_Muon_GenParamSingle_Pbp_CMSL_16s.C",
+    "HF_Muon_GenParamSingle_pPb_CMSH_16r.C",
+    "HF_Muon_GenParamSingle_Pbp_CMSH_16s.C",
+//     "HF_Muon_GenParamJpsi_PbPb5TeV_1.C",
+//     "HF_Muon_GenParamCharmonia_pp13TeV_1.C"
+  };
+  TString optSetups[nParametrizations] = {
+    "",
+    "",
+    "",
+    ""
+  };
+
+  TString macro = "";
+  TString setup = "";
+  if (opt == "custom") {
+    macro += "MuonGeneratorCustom.C";
+  } else {
+    for (Int_t iPar = 0; iPar < nParametrizations; iPar++) {
+      if (opt == optNames[iPar]) {
+        macro += "$ALIDPG_ROOT/MC/CustomGenerators/PWGHF/";
+        macro += optMacros[iPar];
+        setup += optSetups[iPar];
+        break;
+      }
+    }
+  }
+
+  if(macro.IsNull()) {
+    Printf("=============================================");
+    Printf("Error: Parametrization %s not known -> Abort",opt.Data());
+    Printf("=============================================");
+    abort();
+    return NULL;
+  }
+
+  Printf("=============================================");
+  Printf("Info:  Selected parametrization %s",opt.Data());
+  Printf("=============================================");
+
+  if (setup.Contains("EvtGen")) {
+    // load libraries to use EvtGen
+    gSystem->Load("libPhotos");
+    gSystem->Load("libEvtGen");
+    gSystem->Load("libTEvtGen");
+  }
+
+  // set external decayer (needed for AliGenParam)
+  TVirtualMCDecayer* decayer = new AliDecayerPythia();
+  decayer->SetForceDecay(kAll);
+  decayer->Init();
+  gMC->SetExternalDecayer(decayer);
+
+  // compile the macro in the current directory (needed to use precompiled functions)
+  gSystem->Exec(TString::Format("ln -s %s MuonGenerator.C", macro.Data()));
+  gSystem->AddIncludePath("-I$ALICE_ROOT/include");
+  if (gROOT->LoadMacro("MuonGenerator.C+") != 0) {
+    Printf("ERROR: cannot find %s\n", macro.Data());
+    abort();
+    return NULL;
+  }
+
+  return reinterpret_cast<AliGenerator*>(gROOT->ProcessLineSync("CreateGenerator()"));
+}
diff --git a/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_Pbp_CMSH_16s.C b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_Pbp_CMSH_16s.C
new file mode 100644
index 00000000..51f9d0c2
--- /dev/null
+++ b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_Pbp_CMSH_16s.C
@@ -0,0 +1,126 @@
+/*
+ *  MuonGenerator.C
+ *  aliroot_dev
+ *
+ *  Created by philippe pillot on 05/03/13.
+ *  Copyright 2013 SUBATECH. All rights reserved.
+ *
+ */
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include "TRandom.h"
+#include "AliGenerator.h"
+#include "AliGenParam.h"
+#endif
+
+
+static Int_t IpMuon( TRandom *ran);
+static Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ );
+static Double_t YMuon( const Double_t *py,  const Double_t * /*dummy*/ );
+static Double_t V2Muon( const Double_t *pv, const Double_t * /*dummy*/ );
+
+
+//-------------------------------------------------------------------------
+AliGenerator* CreateGenerator()
+{
+  AliGenParam *singleMu = new AliGenParam(1,-1,PtMuon,YMuon,V2Muon,IpMuon);
+  singleMu->SetMomentumRange(0,1e6);
+  singleMu->SetPtRange(4.,40.);
+  singleMu->SetYRange(-4.5,-2.);
+  singleMu->SetPhiRange(0., 360.);
+  singleMu->SetForceDecay(kNoDecay);
+  singleMu->SetTrackingFlag(1);
+  return singleMu;
+}
+
+//-------------------------------------------------------------------------
+Int_t IpMuon(TRandom *ran)
+{
+  // muon composition
+  
+  if (ran->Rndm() < 0.5 )
+  {
+    return 13;
+  }
+  else
+  {
+    return -13;
+  }
+}
+
+//-------------------------------------------------------------------------
+Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ )
+{
+  // muon pT
+  
+  Double_t x=*px;
+  if( x < 1. ) return 0.;
+  Float_t p0,p1,p2,p3, p4, p5, p6, p7;
+  
+  
+  p0 = 108.349;
+  p1 = -0.836925;
+  p2 = -0.345687;
+  p3 = 0.559775;
+  p4 = 3.14749;
+  p5 = 514732;
+  p6 = -17413.1;
+  p7 = 10928.4;
+
+ // norm
+  Float_t norm = p0;
+//slope
+  Float_t slope = p1 * ( 1. - TMath::Exp( p2*x)) + p3;
+//double 
+  Float_t den = 1./TMath::Power(x,p4);
+// pol
+  Float_t pol = p5 + p6*x+p7*x*x;
+  
+  return norm * TMath::Exp( slope*x ) * den * pol;
+}
+
+//-------------------------------------------------------------------------
+Double_t YMuon( const Double_t *py, const Double_t */*dummy*/ )
+{
+  // muon y
+
+//   Float_t yswitch = VAR_GENPARAMCUSTOMSINGLE_Y_RAPIDITY_SWITCH;  
+  Double_t x = *py;
+//   x *= yswitch;
+//pol4 only valid in y= -4;-2.5
+  Float_t p0,p1,p2,p3, p4, p5;
+  
+  p0 = 4.01732e+06;
+  p1 = -233742;
+  p2 = -362585;
+  p3 = 8383.75;
+  p4 = -98.7377;
+  p5 = 1.24522;
+
+ 
+//   p6 = VAR_GENPARAMCUSTOMSINGLE_Y_P6; //
+//   p7 = VAR_GENPARAMCUSTOMSINGLE_Y_P7; //
+//   p8 = VAR_GENPARAMCUSTOMSINGLE_Y_P8; //
+//   return p0 + 
+//         p1* x + 
+//         p2* x*x + 
+//         p3* x*x*x + 
+//         p4* x*x*x*x+
+// 	x*x*x*x*x ;
+//         + 
+//         p5* x*x*x*x*x + 
+//         p6* x*x*x*x*x*x + 
+//         p7* x*x*x*x*x*x*x + 
+//         p8* x*x*x*x*x*x*x*x; 
+// return p0*TMath::Power(x,4) +  p1*TMath::Power(x,3) +  p2*TMath::Power(x,2) + p3*TMath::Power(x,1) + p4;
+return p5*TMath::Power(x,8) + p4*TMath::Power(x,6) +  p3*TMath::Power(x,4) +  p2*TMath::Power(x,2) + p1*TMath::Power(x,1) + p0;
+}
+
+//-------------------------------------------------------------------------
+Double_t V2Muon( const Double_t */*dummy*/, const Double_t */*dummy*/ )
+{
+  //muon v2
+  return 0.;
+}
diff --git a/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_Pbp_CMSL_16s.C b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_Pbp_CMSL_16s.C
new file mode 100644
index 00000000..ca0ffe69
--- /dev/null
+++ b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_Pbp_CMSL_16s.C
@@ -0,0 +1,111 @@
+/*
+ *  MuonGenerator.C
+ *  aliroot_dev
+ *
+ *  Created by philippe pillot on 05/03/13.
+ *  Copyright 2013 SUBATECH. All rights reserved.
+ *
+ */
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include "TRandom.h"
+#include "AliGenerator.h"
+#include "AliGenParam.h"
+#endif
+
+
+static Int_t IpMuon( TRandom *ran);
+static Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ );
+static Double_t YMuon( const Double_t *py,  const Double_t * /*dummy*/ );
+static Double_t V2Muon( const Double_t *pv, const Double_t * /*dummy*/ );
+
+
+//-------------------------------------------------------------------------
+AliGenerator* CreateGenerator()
+{
+  AliGenParam *singleMu = new AliGenParam(1,-1,PtMuon,YMuon,V2Muon,IpMuon);
+  singleMu->SetMomentumRange(0,1e6);
+  singleMu->SetPtRange(0.,40.);
+  singleMu->SetYRange(-4.5,-2.);
+  singleMu->SetPhiRange(0., 360.);
+  singleMu->SetForceDecay(kNoDecay);
+  singleMu->SetTrackingFlag(1);
+  return singleMu;
+}
+
+//-------------------------------------------------------------------------
+Int_t IpMuon(TRandom *ran)
+{
+  // muon composition
+  
+  if (ran->Rndm() < 0.5 )
+  {
+    return 13;
+  }
+  else
+  {
+    return -13;
+  }
+}
+
+//-------------------------------------------------------------------------
+Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ )
+{
+  // muon pT
+  
+  Double_t x=*px;
+  if( x < 1. ) return 0.;
+  Float_t p0,p1,p2,p3, p4, p5, p6, p7;
+  
+  
+   p0 = 108.349;
+   p1 = -0.836925;
+   p2 = -0.345687;
+   p3 = 0.559775;
+   p4 = 3.14749;
+   p5 = 514732;
+   p6 = -17413.1;
+   p7 = 10928.4;
+
+ // norm
+  Float_t norm = p0;
+//slope
+  Float_t slope = p1 * ( 1. - TMath::Exp( p2*x)) + p3;
+//double 
+  Float_t den = 1./TMath::Power(x,p4);
+// pol
+  Float_t pol = p5 + p6*x+p7*x*x;
+  
+  return norm * TMath::Exp( slope*x ) * den * pol;
+}
+
+//-------------------------------------------------------------------------
+Double_t YMuon( const Double_t *py, const Double_t */*dummy*/ )
+{
+  // muon y
+
+//   Float_t yswitch = VAR_GENPARAMCUSTOMSINGLE_Y_RAPIDITY_SWITCH;  
+  Double_t x = *py;
+//   x *= yswitch;
+//pol4 only valid in y= -4;-2.5
+  Float_t p0,p1,p2,p3, p4, p5;
+  
+  p0 = 748673;
+  p1 = -4.41249e+06;
+  p2 = -1.65965e+06;
+  p3 = 68690.5;
+  p4 = -2465.22;
+  p5 = 40.2871;
+ 
+
+return p5*TMath::Power(x,8) + p4*TMath::Power(x,6) +  p3*TMath::Power(x,4) +  p2*TMath::Power(x,2) + p1*TMath::Power(x,1) + p0;
+}
+
+//-------------------------------------------------------------------------
+Double_t V2Muon( const Double_t */*dummy*/, const Double_t */*dummy*/ )
+{
+  //muon v2
+  return 0.;
+}
diff --git a/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_pPb_CMSH_16r.C b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_pPb_CMSH_16r.C
new file mode 100644
index 00000000..9cf1134b
--- /dev/null
+++ b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_pPb_CMSH_16r.C
@@ -0,0 +1,125 @@
+/*
+ *  MuonGenerator.C
+ *  aliroot_dev
+ *
+ *  Created by philippe pillot on 05/03/13.
+ *  Copyright 2013 SUBATECH. All rights reserved.
+ *
+ */
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include "TRandom.h"
+#include "AliGenerator.h"
+#include "AliGenParam.h"
+#endif
+
+
+static Int_t IpMuon( TRandom *ran);
+static Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ );
+static Double_t YMuon( const Double_t *py,  const Double_t * /*dummy*/ );
+static Double_t V2Muon( const Double_t *pv, const Double_t * /*dummy*/ );
+
+
+//-------------------------------------------------------------------------
+AliGenerator* CreateGenerator()
+{
+  AliGenParam *singleMu = new AliGenParam(1,-1,PtMuon,YMuon,V2Muon,IpMuon);
+  singleMu->SetMomentumRange(0,1e6);
+  singleMu->SetPtRange(0.,40.);
+  singleMu->SetYRange(-4.5,-2.);
+  singleMu->SetPhiRange(0., 360.);
+  singleMu->SetForceDecay(kNoDecay);
+  singleMu->SetTrackingFlag(1);
+  return singleMu;
+}
+
+//-------------------------------------------------------------------------
+Int_t IpMuon(TRandom *ran)
+{
+  // muon composition
+  
+  if (ran->Rndm() < 0.5 )
+  {
+    return 13;
+  }
+  else
+  {
+    return -13;
+  }
+}
+
+//-------------------------------------------------------------------------
+Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ )
+{
+  // muon pT
+  
+  Double_t x=*px;
+  if( x < 1. ) return 0.;
+  Float_t p0,p1,p2,p3, p4, p5, p6, p7;
+  
+   p0 = 240.224;
+   p1 = -0.695575;
+   p2 = -0.281131;
+   p3 = 0.486148;
+   p4 = 3.08035;
+   p5 = 309864;
+   p6 = -17990.7;
+   p7 = 4423.42;
+
+ // norm
+  Float_t norm = p0;
+//slope
+  Float_t slope = p1 * ( 1. - TMath::Exp( p2*x)) + p3;
+//double 
+  Float_t den = 1./TMath::Power(x,p4);
+// pol
+  Float_t pol = p5 + p6*x+p7*x*x;
+  
+  return norm * TMath::Exp( slope*x ) * den * pol;
+}
+
+//-------------------------------------------------------------------------
+Double_t YMuon( const Double_t *py, const Double_t */*dummy*/ )
+{
+  // muon y
+
+//   Float_t yswitch = VAR_GENPARAMCUSTOMSINGLE_Y_RAPIDITY_SWITCH;  
+  Double_t x = *py;
+//   x *= yswitch;
+//pol4 only valid in y= -4;-2.5
+  Float_t p0,p1,p2,p3, p4, p5;
+  
+p0 = 749110;
+p1 = -4.41199e+06; 
+p2 = -1.65948e+06; 
+p3 = 68682.3;
+p4 = -2464.92;
+p5 = 40.2823;
+
+ 
+//   p6 = VAR_GENPARAMCUSTOMSINGLE_Y_P6; //
+//   p7 = VAR_GENPARAMCUSTOMSINGLE_Y_P7; //
+//   p8 = VAR_GENPARAMCUSTOMSINGLE_Y_P8; //
+//   return p0 + 
+//         p1* x + 
+//         p2* x*x + 
+//         p3* x*x*x + 
+//         p4* x*x*x*x+
+// 	x*x*x*x*x ;
+//         + 
+//         p5* x*x*x*x*x + 
+//         p6* x*x*x*x*x*x + 
+//         p7* x*x*x*x*x*x*x + 
+//         p8* x*x*x*x*x*x*x*x; 
+// return p0*TMath::Power(x,4) +  p1*TMath::Power(x,3) +  p2*TMath::Power(x,2) + p3*TMath::Power(x,1) + p4;
+return p5*TMath::Power(x,8) + p4*TMath::Power(x,6) +  p3*TMath::Power(x,4) +  p2*TMath::Power(x,2) + p1*TMath::Power(x,1) + p0;
+}
+
+//-------------------------------------------------------------------------
+Double_t V2Muon( const Double_t */*dummy*/, const Double_t */*dummy*/ )
+{
+  //muon v2
+  return 0.;
+}
diff --git a/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_pPb_CMSL_16r.C b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_pPb_CMSL_16r.C
new file mode 100644
index 00000000..44111469
--- /dev/null
+++ b/MC/CustomGenerators/PWGHF/HF_Muon_GenParamSingle_pPb_CMSL_16r.C
@@ -0,0 +1,125 @@
+/*
+ *  MuonGenerator.C
+ *  aliroot_dev
+ *
+ *  Created by philippe pillot on 05/03/13.
+ *  Copyright 2013 SUBATECH. All rights reserved.
+ *
+ */
+
+
+#if !defined(__CINT__) || defined(__MAKECINT__)
+#include <Riostream.h>
+#include "TRandom.h"
+#include "AliGenerator.h"
+#include "AliGenParam.h"
+#endif
+
+
+static Int_t IpMuon( TRandom *ran);
+static Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ );
+static Double_t YMuon( const Double_t *py,  const Double_t * /*dummy*/ );
+static Double_t V2Muon( const Double_t *pv, const Double_t * /*dummy*/ );
+
+
+//-------------------------------------------------------------------------
+AliGenerator* CreateGenerator()
+{
+  AliGenParam *singleMu = new AliGenParam(1,-1,PtMuon,YMuon,V2Muon,IpMuon);
+  singleMu->SetMomentumRange(0,1e6);
+  singleMu->SetPtRange(0.,40.);
+  singleMu->SetYRange(-4.5,-2.);
+  singleMu->SetPhiRange(0., 360.);
+  singleMu->SetForceDecay(kNoDecay);
+  singleMu->SetTrackingFlag(1);
+  return singleMu;
+}
+
+//-------------------------------------------------------------------------
+Int_t IpMuon(TRandom *ran)
+{
+  // muon composition
+  
+  if (ran->Rndm() < 0.5 )
+  {
+    return 13;
+  }
+  else
+  {
+    return -13;
+  }
+}
+
+//-------------------------------------------------------------------------
+Double_t PtMuon( const Double_t *px, const Double_t * /*dummy*/ )
+{
+  // muon pT
+  
+  Double_t x=*px;
+  if( x < 1. ) return 0.;
+  Float_t p0,p1,p2,p3, p4, p5, p6, p7;
+  
+  p0 = 240.224;
+  p1 = -0.695575;
+  p2 = -0.281131;
+  p3 = 0.486148;
+  p4 = 3.08035;
+  p5 = 309864;
+  p6 = -17990.7;
+  p7 = 4423.42;
+  
+ // norm
+  Float_t norm = p0;
+//slope
+  Float_t slope = p1 * ( 1. - TMath::Exp( p2*x)) + p3;
+//double 
+  Float_t den = 1./TMath::Power(x,p4);
+// pol
+  Float_t pol = p5 + p6*x+p7*x*x;
+  
+  return norm * TMath::Exp( slope*x ) * den * pol;
+}
+
+//-------------------------------------------------------------------------
+Double_t YMuon( const Double_t *py, const Double_t */*dummy*/ )
+{
+  // muon y
+
+//   Float_t yswitch = VAR_GENPARAMCUSTOMSINGLE_Y_RAPIDITY_SWITCH;  
+  Double_t x = *py;
+//   x *= yswitch;
+//pol4 only valid in y= -4;-2.5
+  Float_t p0,p1,p2,p3, p4, p5;
+  
+ p0 = 749110;
+ p1 = -4.41199e+06; 
+ p2 = -1.65948e+06; 
+ p3 = 68682.3;
+ p4 = -2464.92;
+ p5 = 40.2823;
+
+ 
+//   p6 = VAR_GENPARAMCUSTOMSINGLE_Y_P6; //
+//   p7 = VAR_GENPARAMCUSTOMSINGLE_Y_P7; //
+//   p8 = VAR_GENPARAMCUSTOMSINGLE_Y_P8; //
+//   return p0 + 
+//         p1* x + 
+//         p2* x*x + 
+//         p3* x*x*x + 
+//         p4* x*x*x*x+
+// 	x*x*x*x*x ;
+//         + 
+//         p5* x*x*x*x*x + 
+//         p6* x*x*x*x*x*x + 
+//         p7* x*x*x*x*x*x*x + 
+//         p8* x*x*x*x*x*x*x*x; 
+// return p0*TMath::Power(x,4) +  p1*TMath::Power(x,3) +  p2*TMath::Power(x,2) + p3*TMath::Power(x,1) + p4;
+return p5*TMath::Power(x,8) + p4*TMath::Power(x,6) +  p3*TMath::Power(x,4) +  p2*TMath::Power(x,2) + p1*TMath::Power(x,1) + p0;
+}
+
+//-------------------------------------------------------------------------
+Double_t V2Muon( const Double_t */*dummy*/, const Double_t */*dummy*/ )
+{
+  //muon v2
+  return 0.;
+}