[u/mrichter/AliRoot.git] / PYTHIA6 / AliPythia.cxx

/**************************************************************************
 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
 *                                                                        *
 * Author: The ALICE Off-line Project.                                    *
 * Contributors are mentioned in the code where appropriate.              *
 *                                                                        *
 * Permission to use, copy, modify and distribute this software and its   *
 * documentation strictly for non-commercial purposes is hereby granted   *
 * without fee, provided that the above copyright notice appears in all   *
 * copies and that both the copyright notice and this permission notice   *
 * appear in the supporting documentation. The authors make no claims     *
 * about the suitability of this software for any purpose. It is          *
 * provided "as is" without express or implied warranty.                  *
 **************************************************************************/

/* $Id$ */

#include "AliPythia.h"
#include "AliPythiaRndm.h"

ClassImp(AliPythia)

#ifndef WIN32
# define pyclus pyclus_
# define pycell pycell_
# define type_of_call
#else
# define pyclus PYCLUS
# define pycell PYCELL
# define type_of_call _stdcall
#endif

extern "C" void type_of_call pyclus(Int_t & );
extern "C" void type_of_call pycell(Int_t & );

//_____________________________________________________________________________

AliPythia* AliPythia::fgAliPythia=NULL;

AliPythia::AliPythia()
{
// Default Constructor
//
//  Set random number
    if (!AliPythiaRndm::GetPythiaRandom()) 
      AliPythiaRndm::SetPythiaRandom(GetRandom());

}

