[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 pyshow pyshow_
# define pyrobo pyrobo_
# define type_of_call
#else
# define pyclus PYCLUS
# define pycell PYCELL
# define pyrobo PYROBO
# define type_of_call _stdcall
#endif

extern "C" void type_of_call pyclus(Int_t & );
extern "C" void type_of_call pycell(Int_t & );
extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_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

//
// ATLAS Tuning
//
        SetMSTP(51, 7);            // CTEQ5L pdf
	SetMSTP(81,1);             // Multiple Interactions ON
	SetMSTP(82,4);             // Double Gaussian Model

	SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
	SetPARP(89,1000.);         // [GeV]   Ref. energy
	SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
	SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
	SetPARP(84,0.5);           // Core radius
	SetPARP(85,0.33);          // Regulates gluon prod. mechanism
	SetPARP(86,0.66);          // Regulates gluon prod. mechanism
	SetPARP(67,1);             // Regulates Initial State Radiation
	break;
    case kPyMbNonDiffr:
// Minimum Bias pp-Collisions
//
//   
//      select Pythia min. bias model
	SetMSEL(0);
	SetMSUB(95,1);	           // low pt production
	
	SetMSTP(51, 7);            // CTEQ5L pdf
	SetMSTP(81,1);             // Multiple Interactions ON
	SetMSTP(82,4);             // Double Gaussian Model

	SetPARP(82,1.8);           // [GeV]    PT_min at Ref. energy
	SetPARP(89,1000.);         // [GeV]   Ref. energy
	SetPARP(90,0.16);          // 2*epsilon (exponent in power law)
	SetPARP(83,0.5);           // Core density in proton matter distribution (def.value)
	SetPARP(84,0.5);           // Core radius
	SetPARP(85,0.33);          // Regulates gluon prod. mechanism
	SetPARP(86,0.66);          // Regulates gluon prod. mechanism
	SetPARP(67,1);             // Regulates Initial State Radiation
	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;
    char*   name = new char[16];    
    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);

		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);
}

void  AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
{
//  Call Pythia jet reconstruction algorithm
//
    Int_t numpart   = fPyjets->N;
    for (Int_t i = 0; i < numpart; i++)
    {
	if (fPyjets->K[2][i] == 7) ip1 = i+1;
	if (fPyjets->K[2][i] == 8) ip2 = i+1;
    }
    

    qmax = 2. * GetVINT(51);
    printf("Pyshow %d %d %f", ip1, ip2, qmax);
    
    pyshow(ip1, ip2, qmax);
}

void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
{
    pyrobo(imi, ima, the, phi, bex, bey, bez);
}


