[u/mrichter/AliRoot.git] / PYTHIA6 / AliGenPythiaJets.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$ */

//
// Generator using the TPythia interface (via AliPythia)
// to generate jets in pp collisions.
// Using SetNuclei() also nuclear modifications to the structure functions
// can be taken into account. 
// Using SetQuenchingFactor(f) quenched jets can be modelled by superimposing 
// two jets with energies e * f and e * (1-f)
//
// andreas.morsch@cern.ch
//


#include "AliGenPythiaJets.h"
#include "AliRun.h"
#include <TParticle.h>

ClassImp(AliGenPythiaJets)

AliGenPythiaJets::AliGenPythiaJets()
                 :AliGenPythia()
{
// Default Constructor
}

AliGenPythiaJets::AliGenPythiaJets(Int_t npart)
                 :AliGenPythia(npart)
{
    fName = "PythiaJets";
    fTitle= "Jet Generator using PYTHIA";
}

AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia):
    AliGenPythia(Pythia)
{
// copy constructor
    Pythia.Copy(*this);
}

AliGenPythiaJets::~AliGenPythiaJets()
{
// Destructor
}

void AliGenPythiaJets::Init()
{
// Initialization
//
    printf("AliGenPythiaJets::Init() \n");
    
    AliGenPythia::Init();
    
    if (fQuench > 0.) {
	 fEtMinJetQ[0]  = fEtMinJet  * fQuench;  
	 fEtMaxJetQ[0]  = fEtMaxJet  * fQuench;  
	 fEtMinJetQ[1]  = fEtMinJet  * (1. - fQuench);  
	 fEtMaxJetQ[1]  = fEtMaxJet  * (1. - fQuench);
	 fPtHardMinQ[0] = fPtHardMin * fQuench;  
	 fPtHardMaxQ[0] = fPtHardMax * fQuench;  
	 fPtHardMinQ[1] = fPtHardMin * (1. - fQuench);  
	 fPtHardMaxQ[1] = fPtHardMax * (1. - fQuench);  
    }
}

void AliGenPythiaJets::Generate()
{
// Generate one event
    fDecayer->ForceDecay();

    Float_t polar[3]   =   {0,0,0};
    Float_t origin[3]  =   {0,0,0};
    Float_t p[3];
//  converts from mm/c to s
    const Float_t kconv=0.001/2.999792458e8;
//
    Int_t nt  = 0;
    Int_t nc  = 0;
    Int_t jev = 0;
    Int_t j, kf, iparent;
    fTrials=0;
//
//  Set collision vertex position 
    if(fVertexSmear==kPerEvent) {
	fPythia->SetMSTP(151,1);
	for (j=0;j<3;j++) {
	    fPythia->SetPARP(151+j, fOsigma[j]*10.);
	}
    } else if (fVertexSmear==kPerTrack) {
	fPythia->SetMSTP(151,0);
    }
//  Event loop    
    while(1)
    {
	if (fQuench > 0.) {
	    fPythia->SetCKIN(3,fPtHardMinQ[jev]);
	    fPythia->SetCKIN(4,fPtHardMaxQ[jev]);
	    fEtMinJet = fEtMinJetQ[jev];
	    fEtMaxJet = fEtMaxJetQ[jev];	    
	}
	
	fPythia->Pyevnt();
	if (gAlice->GetEvNumber()>=fDebugEventFirst &&
	    gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
	fTrials++;
        //
	// Has this jet triggered
        //
	if ((fEtMinJet != -1) && ! CheckTrigger()) continue;
//
	fPythia->ImportParticles(fParticles,"All");
	Int_t i;
	Int_t np = fParticles->GetEntriesFast();
	if (np == 0 ) continue;
// Get event vertex and discard the event if the z coord. is too big	
	TParticle *iparticle = (TParticle *) fParticles->At(0);
	Float_t distz = iparticle->Vz()/10.;
	if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
//
//
	fVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
	fVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
	fVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);

	Int_t* pParent = new Int_t[np];
	for (i=0; i< np; i++) pParent[i] = -1;
	

	//
	for (i = 0; i<np; i++) {
	    Int_t trackIt = 0;
	    TParticle *  iparticle = (TParticle *) fParticles->At(i);
	    kf = CheckPDGCode(iparticle->GetPdgCode());
	    Int_t ks = iparticle->GetStatusCode();
	    Int_t km = iparticle->GetFirstMother();
	    if ((ks == 1  && kf !=0 && KinematicSelection(iparticle, 0)) ||
		(ks != 1) ||
		(fProcess == kPyJets && ks == 21 && km == 0 && i > 1)) {
		nc++;
		if (ks == 1) trackIt = 1;
		Int_t ipa = iparticle->GetFirstMother() - 1;
		iparent = (ipa > -1) ? pParent[ipa] : -1;
//
// Store track information
//
		p[0] = iparticle->Px();
		p[1] = iparticle->Py();
		p[2] = iparticle->Pz();
		origin[0] = fOrigin[0]+iparticle->Vx()/10.;
		origin[1] = fOrigin[1]+iparticle->Vy()/10.;
		origin[2] = fOrigin[2]+iparticle->Vz()/10.;
		Float_t tof=kconv*iparticle->T();
		PushTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
			 tof, kPPrimary, nt, 1., ks);
		KeepTrack(nt);
		pParent[i] = nt;
	    } // select particle
	} // particle loop 
	
	if (pParent) delete[] pParent;
	printf("\n AliGenPythiaJets: I've put %i particles on the stack \n",nc);
	if (nc > 0) {
	    jev += 1;
	    if ((fQuench <= 0.) || (fQuench > 0. && jev == 2)) {
		fKineBias=Float_t(fNpart)/Float_t(fTrials);
		printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
		fNev++;
		MakeHeader();
		break;
	    }
	}
    }
    SetHighWaterMark(nt);