void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
{
// Initialise the process to generate 
    if (!AliPythiaRndm::GetPythiaRandom()) 
      AliPythiaRndm::SetPythiaRandom(GetRandom());
    
    fProcess = process;
    fEcms = energy;
    fStrucFunc = strucfunc;
//  don't decay p0
    SetMDCY(Pycomp(111),1,0);
//  select structure function 
    SetMSTP(52,2);
    SetMSTP(51,strucfunc);
//
// Pythia initialisation for selected processes//
//
// Make MSEL clean
//
    for (Int_t i=1; i<= 200; i++) {
	SetMSUB(i,0);
    }
//  select charm production
    switch (process) 
    {
    case kPyCharm:
	SetMSEL(4);
//
//  heavy quark masses

	SetPMAS(4,1,1.2);
	SetMSTU(16,2);
//
//    primordial pT
	SetMSTP(91,1);
	SetPARP(91,1.);
	SetPARP(93,5.);
//
	break;
    case kPyBeauty:
	SetMSEL(5);
	SetPMAS(5,1,4.75);
	SetMSTU(16,2);
	break;
    case kPyJpsi:
	SetMSEL(0);
// gg->J/Psi g
	SetMSUB(86,1);
	break;
    case kPyJpsiChi:
	SetMSEL(0);
// gg->J/Psi g
	SetMSUB(86,1);
// gg-> chi_0c g
	SetMSUB(87,1);
// gg-> chi_1c g
	SetMSUB(88,1);
// gg-> chi_2c g
	SetMSUB(89,1);	
	break;
    case kPyCharmUnforced:
	SetMSEL(0);
// gq->qg   
	SetMSUB(28,1);
// gg->qq
	SetMSUB(53,1);
// gg->gg
	SetMSUB(68,1);
	break;
    case kPyBeautyUnforced:
	SetMSEL(0);
// gq->qg   
	SetMSUB(28,1);
// gg->qq
	SetMSUB(53,1);
// gg->gg
	SetMSUB(68,1);
	break;
    case kPyMb:
// Minimum Bias pp-Collisions
//
//   
//      select Pythia min. bias model
	SetMSEL(0);
	SetMSUB(92,1);      // single diffraction AB-->XB
	SetMSUB(93,1);      // single diffraction AB-->AX
	SetMSUB(94,1);      // double diffraction
	SetMSUB(95,1);	    // low pt production
	SetMSTP(81,1);      // multiple interactions switched on
	SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
	SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
                            // multiple interaction at a reference energy = 14000 GeV
	SetPARP(89,14000.); // reference energy for the above parameter
	SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
    case kPyMbNonDiffr:
// Minimum Bias pp-Collisions
//
//   
//      select Pythia min. bias model
	SetMSEL(0);
	SetMSUB(95,1);	    // low pt production
	SetMSTP(81,1);      // multiple interactions switched on
	SetMSTP(82,3);      // model with varying impact param. & a single Gaussian
	SetPARP(82,3.47);   // set value pT_0  for turn-off of the cross section of                  
                            // multiple interaction at a reference energy = 14000 GeV
	SetPARP(89,14000.); // reference energy for the above parameter
	SetPARP(90,0.174);  // set exponent for energy dependence of pT_0
 
	break;
    case kPyJets:
//
//  QCD Jets
//
	SetMSEL(1);
	break;
    case kPyDirectGamma:
	SetMSEL(10);
	break;
    case kPyCharmPbPbMNR:
    case kPyD0PbPbMNR:
      // Tuning of Pythia parameters aimed to get a resonable agreement
      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
      // c-cbar single inclusive and double differential distributions.
      // This parameter settings are meant to work with Pb-Pb collisions
      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.

      // All QCD processes
      SetMSEL(1);

      // No multiple interactions
      SetMSTP(81,0);
      SetPARP(81,0.0);
      SetPARP(82,0.0);

      // Initial/final parton shower on (Pythia default)
      SetMSTP(61,1);
      SetMSTP(71,1);

      // 2nd order alpha_s
      SetMSTP(2,2);

      // QCD scales
      SetMSTP(32,2);
      SetPARP(34,1.0);

      // Intrinsic <kT>
      SetMSTP(91,1);
      SetPARP(91,1.304);
      SetPARP(93,6.52);

      // Set c-quark mass
      SetPMAS(4,1,1.2);

      break;
    case kPyCharmpPbMNR:
    case kPyD0pPbMNR:
      // Tuning of Pythia parameters aimed to get a resonable agreement
      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
      // c-cbar single inclusive and double differential distributions.
      // This parameter settings are meant to work with p-Pb collisions
      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.

      // All QCD processes
      SetMSEL(1);

      // No multiple interactions
      SetMSTP(81,0);
      SetPARP(81,0.0);
      SetPARP(82,0.0);

      // Initial/final parton shower on (Pythia default)
      SetMSTP(61,1);
      SetMSTP(71,1);

      // 2nd order alpha_s
      SetMSTP(2,2);

      // QCD scales
      SetMSTP(32,2);
      SetPARP(34,1.0);

      // Intrinsic <kT>
      SetMSTP(91,1);
      SetPARP(91,1.16);
      SetPARP(93,5.8);

      // Set c-quark mass
      SetPMAS(4,1,1.2);

      break;
    case kPyCharmppMNR:
    case kPyD0ppMNR:
      // Tuning of Pythia parameters aimed to get a resonable agreement
      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
      // c-cbar single inclusive and double differential distributions.
      // This parameter settings are meant to work with pp collisions
      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
      // has to be set to 2.1GeV. Example in ConfigCharmPPR.C.

      // All QCD processes
      SetMSEL(1);

      // No multiple interactions
      SetMSTP(81,0);
      SetPARP(81,0.0);
      SetPARP(82,0.0);

      // Initial/final parton shower on (Pythia default)
      SetMSTP(61,1);
      SetMSTP(71,1);

      // 2nd order alpha_s
      SetMSTP(2,2);

      // QCD scales
      SetMSTP(32,2);
      SetPARP(34,1.0);

      // Intrinsic <kT^2>
      SetMSTP(91,1);
      SetPARP(91,1.);
      SetPARP(93,5.);

      // Set c-quark mass
      SetPMAS(4,1,1.2);

      break;
    case kPyBeautyPbPbMNR:
      // Tuning of Pythia parameters aimed to get a resonable agreement
      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
      // b-bbar single inclusive and double differential distributions.
      // This parameter settings are meant to work with Pb-Pb collisions
      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.

      // All QCD processes
      SetMSEL(1);

      // No multiple interactions
      SetMSTP(81,0);
      SetPARP(81,0.0);
      SetPARP(82,0.0);

      // Initial/final parton shower on (Pythia default)
      SetMSTP(61,1);
      SetMSTP(71,1);

      // 2nd order alpha_s
      SetMSTP(2,2);

      // QCD scales
      SetMSTP(32,2);
      SetPARP(34,1.0);
      SetPARP(67,1.0);
      SetPARP(71,1.0);

      // Intrinsic <kT>
      SetMSTP(91,1);
      SetPARP(91,2.035);
      SetPARP(93,10.17);

      // Set b-quark mass
      SetPMAS(5,1,4.75);

      break;
    case kPyBeautypPbMNR:
      // Tuning of Pythia parameters aimed to get a resonable agreement
      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
      // b-bbar single inclusive and double differential distributions.
      // This parameter settings are meant to work with p-Pb collisions
      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.

      // All QCD processes
      SetMSEL(1);

      // No multiple interactions
      SetMSTP(81,0);
      SetPARP(81,0.0);
      SetPARP(82,0.0);

      // Initial/final parton shower on (Pythia default)
      SetMSTP(61,1);
      SetMSTP(71,1);

      // 2nd order alpha_s
      SetMSTP(2,2);

      // QCD scales
      SetMSTP(32,2);
      SetPARP(34,1.0);
      SetPARP(67,1.0);
      SetPARP(71,1.0);

      // Intrinsic <kT>
      SetMSTP(91,1);
      SetPARP(91,1.60);
      SetPARP(93,8.00);

      // Set b-quark mass
      SetPMAS(5,1,4.75);

      break;
    case kPyBeautyppMNR:
      // Tuning of Pythia parameters aimed to get a resonable agreement
      // between with the NLO calculation by Mangano, Nason, Ridolfi for the
      // b-bbar single inclusive and double differential distributions.
      // This parameter settings are meant to work with pp collisions
      // (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
      // To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
      // has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.

      // All QCD processes
      SetMSEL(1);

      // No multiple interactions
      SetMSTP(81,0);
      SetPARP(81,0.0);
      SetPARP(82,0.0);

      // Initial/final parton shower on (Pythia default)
      SetMSTP(61,1);
      SetMSTP(71,1);

      // 2nd order alpha_s
      SetMSTP(2,2);

      // QCD scales
      SetMSTP(32,2);
      SetPARP(34,1.0);
      SetPARP(67,1.0);
      SetPARP(71,1.0);

      // Intrinsic <kT>
      SetMSTP(91,1);
      SetPARP(91,1.);
      SetPARP(93,5.);

      // Set b-quark mass
      SetPMAS(5,1,4.75);

      break;
    }
//
//  Initialize PYTHIA
    SetMSTP(41,1);   // all resonance decays switched on

    Initialize("CMS","p","p",fEcms);

}