void  AliPythia::Quench()
{
//
//
//  Simple Jet Quenching routine:
//  =============================
//  The jet formed by all final state partons radiated by the parton created 
//  in the hard collisions is quenched by a factor z using:
//  (E + p_z)new = (1-z) (E + p_z)old
//
//   The lost momentum is first balanced by one gluon with virtuality > 0.   
//   Subsequently the gluon splits to yield two gluons with E = p.
//
    Float_t p0[2][5];
    Float_t p1[2][5];
    Float_t p2[2][5];
    Int_t   klast[2] = {-1, -1};
    Int_t   kglu[2];
    for (Int_t i = 0; i < 4; i++)
    {
	p0[0][i] = 0.;
	p0[1][i] = 0.;
	p1[0][i] = 0.;
	p1[1][i] = 0.;
	p2[0][i] = 0.;
	p2[1][i] = 0.;
    }

    Int_t numpart   = fPyjets->N;

    for (Int_t i = 0; i < numpart; i++)
    {
	Int_t imo =  fPyjets->K[2][i];
	Int_t kst =  fPyjets->K[0][i];
	Int_t pdg =  fPyjets->K[1][i];
	
//      Quarks and gluons only
	if (pdg != 21 && TMath::Abs(pdg) > 6) continue;

//      Particles from hard scattering only


	Float_t px = fPyjets->P[0][i];
	Float_t py = fPyjets->P[1][i];
	Float_t pz = fPyjets->P[2][i];
	Float_t e  = fPyjets->P[3][i];
	Float_t m  = fPyjets->P[4][i];
	Float_t pt = TMath::Sqrt(px * px + py * py);
//      Skip comment lines
	if (kst != 1 && kst != 2) continue;

	Float_t mt = TMath::Sqrt(px * px + py * py + m * m);

	//
	// Some cuts to be in a save kinematic region
	//
	if (imo != 7 && imo != 8) continue;
	Int_t index = imo - 7;
	klast[index] = i;
	
	p0[index][0] += px;
	p0[index][1] += py;
	p0[index][2] += pz;
	p0[index][3] += e;
//
// Fix z
//
	
	Float_t z       = 0.2;
	Float_t eppzOld = e + pz;
	Float_t empzOld = e - pz;


	//
	// Kinematics of the original parton
	// 

	Float_t eppzNew = (1. - z) * eppzOld;
	Float_t empzNew = empzOld - mt * mt * z / eppzOld;
	Float_t eNew0    = 0.5 * (eppzNew + empzNew);
	Float_t pzNew0   = 0.5 * (eppzNew - empzNew);
	//
        // Skip if pt too small
	// 
	if (m * m > eppzNew * empzNew) continue;
	Float_t ptNew    = TMath::Sqrt(eppzNew * empzNew - m * m);
	Float_t pxNew0   = ptNew / pt * px;
	Float_t pyNew0   = ptNew / pt * py;	

	p1[index][0] += pxNew0;
	p1[index][1] += pyNew0;
	p1[index][2] += pzNew0;
	p1[index][3] += eNew0;	
	//
	// Update event record
	//
	fPyjets->P[0][i] = pxNew0;
	fPyjets->P[1][i] = pyNew0;
	fPyjets->P[2][i] = pzNew0;
	fPyjets->P[3][i] = eNew0;
	
    }

    //
    // Gluons
    // 
    
    for (Int_t k = 0; k < 2; k++) 
    {
	for (Int_t j = 0; j < 4; j++) 
	{
	    p2[k][j] = p0[k][j] - p1[k][j];
	}
	p2[k][4] = p2[k][3] * p2[k][3] - p2[k][0] * p2[k][0] - p2[k][1] * p2[k][1] - p2[k][2] * p2[k][2];

	if (p2[k][4] > 0.)
	{

	    //
	    // Bring gluon back to mass shell via momentum scaling 
	    // (momentum will not be conserved, but energy)
	    //
	    // not used anymore
/*
	    Float_t psq   = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
	    Float_t fact  = TMath::Sqrt(1. + p2[k][4] / psq);
	    p2[k][0] *= fact;
	    p2[k][1] *= fact;
	    p2[k][2] *= fact;
	    p2[k][3]  = TMath::Sqrt(psq) * fact;
	    p2[k][4]  = 0.;
*/
	}
    }

    if (p2[0][4] > 0.) {
	p2[0][4] = TMath::Sqrt(p2[0][4]);
    } else {
	printf("Warning negative mass squared ! \n");
    }
    
    if (p2[1][4] > 0.) {
	p2[1][4] = TMath::Sqrt(p2[1][4]);
    } else {
	printf("Warning negative mass squared ! \n");
    }
    
    //
    // Add the gluons
    //
    
    
    for (Int_t i = 0; i < 2; i++) {
	Int_t ish, jmin, jmax, iGlu, iNew;
	Int_t in = klast[i];
	ish = 0;
	
	if (in == -1) continue;
	if (i == 1 && klast[1] > klast[0]) in += ish;
	
	jmin = in - 1;
	ish  = 1;

	if (p2[i][4] > 0) ish = 2;
	
	iGlu = in;
	iNew = in + ish;
	jmax = numpart + ish - 1;
 
	if (fPyjets->K[0][in-1] == 1 || fPyjets->K[0][in-1] == 21 || fPyjets->K[0][in-1] == 11) {
	    jmin = in;
	    iGlu = in + 1;
	    iNew = in;
	}

	kglu[i] = iGlu;
	
	for (Int_t j = jmax; j > jmin; j--)
	{
	    for (Int_t k = 0; k < 5; k++) {
		fPyjets->K[k][j] =  fPyjets->K[k][j-ish];
		fPyjets->P[k][j] =  fPyjets->P[k][j-ish];
		fPyjets->V[k][j] =  fPyjets->V[k][j-ish];
	    }
	} // end shifting
	numpart += ish;
	(fPyjets->N) += ish;
	
	if (ish == 1) {
	    fPyjets->P[0][iGlu] = p2[i][0];
	    fPyjets->P[1][iGlu] = p2[i][1];
	    fPyjets->P[2][iGlu] = p2[i][2];
	    fPyjets->P[3][iGlu] = p2[i][3];
	    fPyjets->P[4][iGlu] = p2[i][4];
	    
	    fPyjets->K[0][iGlu] = 2;
	    fPyjets->K[1][iGlu] = 21;	
	    fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
	    fPyjets->K[3][iGlu] = -1;	
	    fPyjets->K[4][iGlu] = -1;
	} else {
	    //
	    // Split gluon in rest frame.
	    //
	    Double_t bx =  p2[i][0] / p2[i][3];
	    Double_t by =  p2[i][1] / p2[i][3];
	    Double_t bz =  p2[i][2] / p2[i][3];
	    
	    Float_t pst  = p2[i][4] / 2.;
	    
	    Float_t cost = 2. * gRandom->Rndm() - 1.;
	    Float_t sint = TMath::Sqrt(1. - cost * cost);
	    Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
	    
	    Float_t pz1 =   pst * cost;
	    Float_t pz2 =  -pst * cost;
	    Float_t pt1 =   pst * sint;
	    Float_t pt2 =  -pst * sint;
	    Float_t px1 =   pt1 * TMath::Cos(phi);
	    Float_t py1 =   pt1 * TMath::Sin(phi);	    
	    Float_t px2 =   pt2 * TMath::Cos(phi);
	    Float_t py2 =   pt2 * TMath::Sin(phi);	    
	    
	    fPyjets->P[0][iGlu] = px1;
	    fPyjets->P[1][iGlu] = py1;
	    fPyjets->P[2][iGlu] = pz1;
	    fPyjets->P[3][iGlu] = pst;
	    fPyjets->P[4][iGlu] = 0.;
	    
	    fPyjets->K[0][iGlu] = 2;
	    fPyjets->K[1][iGlu] = 21;	
	    fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
	    fPyjets->K[3][iGlu] = -1;	
	    fPyjets->K[4][iGlu] = -1;

	    fPyjets->P[0][iGlu+1] = px2;
	    fPyjets->P[1][iGlu+1] = py2;
	    fPyjets->P[2][iGlu+1] = pz2;
	    fPyjets->P[3][iGlu+1] = pst;
	    fPyjets->P[4][iGlu+1] = 0.;
	    
	    fPyjets->K[0][iGlu+1] = 2;
	    fPyjets->K[1][iGlu+1] = 21;	
	    fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew];
	    fPyjets->K[3][iGlu+1] = -1;	
	    fPyjets->K[4][iGlu+1] = -1;
	    SetMSTU(1,0);
	    SetMSTU(2,0);

	    //
	    // Boost back
	    //
	    Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);

	}
    } // end adding gluons
} // end quench
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_
452af8c7	26	# define pyshow pyshow_
452af8c7	27	# define pyrobo pyrobo_
8d2cd130	28	# define type_of_call
	29	#else
	30	# define pyclus PYCLUS
	31	# define pycell PYCELL