//  Get cross-section
    fXsection=fPythia->GetPARI(1);
}

Bool_t AliGenPythiaJets::CheckTrigger()
{
// Check the kinematic trigger condition
//
    Bool_t   triggered = kFALSE;
    Int_t njets = 0;
    Int_t ntrig = 0;
    Float_t jets[4][10];
//
// Use Pythia clustering on parton level to determine jet axis
//
    GetJets(njets, ntrig, jets);
    
    if (ntrig) {
	triggered = kTRUE;
	Float_t px   = jets[0][0];
	Float_t py   = jets[1][0];
	Float_t pz   = jets[2][0];
	Float_t e    = jets[3][0];
	Float_t beta = pz/e;
	Float_t phi  =  TMath::Pi()+TMath::ATan2(-py,-px);
	TransformEvent(beta, -2. * TMath::Pi() / 3. + phi);
    }
    return triggered;
}
	  
AliGenPythiaJets& AliGenPythiaJets::operator=(const  AliGenPythiaJets& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return *this;
}

void  AliGenPythiaJets::TransformEvent(Float_t beta, Float_t phi)
{
//
// Perform Lorentz Transformation and Rotation
//
    Float_t gamma = 1./TMath::Sqrt(1. - beta * beta);
    Int_t npart = (fPythia->GetPyjets())->N;

    for (Int_t part = 0; part < npart; part++) {
	Float_t px   = 	(fPythia->GetPyjets())->P[0][part];
	Float_t py   =  (fPythia->GetPyjets())->P[1][part];
	Float_t pz   =  (fPythia->GetPyjets())->P[2][part];
	Float_t e    =  (fPythia->GetPyjets())->P[3][part];
	//
	// Lorentz Transform
	//
	Float_t pzt =  gamma * pz        - gamma * beta * e;
	Float_t et  = -gamma * beta * pz + gamma        * e;
	//
	// Rotation
	//
	Float_t pxt =   TMath::Cos(phi) * px +  TMath::Sin(phi) * py;
	Float_t pyt = - TMath::Sin(phi) * px +  TMath::Cos(phi) * py;
	//
	//
	(fPythia->GetPyjets())->P[0][part] = pxt;
	(fPythia->GetPyjets())->P[1][part] = pyt;
	(fPythia->GetPyjets())->P[2][part] = pzt;
	(fPythia->GetPyjets())->P[3][part] = et;
    }
}
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
88cb7938	16	/* $Id$ */
8d2cd130	17
	18	//
	19	// Generator using the TPythia interface (via AliPythia)
	20	// to generate jets in pp collisions.
	21	// Using SetNuclei() also nuclear modifications to the structure functions
	22	// can be taken into account.
	23	// Using SetQuenchingFactor(f) quenched jets can be modelled by superimposing
	24	// two jets with energies e * f and e * (1-f)
	25	//
	26	// andreas.morsch@cern.ch
	27	//
	28
	29
	30	#include "AliGenPythiaJets.h"
	31	#include "AliRun.h"
	32	#include <TParticle.h>
	33
	34	ClassImp(AliGenPythiaJets)
	35
	36	AliGenPythiaJets::AliGenPythiaJets()
	37	:AliGenPythia()
	38	{
	39	// Default Constructor
	40	}
	41
	42	AliGenPythiaJets::AliGenPythiaJets(Int_t npart)
	43	:AliGenPythia(npart)
	44	{
	45	fName = "PythiaJets";
	46	fTitle= "Jet Generator using PYTHIA";
	47	}
	48