Int_t AliPythia::CheckedLuComp(Int_t kf)
{
// Check Lund particle code (for debugging)
    Int_t kc=Pycomp(kf);
    printf("\n Lucomp kf,kc %d %d",kf,kc);
    return kc;
}

void AliPythia::SetNuclei(Int_t a1, Int_t a2)
{
// Treat protons as inside nuclei with mass numbers a1 and a2  
//    The MSTP array in the PYPARS common block is used to enable and 
//    select the nuclear structure functions. 
//    MSTP(52)  : (D=1) choice of proton and nuclear structure-function library
//            =1: internal PYTHIA acording to MSTP(51) 
//            =2: PDFLIB proton  s.f., with MSTP(51)  = 1000xNGROUP+NSET
//    If the following mass number both not equal zero, nuclear corrections of the stf are used.
//    MSTP(192) : Mass number of nucleus side 1
//    MSTP(193) : Mass number of nucleus side 2
    SetMSTP(52,2);
    SetMSTP(192, a1);
    SetMSTP(193, a2);  
}
	

AliPythia* AliPythia::Instance()
{ 
// Set random number generator 
    if (fgAliPythia) {
	return fgAliPythia;
    } else {
	fgAliPythia = new AliPythia();
	return fgAliPythia;
    }
}

void AliPythia::PrintParticles()
{ 
// Print list of particl properties
    Int_t np = 0;
    
    for (Int_t kf=0; kf<1000000; kf++) {
	for (Int_t c = 1;  c > -2; c-=2) {
	    
	    Int_t kc = Pycomp(c*kf);
	    if (kc) {
		Float_t mass  = GetPMAS(kc,1);
		Float_t width = GetPMAS(kc,2);	
		Float_t tau   = GetPMAS(kc,4);
		
		char*   name = new char[8];
		Pyname(kf,name);
	
		np++;
		
		printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e", 
		       c*kf, name, mass, width, tau);
	    }
	}
    }
    printf("\n Number of particles %d \n \n", np);
}