452af8c7	32	# define pyrobo PYROBO
8d2cd130	33	# define type_of_call _stdcall
	34	#endif
	35
	36	extern "C" void type_of_call pyclus(Int_t & );
	37	extern "C" void type_of_call pycell(Int_t & );
452af8c7	38	extern "C" void type_of_call pyshow(Int_t &, Int_t &, Double_t &);
452af8c7	39	extern "C" void type_of_call pyrobo(Int_t &, Int_t &, Double_t &, Double_t &, Double_t &, Double_t &, Double_t &);
8d2cd130	40
	41	//_____________________________________________________________________________
	42
	43	AliPythia* AliPythia::fgAliPythia=NULL;
	44
	45	AliPythia::AliPythia()
	46	{
	47	// Default Constructor
	48	//
	49	// Set random number
7cdba479	50	if (!AliPythiaRndm::GetPythiaRandom())
7cdba479	51	AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130	52
	53	}
	54
	55	void AliPythia::ProcInit(Process_t process, Float_t energy, StrucFunc_t strucfunc)
	56	{
	57	// Initialise the process to generate
7cdba479	58	if (!AliPythiaRndm::GetPythiaRandom())
7cdba479	59	AliPythiaRndm::SetPythiaRandom(GetRandom());
8d2cd130	60
	61	fProcess = process;
	62	fEcms = energy;
	63	fStrucFunc = strucfunc;
	64	// don't decay p0
	65	SetMDCY(Pycomp(111),1,0);
	66	// select structure function
	67	SetMSTP(52,2);
	68	SetMSTP(51,strucfunc);
	69	//
	70	// Pythia initialisation for selected processes//
	71	//
	72	// Make MSEL clean
	73	//
	74	for (Int_t i=1; i<= 200; i++) {
	75	SetMSUB(i,0);
	76	}
	77	// select charm production
	78	switch (process)
	79	{
	80	case kPyCharm:
	81	SetMSEL(4);
	82	//
	83	// heavy quark masses
	84
	85	SetPMAS(4,1,1.2);
	86	SetMSTU(16,2);
	87	//
	88	// primordial pT
	89	SetMSTP(91,1);
	90	SetPARP(91,1.);
	91	SetPARP(93,5.);
	92	//
	93	break;
	94	case kPyBeauty:
	95	SetMSEL(5);
	96	SetPMAS(5,1,4.75);
	97	SetMSTU(16,2);
	98	break;
	99	case kPyJpsi:
	100	SetMSEL(0);
	101	// gg->J/Psi g
	102	SetMSUB(86,1);
	103	break;
	104	case kPyJpsiChi:
	105	SetMSEL(0);
	106	// gg->J/Psi g
	107	SetMSUB(86,1);
	108	// gg-> chi_0c g
	109	SetMSUB(87,1);
	110	// gg-> chi_1c g
	111	SetMSUB(88,1);
	112	// gg-> chi_2c g
	113	SetMSUB(89,1);
	114	break;
	115	case kPyCharmUnforced:
	116	SetMSEL(0);
	117	// gq->qg
	118	SetMSUB(28,1);
	119	// gg->qq
	120	SetMSUB(53,1);
	121	// gg->gg
	122	SetMSUB(68,1);
	123	break;
124	case kPyBeautyUnforced:
125	SetMSEL(0);
126	// gq->qg
127	SetMSUB(28,1);
128	// gg->qq
129	SetMSUB(53,1);
130	// gg->gg
131	SetMSUB(68,1);
132	break;
133	case kPyMb:
134	// Minimum Bias pp-Collisions
135	//
136	//
137	// select Pythia min. bias model
138	SetMSEL(0);
511db649	139	SetMSUB(92,1); // single diffraction AB-->XB
	140	SetMSUB(93,1); // single diffraction AB-->AX
	141	SetMSUB(94,1); // double diffraction
	142	SetMSUB(95,1); // low pt production
	143
	144	//
	145	// ATLAS Tuning
	146	//
	147	SetMSTP(51, 7); // CTEQ5L pdf
	148	SetMSTP(81,1); // Multiple Interactions ON
	149	SetMSTP(82,4); // Double Gaussian Model
	150
	151	SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
	152	SetPARP(89,1000.); // [GeV] Ref. energy
	153	SetPARP(90,0.16); // 2*epsilon (exponent in power law)
	154	SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
	155	SetPARP(84,0.5); // Core radius
	156	SetPARP(85,0.33); // Regulates gluon prod. mechanism
	157	SetPARP(86,0.66); // Regulates gluon prod. mechanism
	158	SetPARP(67,1); // Regulates Initial State Radiation
	159	break;
8d2cd130	160	case kPyMbNonDiffr:
	161	// Minimum Bias pp-Collisions
	162	//
	163	//
	164	// select Pythia min. bias model
	165	SetMSEL(0);
511db649	166	SetMSUB(95,1); // low pt production
	167
	168	SetMSTP(51, 7); // CTEQ5L pdf
	169	SetMSTP(81,1); // Multiple Interactions ON
	170	SetMSTP(82,4); // Double Gaussian Model
	171
	172	SetPARP(82,1.8); // [GeV] PT_min at Ref. energy
	173	SetPARP(89,1000.); // [GeV] Ref. energy
	174	SetPARP(90,0.16); // 2*epsilon (exponent in power law)
	175	SetPARP(83,0.5); // Core density in proton matter distribution (def.value)
	176	SetPARP(84,0.5); // Core radius
	177	SetPARP(85,0.33); // Regulates gluon prod. mechanism
	178	SetPARP(86,0.66); // Regulates gluon prod. mechanism
	179	SetPARP(67,1); // Regulates Initial State Radiation
8d2cd130	180	break;
	181	case kPyJets:
	182	//
	183	// QCD Jets
	184	//
	185	SetMSEL(1);
	186	break;
	187	case kPyDirectGamma:
	188	SetMSEL(10);
	189	break;
adf4d898	190	case kPyCharmPbPbMNR:
adf4d898	191	case kPyD0PbPbMNR:
8d2cd130	192	// Tuning of Pythia parameters aimed to get a resonable agreement
	193	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	194	// c-cbar single inclusive and double differential distributions.
	195	// This parameter settings are meant to work with Pb-Pb collisions
adf4d898	196	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
8d2cd130	197	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	198	// has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
	199
	200	// All QCD processes
	201	SetMSEL(1);
	202
	203	// No multiple interactions
	204	SetMSTP(81,0);
	205	SetPARP(81,0.0);
	206	SetPARP(82,0.0);
	207
	208	// Initial/final parton shower on (Pythia default)
	209	SetMSTP(61,1);
	210	SetMSTP(71,1);
	211
	212	// 2nd order alpha_s
	213	SetMSTP(2,2);
	214
	215	// QCD scales
	216	SetMSTP(32,2);
	217	SetPARP(34,1.0);
	218
adf4d898	219	// Intrinsic <kT>
8d2cd130	220	SetMSTP(91,1);
	221	SetPARP(91,1.304);
	222	SetPARP(93,6.52);
	223
	224	// Set c-quark mass
	225	SetPMAS(4,1,1.2);
	226
	227	break;
adf4d898	228	case kPyCharmpPbMNR:
	229	case kPyD0pPbMNR:
	230	// Tuning of Pythia parameters aimed to get a resonable agreement
	231	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	232	// c-cbar single inclusive and double differential distributions.
	233	// This parameter settings are meant to work with p-Pb collisions
	234	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	235	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	236	// has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
	237
	238	// All QCD processes
	239	SetMSEL(1);
	240
	241	// No multiple interactions
	242	SetMSTP(81,0);
	243	SetPARP(81,0.0);
	244	SetPARP(82,0.0);
	245
	246	// Initial/final parton shower on (Pythia default)
	247	SetMSTP(61,1);
	248	SetMSTP(71,1);
	249
	250	// 2nd order alpha_s
	251	SetMSTP(2,2);
	252
	253	// QCD scales
	254	SetMSTP(32,2);
	255	SetPARP(34,1.0);
	256
	257	// Intrinsic <kT>
	258	SetMSTP(91,1);
	259	SetPARP(91,1.16);
	260	SetPARP(93,5.8);
	261
	262	// Set c-quark mass
	263	SetPMAS(4,1,1.2);
	264
	265	break;
	266	case kPyCharmppMNR:
	267	case kPyD0ppMNR:
	268	// Tuning of Pythia parameters aimed to get a resonable agreement
	269	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	270	// c-cbar single inclusive and double differential distributions.
	271	// This parameter settings are meant to work with pp collisions
	272	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	273	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	274	// has to be set to 2.1GeV. Example in ConfigCharmPPR.C.
	275
	276	// All QCD processes
	277	SetMSEL(1);
	278
	279	// No multiple interactions
	280	SetMSTP(81,0);
	281	SetPARP(81,0.0);
	282	SetPARP(82,0.0);
	283
	284	// Initial/final parton shower on (Pythia default)
	285	SetMSTP(61,1);
	286	SetMSTP(71,1);
	287
	288	// 2nd order alpha_s
	289	SetMSTP(2,2);
	290
	291	// QCD scales
292	SetMSTP(32,2);
293	SetPARP(34,1.0);
294
295	// Intrinsic <kT^2>
296	SetMSTP(91,1);
297	SetPARP(91,1.);
298	SetPARP(93,5.);
299
300	// Set c-quark mass
301	SetPMAS(4,1,1.2);
302
303	break;
304	case kPyBeautyPbPbMNR:
8d2cd130	305	// Tuning of Pythia parameters aimed to get a resonable agreement
	306	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	307	// b-bbar single inclusive and double differential distributions.
	308	// This parameter settings are meant to work with Pb-Pb collisions
	309	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	310	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	311	// has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
	312
	313	// All QCD processes
	314	SetMSEL(1);
	315
	316	// No multiple interactions
	317	SetMSTP(81,0);
	318	SetPARP(81,0.0);
	319	SetPARP(82,0.0);
	320
	321	// Initial/final parton shower on (Pythia default)
	322	SetMSTP(61,1);
	323	SetMSTP(71,1);
	324
	325	// 2nd order alpha_s
	326	SetMSTP(2,2);
	327
	328	// QCD scales
	329	SetMSTP(32,2);
	330	SetPARP(34,1.0);
	331	SetPARP(67,1.0);
	332	SetPARP(71,1.0);
	333
adf4d898	334	// Intrinsic <kT>
8d2cd130	335	SetMSTP(91,1);
	336	SetPARP(91,2.035);
	337	SetPARP(93,10.17);
	338
	339	// Set b-quark mass
	340	SetPMAS(5,1,4.75);
	341
adf4d898	342	break;
	343	case kPyBeautypPbMNR:
	344	// Tuning of Pythia parameters aimed to get a resonable agreement
	345	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	346	// b-bbar single inclusive and double differential distributions.
	347	// This parameter settings are meant to work with p-Pb collisions
	348	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	349	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	350	// has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
	351
	352	// All QCD processes
	353	SetMSEL(1);
	354
	355	// No multiple interactions
	356	SetMSTP(81,0);
	357	SetPARP(81,0.0);
	358	SetPARP(82,0.0);
	359
	360	// Initial/final parton shower on (Pythia default)
	361	SetMSTP(61,1);
	362	SetMSTP(71,1);
	363
	364	// 2nd order alpha_s
	365	SetMSTP(2,2);
	366
	367	// QCD scales
	368	SetMSTP(32,2);
	369	SetPARP(34,1.0);
	370	SetPARP(67,1.0);
	371	SetPARP(71,1.0);
	372
	373	// Intrinsic <kT>
	374	SetMSTP(91,1);
	375	SetPARP(91,1.60);
	376	SetPARP(93,8.00);
	377
	378	// Set b-quark mass
	379	SetPMAS(5,1,4.75);
	380
	381	break;
	382	case kPyBeautyppMNR:
	383	// Tuning of Pythia parameters aimed to get a resonable agreement
	384	// between with the NLO calculation by Mangano, Nason, Ridolfi for the
	385	// b-bbar single inclusive and double differential distributions.
	386	// This parameter settings are meant to work with pp collisions
	387	// (AliGenPythia::SetNuclei) and with kCTEQ4L PDFs.
	388	// To get a good agreement the minimum ptHard (AliGenPythia::SetPtHard)
	389	// has to be set to 2.75GeV. Example in ConfigBeautyPPR.C.
	390
	391	// All QCD processes
	392	SetMSEL(1);
	393
	394	// No multiple interactions
	395	SetMSTP(81,0);
	396	SetPARP(81,0.0);
	397	SetPARP(82,0.0);
	398
	399	// Initial/final parton shower on (Pythia default)
	400	SetMSTP(61,1);
	401	SetMSTP(71,1);
	402
	403	// 2nd order alpha_s
	404	SetMSTP(2,2);
	405
406	// QCD scales
407	SetMSTP(32,2);
408	SetPARP(34,1.0);
409	SetPARP(67,1.0);
410	SetPARP(71,1.0);
411
412	// Intrinsic <kT>
413	SetMSTP(91,1);
414	SetPARP(91,1.);
415	SetPARP(93,5.);
416
417	// Set b-quark mass
418	SetPMAS(5,1,4.75);
419
8d2cd130	420	break;
	421	}
	422	//
	423	// Initialize PYTHIA
	424	SetMSTP(41,1); // all resonance decays switched on
	425
	426	Initialize("CMS","p","p",fEcms);
	427
	428	}
	429
	430	Int_t AliPythia::CheckedLuComp(Int_t kf)
	431	{
	432	// Check Lund particle code (for debugging)
	433	Int_t kc=Pycomp(kf);
	434	printf("\n Lucomp kf,kc %d %d",kf,kc);
	435	return kc;
	436	}
	437
	438	void AliPythia::SetNuclei(Int_t a1, Int_t a2)
	439	{
	440	// Treat protons as inside nuclei with mass numbers a1 and a2
	441	// The MSTP array in the PYPARS common block is used to enable and
	442	// select the nuclear structure functions.
	443	// MSTP(52) : (D=1) choice of proton and nuclear structure-function library
	444	// =1: internal PYTHIA acording to MSTP(51)
	445	// =2: PDFLIB proton s.f., with MSTP(51) = 1000xNGROUP+NSET
	446	// If the following mass number both not equal zero, nuclear corrections of the stf are used.
	447	// MSTP(192) : Mass number of nucleus side 1
	448	// MSTP(193) : Mass number of nucleus side 2
	449	SetMSTP(52,2);
	450	SetMSTP(192, a1);
	451	SetMSTP(193, a2);
	452	}
	453
	454
	455	AliPythia* AliPythia::Instance()
	456	{
	457	// Set random number generator
	458	if (fgAliPythia) {
	459	return fgAliPythia;
	460	} else {
	461	fgAliPythia = new AliPythia();
	462	return fgAliPythia;
	463	}
	464	}
	465
	466	void AliPythia::PrintParticles()
	467	{
	468	// Print list of particl properties
	469	Int_t np = 0;
c31f1d37	470	char* name = new char[16];
8d2cd130	471	for (Int_t kf=0; kf<1000000; kf++) {
8d2cd130	472	for (Int_t c = 1; c > -2; c-=2) {
8d2cd130	473	Int_t kc = Pycomp(c*kf);
	474	if (kc) {
	475	Float_t mass = GetPMAS(kc,1);
	476	Float_t width = GetPMAS(kc,2);
	477	Float_t tau = GetPMAS(kc,4);
c31f1d37	478
8d2cd130	479	Pyname(kf,name);
	480
	481	np++;
	482
	483	printf("\n mass, width, tau: %6d %s %10.3f %10.3e %10.3e",
	484	c*kf, name, mass, width, tau);
	485	}
	486	}
	487	}
	488	printf("\n Number of particles %d \n \n", np);
	489	}
	490
	491	void AliPythia::ResetDecayTable()
	492	{
	493	// Set default values for pythia decay switches
	494	Int_t i;
	495	for (i = 1; i < 501; i++) SetMDCY(i,1,fDefMDCY[i]);
	496	for (i = 1; i < 2001; i++) SetMDME(i,1,fDefMDME[i]);
	497	}
	498
	499	void AliPythia::SetDecayTable()
	500	{
	501	// Set default values for pythia decay switches
	502	//
	503	Int_t i;
	504	for (i = 1; i < 501; i++) fDefMDCY[i] = GetMDCY(i,1);
	505	for (i = 1; i < 2001; i++) fDefMDME[i] = GetMDME(i,1);
	506	}
	507
	508	void AliPythia::Pyclus(Int_t& njet)
	509	{
	510	// Call Pythia clustering algorithm
	511	//
	512	pyclus(njet);
	513	}
	514
	515	void AliPythia::Pycell(Int_t& njet)
	516	{
	517	// Call Pythia jet reconstruction algorithm
	518	//
	519	pycell(njet);
	520	}
	521
