From a36c494080816764ef9f36d57e0c6d9b437cece3 Mon Sep 17 00:00:00 2001
From: ppillot <pillot@subatech.in2p3.fr>
Date: Wed, 21 Nov 2018 14:31:49 +0100
Subject: [PATCH] new single muon parameterisation for PbPb at 5TeV

---
 MC/CustomGenerators/PWGDQ/Muon_GenParam.C     |  5 +-
 .../PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C    | 88 +++++++++++++++++++
 2 files changed, 92 insertions(+), 1 deletion(-)
 create mode 100644 MC/CustomGenerators/PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C

diff --git a/MC/CustomGenerators/PWGDQ/Muon_GenParam.C b/MC/CustomGenerators/PWGDQ/Muon_GenParam.C
index 78c2f3ad..a986e31c 100644
--- a/MC/CustomGenerators/PWGDQ/Muon_GenParam.C
+++ b/MC/CustomGenerators/PWGDQ/Muon_GenParam.C
@@ -8,18 +8,21 @@ AliGenerator* GeneratorCustom(TString opt = "")
   // 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 = 3;
+  const Int_t nParametrizations = 4;
   TString optNames[nParametrizations] = {
     "single_pp13TeV_1",
+    "single_PbPb5TeV_1",
     "upsilon_pPb8TeV_1",
     "upsilon_Pbp8TeV_1"
   };
   TString optMacros[nParametrizations] = {
     "Muon_GenParamSingle_pp13TeV_1.C",
+    "Muon_GenParamSingle_PbPb5TeV_1.C",
     "Muon_GenParamUpsilon_pPb8TeV_1.C",
     "Muon_GenParamUpsilon_Pbp8TeV_1.C"
   };
   TString optSetups[nParametrizations] = {
+    "",
     "",
     "EvtGen",
     "EvtGen"
diff --git a/MC/CustomGenerators/PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C b/MC/CustomGenerators/PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C
new file mode 100644
index 00000000..43a0b3ac
--- /dev/null
+++ b/MC/CustomGenerators/PWGDQ/Muon_GenParamSingle_PbPb5TeV_1.C
@@ -0,0 +1,88 @@
+//==========================================//
+// Parameterized generator for single muons //
+//==========================================//
+// Info about this particular tunning:      //
+// - PbPb@5.02TeV data (LHC15o)             //
+// - CMSL7 events with physics selection    //
+// - 50-90% centrality class (V0M)          //
+// - Standard muon selection with pDCA      //
+// - pT > 1GeV/c                            //
+//==========================================//
+
+#include <Riostream.h>
+#include "TRandom.h"
+#include "AliGenerator.h"
+#include "AliGenParam.h"
+
+Int_t IpMuon(TRandom *ran);
+Double_t PtMuon(const Double_t *px, const Double_t */*dummy*/);
+Double_t YMuon(const Double_t *py, const Double_t */*dummy*/);
+Double_t V2Muon(const Double_t */*dummy*/, const Double_t */*dummy*/);
+
+//-------------------------------------------------------------------------
+AliGenerator* CreateGenerator()
+{
+  AliGenParam *singleMu = new AliGenParam(2,-1,PtMuon,YMuon,V2Muon,IpMuon);
+  singleMu->SetMomentumRange(0., 1.e6);
+  singleMu->SetPtRange(0.8, 999.);
+  singleMu->SetYRange(-4.3, -2.3);
+  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;
+
+  Double_t p0 = 811.367;
+  Double_t p1 = 0.804372;
+  Double_t p2 = 0.614056;
+  Double_t p3 = 10.4864;
+  Double_t p4 = -0.000650586;
+  Double_t p5 = -1.72877;
+
+  return p0 * (1. / TMath::Power(p1 + TMath::Power(x,p2), p3) + p4 * TMath::Exp(p5*x));
+}
+
+//-------------------------------------------------------------------------
+Double_t YMuon(const Double_t *py, const Double_t */*dummy*/)
+{
+  // muon y
+
+  Double_t y = *py;
+
+  Double_t p0 = 1.81086;
+  Double_t p1 = 0.;
+  Double_t p2 = -0.0983371;
+  Double_t p3 = 0.;
+  Double_t p4 = 0.00270169;
+
+  return p0 * (1. + p1*y + p2*y*y + p3*y*y*y + p4*y*y*y*y);
+}
+
+//-------------------------------------------------------------------------
+Double_t V2Muon(const Double_t */*dummy*/, const Double_t */*dummy*/)
+{
+  // muon v2
+  return 0.;
+}