void  AliPythia::ResetDecayTable()
{
//  Set default values for pythia decay switches
    Int_t i;
    for (i = 1; i <  501; i++) SetMDCY(i,1,fDefMDCY[i]);
    for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
}

void  AliPythia::SetDecayTable()
{
//  Set default values for pythia decay switches
//
    Int_t i;
    for (i = 1; i <  501; i++) fDefMDCY[i] = GetMDCY(i,1);
    for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
}

void  AliPythia::Pyclus(Int_t& njet)
{
//  Call Pythia clustering algorithm
//
    pyclus(njet);
}

void  AliPythia::Pycell(Int_t& njet)
{
//  Call Pythia jet reconstruction algorithm
//
    pycell(njet);
}
Commit	Line	Data
8d2cd130	1	/**************************************************************************
	2	* Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
	3	* *
	4	* Author: The ALICE Off-line Project. *
	5	* Contributors are mentioned in the code where appropriate. *
	6	* *
	7	* Permission to use, copy, modify and distribute this software and its *
	8	* documentation strictly for non-commercial purposes is hereby granted *
	9	* without fee, provided that the above copyright notice appears in all *
	10	* copies and that both the copyright notice and this permission notice *
	11	* appear in the supporting documentation. The authors make no claims *
	12	* about the suitability of this software for any purpose. It is *
	13	* provided "as is" without express or implied warranty. *
	14	**************************************************************************/
	15