452af8c7	522	void AliPythia::Pyshow(Int_t ip1, Int_t ip2, Double_t qmax)
	523	{
	524	// Call Pythia jet reconstruction algorithm
	525	//
	526	Int_t numpart = fPyjets->N;
	527	for (Int_t i = 0; i < numpart; i++)
	528	{
	529	if (fPyjets->K[2][i] == 7) ip1 = i+1;
	530	if (fPyjets->K[2][i] == 8) ip2 = i+1;
	531	}
	532
	533
	534	qmax = 2. * GetVINT(51);
	535	printf("Pyshow %d %d %f", ip1, ip2, qmax);
	536
	537	pyshow(ip1, ip2, qmax);
	538	}
	539
	540	void AliPythia::Pyrobo(Int_t imi, Int_t ima, Double_t the, Double_t phi, Double_t bex, Double_t bey, Double_t bez)
	541	{
	542	pyrobo(imi, ima, the, phi, bex, bey, bez);
	543	}
	544
	545
	546
	547	void AliPythia::Quench()
	548	{
	549	//
	550	//
	551	// Simple Jet Quenching routine:
	552	// =============================
	553	// The jet formed by all final state partons radiated by the parton created
	554	// in the hard collisions is quenched by a factor z using:
	555	// (E + p_z)new = (1-z) (E + p_z)old
	556	//
	557	// The lost momentum is first balanced by one gluon with virtuality > 0.
	558	// Subsequently the gluon splits to yield two gluons with E = p.
	559	//
	560	Float_t p0[2][5];
	561	Float_t p1[2][5];
	562	Float_t p2[2][5];
	563	Int_t klast[2] = {-1, -1};
	564	Int_t kglu[2];
	565	for (Int_t i = 0; i < 4; i++)
	566	{
	567	p0[0][i] = 0.;
	568	p0[1][i] = 0.;
	569	p1[0][i] = 0.;
	570	p1[1][i] = 0.;
	571	p2[0][i] = 0.;
	572	p2[1][i] = 0.;
	573	}
	574
	575	Int_t numpart = fPyjets->N;
	576
	577	for (Int_t i = 0; i < numpart; i++)
	578	{
	579	Int_t imo = fPyjets->K[2][i];
	580	Int_t kst = fPyjets->K[0][i];
	581	Int_t pdg = fPyjets->K[1][i];
	582
	583	// Quarks and gluons only
	584	if (pdg != 21 && TMath::Abs(pdg) > 6) continue;
	585
586	// Particles from hard scattering only
587
588
589	Float_t px = fPyjets->P[0][i];
590	Float_t py = fPyjets->P[1][i];
591	Float_t pz = fPyjets->P[2][i];
592	Float_t e = fPyjets->P[3][i];
593	Float_t m = fPyjets->P[4][i];
594	Float_t pt = TMath::Sqrt(px * px + py * py);
595	// Skip comment lines
596	if (kst != 1 && kst != 2) continue;
597
598	Float_t mt = TMath::Sqrt(px * px + py * py + m * m);
599
600	//
601	// Some cuts to be in a save kinematic region
602	//
603	if (imo != 7 && imo != 8) continue;
604	Int_t index = imo - 7;
605	klast[index] = i;
606
607	p0[index][0] += px;
608	p0[index][1] += py;
609	p0[index][2] += pz;
511db649	610	p0[index][3] += e;
	611	//
	612	// Fix z
	613	//
452af8c7	614
452af8c7	615	Float_t z = 0.2;
	616	Float_t eppzOld = e + pz;
	617	Float_t empzOld = e - pz;
	618
	619
	620	//
	621	// Kinematics of the original parton
	622	//
	623
	624	Float_t eppzNew = (1. - z) * eppzOld;
	625	Float_t empzNew = empzOld - mt * mt * z / eppzOld;
	626	Float_t eNew0 = 0.5 * (eppzNew + empzNew);
	627	Float_t pzNew0 = 0.5 * (eppzNew - empzNew);
	628	//
	629	// Skip if pt too small
	630	//
	631	if (m * m > eppzNew * empzNew) continue;
	632	Float_t ptNew = TMath::Sqrt(eppzNew * empzNew - m * m);
	633	Float_t pxNew0 = ptNew / pt * px;
	634	Float_t pyNew0 = ptNew / pt * py;
	635
	636	p1[index][0] += pxNew0;
	637	p1[index][1] += pyNew0;
	638	p1[index][2] += pzNew0;
	639	p1[index][3] += eNew0;
	640	//
	641	// Update event record
	642	//
	643	fPyjets->P[0][i] = pxNew0;
	644	fPyjets->P[1][i] = pyNew0;
	645	fPyjets->P[2][i] = pzNew0;
	646	fPyjets->P[3][i] = eNew0;
	647
	648	}
	649
	650	//
	651	// Gluons
	652	//
	653
	654	for (Int_t k = 0; k < 2; k++)
	655	{
	656	for (Int_t j = 0; j < 4; j++)
	657	{
	658	p2[k][j] = p0[k][j] - p1[k][j];
	659	}
	660	p2[k][4] = p2[k][3] * p2[k][3] - p2[k][0] * p2[k][0] - p2[k][1] * p2[k][1] - p2[k][2] * p2[k][2];
	661
	662	if (p2[k][4] > 0.)
	663	{
	664
	665	//
	666	// Bring gluon back to mass shell via momentum scaling
	667	// (momentum will not be conserved, but energy)
	668	//
	669	// not used anymore
	670	/*
	671	Float_t psq = p2[k][0] * p2[k][0] + p2[k][1] * p2[k][1] + p2[k][2] * p2[k][2];
	672	Float_t fact = TMath::Sqrt(1. + p2[k][4] / psq);
	673	p2[k][0] *= fact;
	674	p2[k][1] *= fact;
	675	p2[k][2] *= fact;
	676	p2[k][3] = TMath::Sqrt(psq) * fact;
	677	p2[k][4] = 0.;
	678	*/
679	}
680	}
681
682	if (p2[0][4] > 0.) {
683	p2[0][4] = TMath::Sqrt(p2[0][4]);
684	} else {
685	printf("Warning negative mass squared ! \n");
686	}
687
688	if (p2[1][4] > 0.) {
689	p2[1][4] = TMath::Sqrt(p2[1][4]);
690	} else {
691	printf("Warning negative mass squared ! \n");
692	}
693
694	//
695	// Add the gluons
696	//
697
698
699	for (Int_t i = 0; i < 2; i++) {
700	Int_t ish, jmin, jmax, iGlu, iNew;
701	Int_t in = klast[i];
702	ish = 0;
703
704	if (in == -1) continue;
705	if (i == 1 && klast[1] > klast[0]) in += ish;
706
707	jmin = in - 1;
708	ish = 1;
709
710	if (p2[i][4] > 0) ish = 2;
711
712	iGlu = in;
713	iNew = in + ish;
714	jmax = numpart + ish - 1;
715
716	if (fPyjets->K[0][in-1] == 1 \|\| fPyjets->K[0][in-1] == 21 \|\| fPyjets->K[0][in-1] == 11) {
717	jmin = in;
718	iGlu = in + 1;
719	iNew = in;
720	}
721
722	kglu[i] = iGlu;
723
724	for (Int_t j = jmax; j > jmin; j--)
725	{
726	for (Int_t k = 0; k < 5; k++) {
727	fPyjets->K[k][j] = fPyjets->K[k][j-ish];
728	fPyjets->P[k][j] = fPyjets->P[k][j-ish];
729	fPyjets->V[k][j] = fPyjets->V[k][j-ish];
730	}
731	} // end shifting
732	numpart += ish;
733	(fPyjets->N) += ish;
734
735	if (ish == 1) {
736	fPyjets->P[0][iGlu] = p2[i][0];
737	fPyjets->P[1][iGlu] = p2[i][1];
738	fPyjets->P[2][iGlu] = p2[i][2];
739	fPyjets->P[3][iGlu] = p2[i][3];
740	fPyjets->P[4][iGlu] = p2[i][4];
741
742	fPyjets->K[0][iGlu] = 2;
743	fPyjets->K[1][iGlu] = 21;
744	fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
745	fPyjets->K[3][iGlu] = -1;
746	fPyjets->K[4][iGlu] = -1;
747	} else {
748	//
749	// Split gluon in rest frame.
750	//
751	Double_t bx = p2[i][0] / p2[i][3];
752	Double_t by = p2[i][1] / p2[i][3];
753	Double_t bz = p2[i][2] / p2[i][3];
754
755	Float_t pst = p2[i][4] / 2.;
756
757	Float_t cost = 2. * gRandom->Rndm() - 1.;
758	Float_t sint = TMath::Sqrt(1. - cost * cost);
759	Float_t phi = 2. * TMath::Pi() * gRandom->Rndm();
760
761	Float_t pz1 = pst * cost;
762	Float_t pz2 = -pst * cost;
763	Float_t pt1 = pst * sint;
764	Float_t pt2 = -pst * sint;
765	Float_t px1 = pt1 * TMath::Cos(phi);
766	Float_t py1 = pt1 * TMath::Sin(phi);
767	Float_t px2 = pt2 * TMath::Cos(phi);
768	Float_t py2 = pt2 * TMath::Sin(phi);
769
770	fPyjets->P[0][iGlu] = px1;
771	fPyjets->P[1][iGlu] = py1;
772	fPyjets->P[2][iGlu] = pz1;
773	fPyjets->P[3][iGlu] = pst;
774	fPyjets->P[4][iGlu] = 0.;
775
776	fPyjets->K[0][iGlu] = 2;
777	fPyjets->K[1][iGlu] = 21;
778	fPyjets->K[2][iGlu] = fPyjets->K[2][iNew];
779	fPyjets->K[3][iGlu] = -1;
780	fPyjets->K[4][iGlu] = -1;
781
782	fPyjets->P[0][iGlu+1] = px2;
783	fPyjets->P[1][iGlu+1] = py2;
784	fPyjets->P[2][iGlu+1] = pz2;
785	fPyjets->P[3][iGlu+1] = pst;
786	fPyjets->P[4][iGlu+1] = 0.;
787
788	fPyjets->K[0][iGlu+1] = 2;
789	fPyjets->K[1][iGlu+1] = 21;
790	fPyjets->K[2][iGlu+1] = fPyjets->K[2][iNew];
791	fPyjets->K[3][iGlu+1] = -1;
792	fPyjets->K[4][iGlu+1] = -1;
793	SetMSTU(1,0);
794	SetMSTU(2,0);
795
796	//
797	// Boost back
798	//
799	Pyrobo(iGlu + 1, iGlu + 2, 0., 0., bx, by, bz);
800
801	}
802	} // end adding gluons
803	} // end quench
804
805