014a9521	49	AliGenPythiaJets::AliGenPythiaJets(const AliGenPythiaJets & Pythia):
014a9521	50	AliGenPythia(Pythia)
8d2cd130	51	{
	52	// copy constructor
	53	Pythia.Copy(*this);
	54	}
	55
	56	AliGenPythiaJets::~AliGenPythiaJets()
	57	{
	58	// Destructor
	59	}
	60
	61	void AliGenPythiaJets::Init()
	62	{
	63	// Initialization
	64	//
	65	printf("AliGenPythiaJets::Init() \n");
	66
	67	AliGenPythia::Init();
	68
	69	if (fQuench > 0.) {
	70	fEtMinJetQ[0] = fEtMinJet * fQuench;
	71	fEtMaxJetQ[0] = fEtMaxJet * fQuench;
	72	fEtMinJetQ[1] = fEtMinJet * (1. - fQuench);
	73	fEtMaxJetQ[1] = fEtMaxJet * (1. - fQuench);
	74	fPtHardMinQ[0] = fPtHardMin * fQuench;
	75	fPtHardMaxQ[0] = fPtHardMax * fQuench;
	76	fPtHardMinQ[1] = fPtHardMin * (1. - fQuench);
	77	fPtHardMaxQ[1] = fPtHardMax * (1. - fQuench);
	78	}
	79	}
	80
	81	void AliGenPythiaJets::Generate()
	82	{
	83	// Generate one event
	84	fDecayer->ForceDecay();
	85
	86	Float_t polar[3] = {0,0,0};
	87	Float_t origin[3] = {0,0,0};
	88	Float_t p[3];
	89	// converts from mm/c to s
	90	const Float_t kconv=0.001/2.999792458e8;
	91	//
	92	Int_t nt = 0;
	93	Int_t nc = 0;
	94	Int_t jev = 0;
	95	Int_t j, kf, iparent;
	96	fTrials=0;
	97	//
	98	// Set collision vertex position
	99	if(fVertexSmear==kPerEvent) {
	100	fPythia->SetMSTP(151,1);
	101	for (j=0;j<3;j++) {
	102	fPythia->SetPARP(151+j, fOsigma[j]*10.);
	103	}
	104	} else if (fVertexSmear==kPerTrack) {
	105	fPythia->SetMSTP(151,0);
	106	}
	107	// Event loop
	108	while(1)
	109	{
	110	if (fQuench > 0.) {
	111	fPythia->SetCKIN(3,fPtHardMinQ[jev]);
	112	fPythia->SetCKIN(4,fPtHardMaxQ[jev]);
	113	fEtMinJet = fEtMinJetQ[jev];
	114	fEtMaxJet = fEtMaxJetQ[jev];
115	}
116
117	fPythia->Pyevnt();
118	if (gAlice->GetEvNumber()>=fDebugEventFirst &&
119	gAlice->GetEvNumber()<=fDebugEventLast) fPythia->Pylist(1);
120	fTrials++;
121	//
122	// Has this jet triggered
123	//
124	if ((fEtMinJet != -1) && ! CheckTrigger()) continue;
125	//
126	fPythia->ImportParticles(fParticles,"All");
127	Int_t i;
128	Int_t np = fParticles->GetEntriesFast();
129	if (np == 0 ) continue;
130	// Get event vertex and discard the event if the z coord. is too big
131	TParticle iparticle = (TParticle ) fParticles->At(0);
132	Float_t distz = iparticle->Vz()/10.;
133	if(TMath::Abs(distz)>fCutVertexZ*fOsigma[2]) continue;
134	//
135	//
d25cfd65	136	fVertex[0] = iparticle->Vx()/10.+fOrigin.At(0);
	137	fVertex[1] = iparticle->Vy()/10.+fOrigin.At(1);
	138	fVertex[2] = iparticle->Vz()/10.+fOrigin.At(2);
8d2cd130	139
	140	Int_t* pParent = new Int_t[np];
	141	for (i=0; i< np; i++) pParent[i] = -1;
	142
	143
	144	//
	145	for (i = 0; i<np; i++) {
	146	Int_t trackIt = 0;
	147	TParticle * iparticle = (TParticle *) fParticles->At(i);
	148	kf = CheckPDGCode(iparticle->GetPdgCode());
	149	Int_t ks = iparticle->GetStatusCode();
	150	Int_t km = iparticle->GetFirstMother();
	151	if ((ks == 1 && kf !=0 && KinematicSelection(iparticle, 0)) \|\|
	152	(ks != 1) \|\|
	153	(fProcess == kPyJets && ks == 21 && km == 0 && i > 1)) {
	154	nc++;
	155	if (ks == 1) trackIt = 1;
	156	Int_t ipa = iparticle->GetFirstMother() - 1;
	157	iparent = (ipa > -1) ? pParent[ipa] : -1;
	158	//
	159	// Store track information
	160	//
	161	p[0] = iparticle->Px();
	162	p[1] = iparticle->Py();
	163	p[2] = iparticle->Pz();
	164	origin[0] = fOrigin[0]+iparticle->Vx()/10.;
	165	origin[1] = fOrigin[1]+iparticle->Vy()/10.;
	166	origin[2] = fOrigin[2]+iparticle->Vz()/10.;
	167	Float_t tof=kconv*iparticle->T();
642f15cf	168	PushTrack(fTrackIt*trackIt, iparent, kf, p, origin, polar,
8d2cd130	169	tof, kPPrimary, nt, 1., ks);
	170	KeepTrack(nt);
	171	pParent[i] = nt;
	172	} // select particle
	173	} // particle loop
	174
	175	if (pParent) delete[] pParent;
	176	printf("\n AliGenPythiaJets: I've put %i particles on the stack \n",nc);
	177	if (nc > 0) {
	178	jev += 1;
	179	if ((fQuench <= 0.) \|\| (fQuench > 0. && jev == 2)) {
	180	fKineBias=Float_t(fNpart)/Float_t(fTrials);
	181	printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
	182	fNev++;
	183	MakeHeader();
	184	break;
	185	}
	186	}
	187	}
	188	SetHighWaterMark(nt);
	189	// Get cross-section
	190	fXsection=fPythia->GetPARI(1);
	191	}
	192
	193	Bool_t AliGenPythiaJets::CheckTrigger()
	194	{
	195	// Check the kinematic trigger condition
	196	//
	197	Bool_t triggered = kFALSE;
	198	Int_t njets = 0;
	199	Int_t ntrig = 0;
	200	Float_t jets[4][10];
	201	//
	202	// Use Pythia clustering on parton level to determine jet axis
	203	//
	204	GetJets(njets, ntrig, jets);
	205
	206	if (ntrig) {
	207	triggered = kTRUE;
	208	Float_t px = jets[0][0];
	209	Float_t py = jets[1][0];
	210	Float_t pz = jets[2][0];
	211	Float_t e = jets[3][0];
	212	Float_t beta = pz/e;
	213	Float_t phi = TMath::Pi()+TMath::ATan2(-py,-px);
	214	TransformEvent(beta, -2. * TMath::Pi() / 3. + phi);
	215	}
	216	return triggered;
	217	}
	218
	219	AliGenPythiaJets& AliGenPythiaJets::operator=(const AliGenPythiaJets& rhs)
	220	{
	221	// Assignment operator
014a9521	222	rhs.Copy(*this);
8d2cd130	223	return *this;
	224	}
	225
	226	void AliGenPythiaJets::TransformEvent(Float_t beta, Float_t phi)
	227	{
	228	//
	229	// Perform Lorentz Transformation and Rotation
	230	//
	231	Float_t gamma = 1./TMath::Sqrt(1. - beta * beta);
	232	Int_t npart = (fPythia->GetPyjets())->N;
	233
	234	for (Int_t part = 0; part < npart; part++) {
	235	Float_t px = (fPythia->GetPyjets())->P[0][part];
	236	Float_t py = (fPythia->GetPyjets())->P[1][part];
	237	Float_t pz = (fPythia->GetPyjets())->P[2][part];
	238	Float_t e = (fPythia->GetPyjets())->P[3][part];
	239	//
	240	// Lorentz Transform
	241	//
	242	Float_t pzt = gamma * pz - gamma * beta * e;
	243	Float_t et = -gamma * beta * pz + gamma * e;
	244	//
	245	// Rotation
	246	//
	247	Float_t pxt = TMath::Cos(phi) * px + TMath::Sin(phi) * py;
	248	Float_t pyt = - TMath::Sin(phi) * px + TMath::Cos(phi) * py;
	249	//
	250	//
	251	(fPythia->GetPyjets())->P[0][part] = pxt;
	252	(fPythia->GetPyjets())->P[1][part] = pyt;
	253	(fPythia->GetPyjets())->P[2][part] = pzt;
	254	(fPythia->GetPyjets())->P[3][part] = et;
	255	}
	256	}
	257