7cdba479	16	/* $Id$ */
8d2cd130	17
8d2cd130	18	#include "AliPythia.h"
7cdba479	19	#include "AliPythiaRndm.h"
8d2cd130	20
	21	ClassImp(AliPythia)
	22
	23	#ifndef WIN32
	24	# define pyclus pyclus_
	25	# define pycell pycell_
	26	# define type_of_call
	27	#else
	28	# define pyclus PYCLUS
	29	# define pycell PYCELL
	30	# define type_of_call _stdcall
	31	#endif
	32
	33	extern "C" void type_of_call pyclus(Int_t & );
	34	extern "C" void type_of_call pycell(Int_t & );
	35
	36	//_____________________________________________________________________________
	37
	38	AliPythia* AliPythia::fgAliPythia=NULL;
	39
	40	AliPythia::AliPythia()
	41	{
	42	// Default Constructor
	43	//
	44	// Set random number
7cdba479	45	if (!AliPythiaRndm::GetPythiaRandom())
7cdba479	46	AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130	47
	48	}
	49
	50	void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
	51	{
	52	// Initialise the process to generate
7cdba479	53	if (!AliPythiaRndm::GetPythiaRandom())
7cdba479	54	AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130	55
	56	fProcess = process;
	57	fEcms = energy;
	58	fStrucFunc = strucfunc;
	59	// don't decay p0
	60	SetMDCY(Pycomp(111),1,0);
	61	// select structure function
	62	SetMSTP(52,2);
	63	SetMSTP(51,strucfunc);
	64	//
	65	// Pythia initialisation for selected processes//
	66	//
	67	// Make MSEL clean
	68	//
	69	for (Int_t i=1; i<= 200; i++) {
	70	SetMSUB(i,0);
	71	}
	72	// select charm production
	73	switch (process)
	74	{
	75	case kPyCharm:
	76	SetMSEL(4);
	77	//
	78	// heavy quark masses
	79
	80	SetPMAS(4,1,1.2);
	81	SetMSTU(16,2);
	82	//
	83	// primordial pT
	84	SetMSTP(91,1);
	85	SetPARP(91,1.);
	86	SetPARP(93,5.);
	87	//
	88	break;
	89	case kPyBeauty:
	90	SetMSEL(5);
	91	SetPMAS(5,1,4.75);
	92	SetMSTU(16,2);
	93	break;
	94	case kPyJpsi:
	95	SetMSEL(0);
	96	// gg->J/Psi g
	97	SetMSUB(86,1);
	98	break;
	99	case kPyJpsiChi:
	100	SetMSEL(0);
	101	// gg->J/Psi g
	102	SetMSUB(86,1);
	103	// gg-> chi_0c g
	104	SetMSUB(87,1);
	105	// gg-> chi_1c g
	106	SetMSUB(88,1);
	107	// gg-> chi_2c g
	108	SetMSUB(89,1);
	109	break;
	110	case kPyCharmUnforced:
	111	SetMSEL(0);
	112	// gq->qg
	113	SetMSUB(28,1);
	114	// gg->qq
	115	SetMSUB(53,1);
	116	// gg->gg
	117	SetMSUB(68,1);
	118	break;
119	case kPyBeautyUnforced:
120	SetMSEL(0);
121	// gq->qg
122	SetMSUB(28,1);
123	// gg->qq
124	SetMSUB(53,1);
125	// gg->gg
126	SetMSUB(68,1);
127	break;
128	case kPyMb:
129	// Minimum Bias pp-Collisions
130	//
131	//
132	// select Pythia min. bias model
133	SetMSEL(0);
134	SetMSUB(92,1); // single diffraction AB-->XB
135	SetMSUB(93,1); // single diffraction AB-->AX
136	SetMSUB(94,1); // double diffraction
137	SetMSUB(95,1); // low pt production
138	SetMSTP(81,1); // multiple interactions switched on
139	SetMSTP(82,3); // model with varying impact param. & a single Gaussian
140	SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
141	// multiple interaction at a reference energy = 14000 GeV
142	SetPARP(89,14000.); // reference energy for the above parameter
143	SetPARP(90,0.174); // set exponent for energy dependence of pT_0
144	case kPyMbNonDiffr:
145	// Minimum Bias pp-Collisions
146	//
147	//
148	// select Pythia min. bias model
149	SetMSEL(0);
150	SetMSUB(95,1); // low pt production
151	SetMSTP(81,1); // multiple interactions switched on
152	SetMSTP(82,3); // model with varying impact param. & a single Gaussian
153	SetPARP(82,3.47); // set value pT_0 for turn-off of the cross section of
154	// multiple interaction at a reference energy = 14000 GeV
155	SetPARP(89,14000.); // reference energy for the above parameter
156	SetPARP(90,0.174); // set exponent for energy dependence of pT_0
157
158	break;
159	case kPyJets:
160	//
161	// QCD Jets
162	//
163	SetMSEL(1);
164	break;
165	case kPyDirectGamma:
166	SetMSEL(10);
167	break;
adf4d898	168	case kPyCharmPbPbMNR:
adf4d898	169	case kPyD0PbPbMNR:
8d2cd130	170	// Tuning of Pythia parameters aimed to get a resonable agreement
	171	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	172	// c-cbar single inclusive and double differential distributions.
	173	// This parameter settings are meant to work with Pb-Pb collisions
adf4d898	174	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
8d2cd130	175	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	176	// has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
	177
	178	// All QCD processes
	179	SetMSEL(1);
	180
	181	// No multiple interactions
	182	SetMSTP(81,0);
	183	SetPARP(81,0.0);
	184	SetPARP(82,0.0);
	185
	186	// Initial/final parton shower on (Pythia default)
	187	SetMSTP(61,1);
	188	SetMSTP(71,1);
	189
	190	// 2nd order alpha_s
	191	SetMSTP(2,2);
	192
	193	// QCD scales
	194	SetMSTP(32,2);
	195	SetPARP(34,1.0);
	196
adf4d898	197	// Intrinsic <kT>
8d2cd130	198	SetMSTP(91,1);
	199	SetPARP(91,1.304);
	200	SetPARP(93,6.52);
	201
	202	// Set c-quark mass
	203	SetPMAS(4,1,1.2);
	204
	205	break;
adf4d898	206	case kPyCharmpPbMNR:
	207	case kPyD0pPbMNR:
	208	// Tuning of Pythia parameters aimed to get a resonable agreement
	209	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	210	// c-cbar single inclusive and double differential distributions.
	211	// This parameter settings are meant to work with p-Pb collisions
	212	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	213	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	214	// has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
	215
	216	// All QCD processes
	217	SetMSEL(1);
	218
	219	// No multiple interactions
	220	SetMSTP(81,0);
	221	SetPARP(81,0.0);
	222	SetPARP(82,0.0);
	223
	224	// Initial/final parton shower on (Pythia default)
	225	SetMSTP(61,1);
	226	SetMSTP(71,1);
	227
	228	// 2nd order alpha_s
	229	SetMSTP(2,2);
	230
	231	// QCD scales
	232	SetMSTP(32,2);
	233	SetPARP(34,1.0);
	234
	235	// Intrinsic <kT>
	236	SetMSTP(91,1);
	237	SetPARP(91,1.16);
	238	SetPARP(93,5.8);
	239
	240	// Set c-quark mass
	241	SetPMAS(4,1,1.2);
	242
	243	break;
	244	case kPyCharmppMNR:
	245	case kPyD0ppMNR:
	246	// Tuning of Pythia parameters aimed to get a resonable agreement
	247	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	248	// c-cbar single inclusive and double differential distributions.
	249	// This parameter settings are meant to work with pp collisions
	250	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	251	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	252	// has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
	253
	254	// All QCD processes
	255	SetMSEL(1);
	256
	257	// No multiple interactions
	258	SetMSTP(81,0);
	259	SetPARP(81,0.0);
	260	SetPARP(82,0.0);
	261
	262	// Initial/final parton shower on (Pythia default)
	263	SetMSTP(61,1);
	264	SetMSTP(71,1);
	265
	266	// 2nd order alpha_s
	267	SetMSTP(2,2);
	268
	269	// QCD scales
270	SetMSTP(32,2);
271	SetPARP(34,1.0);
272
273	// Intrinsic <kT^2>
274	SetMSTP(91,1);
275	SetPARP(91,1.);
276	SetPARP(93,5.);
277
278	// Set c-quark mass
279	SetPMAS(4,1,1.2);
280
281	break;
282	case kPyBeautyPbPbMNR:
8d2cd130	283	// Tuning of Pythia parameters aimed to get a resonable agreement
	284	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	285	// b-bbar single inclusive and double differential distributions.
	286	// This parameter settings are meant to work with Pb-Pb collisions
	287	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	288	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	289	// has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
	290
	291	// All QCD processes
	292	SetMSEL(1);
	293
	294	// No multiple interactions
	295	SetMSTP(81,0);
	296	SetPARP(81,0.0);
	297	SetPARP(82,0.0);
	298
	299	// Initial/final parton shower on (Pythia default)
	300	SetMSTP(61,1);
	301	SetMSTP(71,1);
	302
	303	// 2nd order alpha_s
	304	SetMSTP(2,2);
	305
	306	// QCD scales
	307	SetMSTP(32,2);
	308	SetPARP(34,1.0);
	309	SetPARP(67,1.0);
	310	SetPARP(71,1.0);
	311
adf4d898	312	// Intrinsic <kT>
8d2cd130	313	SetMSTP(91,1);
	314	SetPARP(91,2.035);
	315	SetPARP(93,10.17);
	316
	317	// Set b-quark mass
	318	SetPMAS(5,1,4.75);
	319
adf4d898	320	break;
	321	case kPyBeautypPbMNR:
	322	// Tuning of Pythia parameters aimed to get a resonable agreement
	323	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	324	// b-bbar single inclusive and double differential distributions.
	325	// This parameter settings are meant to work with p-Pb collisions
	326	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	327	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	328	// has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
	329
	330	// All QCD processes
	331	SetMSEL(1);
	332
	333	// No multiple interactions
	334	SetMSTP(81,0);
	335	SetPARP(81,0.0);
	336	SetPARP(82,0.0);
	337
	338	// Initial/final parton shower on (Pythia default)
	339	SetMSTP(61,1);
	340	SetMSTP(71,1);
	341
	342	// 2nd order alpha_s
	343	SetMSTP(2,2);
	344
	345	// QCD scales
	346	SetMSTP(32,2);
	347	SetPARP(34,1.0);
	348	SetPARP(67,1.0);
	349	SetPARP(71,1.0);
	350
	351	// Intrinsic <kT>
	352	SetMSTP(91,1);
	353	SetPARP(91,1.60);
	354	SetPARP(93,8.00);
	355
	356	// Set b-quark mass
	357	SetPMAS(5,1,4.75);
	358
	359	break;
	360	case kPyBeautyppMNR:
	361	// Tuning of Pythia parameters aimed to get a resonable agreement
	362	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	363	// b-bbar single inclusive and double differential distributions.
	364	// This parameter settings are meant to work with pp collisions
	365	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	366	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	367	// has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
	368
	369	// All QCD processes
	370	SetMSEL(1);
	371
	372	// No multiple interactions
	373	SetMSTP(81,0);
	374	SetPARP(81,0.0);
	375	SetPARP(82,0.0);
	376
	377	// Initial/final parton shower on (Pythia default)
	378	SetMSTP(61,1);
	379	SetMSTP(71,1);
	380
	381	// 2nd order alpha_s
	382	SetMSTP(2,2);
	383
384	// QCD scales
385	SetMSTP(32,2);
386	SetPARP(34,1.0);
387	SetPARP(67,1.0);
388	SetPARP(71,1.0);
389
390	// Intrinsic <kT>
391	SetMSTP(91,1);
392	SetPARP(91,1.);
393	SetPARP(93,5.);
394
395	// Set b-quark mass
396	SetPMAS(5,1,4.75);
397
8d2cd130	398	break;
	399	}
	400	//
	401	// Initialize PYTHIA
	402	SetMSTP(41,1); // all resonance decays switched on
	403
	404	Initialize("CMS","p","p",fEcms);
	405
	406	}
	407
	408	Int_t AliPythia::CheckedLuComp(Int_t kf)
	409	{
	410	// Check Lund particle code (for debugging)
	411	Int_t kc=Pycomp(kf);
	412	printf("\n Lucomp kf,kc %d %d",kf,kc);
	413	return kc;
	414	}
	415
	416	void AliPythia::SetNuclei(Int_t a1, Int_t a2)
	417	{
	418	// Treat protons as inside nuclei with mass numbers a1 and a2
	419	// The MSTP array in the PYPARS common block is used to enable and
	420	// select the nuclear structure functions.
	421	// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
	422	// =1: internal PYTHIA acording to MSTP(51)
	423	// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
	424	// If the following mass number both not equal zero, nuclear corrections of the stf are used.
	425	// MSTP(192) : Mass number of nucleus side 1
	426	// MSTP(193) : Mass number of nucleus side 2
	427	SetMSTP(52,2);
	428	SetMSTP(192, a1);
	429	SetMSTP(193, a2);
	430	}
	431
	432
	433	AliPythia* AliPythia::Instance()
	434	{
	435	// Set random number generator
	436	if (fgAliPythia) {
	437	return fgAliPythia;
	438	} else {
	439	fgAliPythia = new AliPythia();
	440	return fgAliPythia;
	441	}
	442	}
	443
	444	void AliPythia::PrintParticles()
	445	{
	446	// Print list of particl properties
	447	Int_t np = 0;
	448
	449	for (Int_t kf=0; kf<1000000; kf++) {
	450	for (Int_t c = 1; c > -2; c-=2) {
	451
	452	Int_t kc = Pycomp(c*kf);
	453	if (kc) {
	454	Float_t mass = GetPMAS(kc,1);
	455	Float_t width = GetPMAS(kc,2);
	456	Float_t tau = GetPMAS(kc,4);
	457
	458	char* name = new char[8];
	459	Pyname(kf,name);
	460
	461	np++;
462
463	printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
464	c*kf, name, mass, width, tau);
465	}
466	}
467	}
468	printf("\n Number of particles %d \n \n", np);
469	}
470
471	void AliPythia::ResetDecayTable()
472	{
473	// Set default values for pythia decay switches
474	Int_t i;
475	for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
476	for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
477	}
478
479	void AliPythia::SetDecayTable()
480	{
481	// Set default values for pythia decay switches
482	//
483	Int_t i;
484	for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
485	for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
486	}
487
488	void AliPythia::Pyclus(Int_t& njet)
489	{
490	// Call Pythia clustering algorithm
491	//
492	pyclus(njet);
493	}
494
495	void AliPythia::Pycell(Int_t& njet)
496	{
497	// Call Pythia jet reconstruction algorithm
498	//
499	pycell(njet);
500	}
501
502
503