[u/mrichter/AliRoot.git] / EVGEN / AliGenHijing.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.                  *
 **************************************************************************/

/*
$Log$
Revision 1.7  2000/10/02 21:28:06  fca
Removal of useless dependecies via forward declarations

Revision 1.6  2000/09/11 13:23:37  morsch
Write last seed to file (fortran lun 50) and reed back from same lun using calls to
luget_hijing and luset_hijing.

Revision 1.5  2000/09/07 16:55:40  morsch
fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)

Revision 1.4  2000/07/11 18:24:56  fca
Coding convention corrections + few minor bug fixes

Revision 1.3  2000/06/30 12:08:36  morsch
In member data: char* replaced by TString, Init takes care of resizing the strings to
8 characters required by Hijing.

Revision 1.2  2000/06/15 14:15:05  morsch
Add possibility for heavy flavor selection: charm and beauty.

Revision 1.1  2000/06/09 20:47:27  morsch
AliGenerator interface class to HIJING using THijing (test version)

*/

#include "AliGenHijing.h"
#include "AliGenHijingEventHeader.h"
#include "AliRun.h"
#include "AliMC.h"

#include <TArrayI.h>
#include <TParticle.h>
#include <THijing.h>

 ClassImp(AliGenHijing)

AliGenHijing::AliGenHijing()
                 :AliGenerator()
{
// Constructor
}

AliGenHijing::AliGenHijing(Int_t npart)
                 :AliGenerator(npart)
{
// Default PbPb collisions at 5. 5 TeV
//
    SetEnergyCMS();
    SetImpactParameterRange();
    SetTarget();
    SetProjectile();
    fKeep=0;
    fQuench=1;
    fShadowing=1;
    fTrigger=0;
    fDecaysOff=1;
    fEvaluate=0;
    fSelectAll=0;
    fFlavor=0;
}

AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
{
// copy constructor
}


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

void AliGenHijing::Init()
{
// Initialisation
    fFrame.Resize(8);
    fTarget.Resize(8);
    fProjectile.Resize(8);
    
    SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget, 
		      fAProjectile, fZProjectile, fATarget, fZTarget, 
		      fMinImpactParam, fMaxImpactParam));

    fHijing=(THijing*) fgMCEvGen;

    fHijing->SetIHPR2(3,  fTrigger);
    fHijing->SetIHPR2(4,  fQuench);
    fHijing->SetIHPR2(6,  fShadowing);
    fHijing->SetIHPR2(12, fDecaysOff);    
    fHijing->SetIHPR2(21, fKeep);
    fHijing->Rluset(50,0);
    fHijing->Initialize();

    
//
    if (fEvaluate) EvaluateCrossSections();
}

void AliGenHijing::Generate()
{
// Generate one event

    Float_t polar[3] =   {0,0,0};
    Float_t origin[3]=   {0,0,0};
    Float_t origin0[3]=  {0,0,0};
    Float_t p[3], random[6];
    Float_t tof;

    static TClonesArray *particles;
//  converts from mm/c to s
    const Float_t kconv=0.001/2.999792458e8;
//
    Int_t nt=0;
    Int_t jev=0;
    Int_t j, kf, ks, imo;
    kf=0;
    
    if(!particles) particles=new TClonesArray("TParticle",10000);
    
    fTrials=0;
    for (j=0;j<3;j++) origin0[j]=fOrigin[j];
    if(fVertexSmear==kPerEvent) {
	gMC->Rndm(random,6);
	for (j=0;j<3;j++) {
	    origin0[j]+=fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
		TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
//	    fHijing->SetMSTP(151,0);
	}
    } else if (fVertexSmear==kPerTrack) {
//	fHijing->SetMSTP(151,0);
	for (j=0;j<3;j++) {
//	    fHijing->SetPARP(151+j, fOsigma[j]*10.);
	}
    }
    while(1)
    {

	fHijing->GenerateEvent();
	fTrials++;
	fHijing->ImportParticles(particles,"All");
	Int_t np = particles->GetEntriesFast();
	printf("\n **************************************************%d\n",np);
	Int_t nc=0;
	if (np == 0 ) continue;
	Int_t i;
	Int_t * newPos = new Int_t[np];
	for (i = 0; i<np; i++) *(newPos+i)=i;
	
	for (i = 0; i<np; i++) {
	    TParticle *  iparticle       = (TParticle *) particles->At(i);

	    Bool_t  hasMother            =  (iparticle->GetFirstMother()   >=0);
	    Bool_t  hasDaughter          =  (iparticle->GetFirstDaughter() >=0);
	    Bool_t  selected             =  kTRUE;
	    Bool_t  hasSelectedDaughters =  kFALSE;


	    kf        = iparticle->GetPdgCode();
	    if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
	    if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
//
// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
//

	    if (selected || hasSelectedDaughters) {
		nc++;
		ks        = iparticle->GetStatusCode();
		p[0]=iparticle->Px();
		p[1]=iparticle->Py();
		p[2]=iparticle->Pz();
		origin[0]=origin0[0]+iparticle->Vx()/10;
		origin[1]=origin0[1]+iparticle->Vy()/10;
		origin[2]=origin0[2]+iparticle->Vz()/10;
		tof=kconv*iparticle->T();
		imo=-1;
		if (hasMother) {
		    imo=iparticle->GetFirstMother();
		    imo=*(newPos+imo);
		}
		
//		printf("\n selected iparent %d %d %d \n",i, kf, imo);
		if (hasDaughter) {
		    gAlice->SetTrack(0,imo,kf,p,origin,polar,
				     tof,"Primary",nt);
		} else {
		    gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
				     tof,"Secondary",nt);
		}
		*(newPos+i)=nt;
	    } // selected
	} // particle loop 
	delete newPos;
	printf("\n I've put %i particles on the stack \n",nc);
	if (nc > 0) {
	    jev+=nc;
	    if (jev >= fNpart || fNpart == -1) {
		fKineBias=Float_t(fNpart)/Float_t(fTrials);
		printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
		break;
	    }
	}
    } // event loop
    fHijing->Rluget(50,-1);
}

Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
{
// Perform kinematic selection
    Float_t px=particle->Px();
    Float_t py=particle->Py();
    Float_t pz=particle->Pz();
    Float_t  e=particle->Energy();

//
//  transverse momentum cut    
    Float_t pt=TMath::Sqrt(px*px+py*py);
    if (pt > fPtMax || pt < fPtMin) 
    {
//	printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
	return kFALSE;
    }
//
// momentum cut
    Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
    if (p > fPMax || p < fPMin) 
    {
//	printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
	return kFALSE;
    }
    
//
// theta cut
    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
    if (theta > fThetaMax || theta < fThetaMin) 
    {
	
//    	printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
	return kFALSE;
    }

//
// rapidity cut
    Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
    if (y > fYMax || y < fYMin)
    {
//	printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
	return kFALSE;
    }

//
// phi cut
    Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
    if (phi > fPhiMax || phi < fPhiMin)
    {
//	printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
	return kFALSE;
    }

    return kTRUE;
}

void AliGenHijing::KeepFullEvent()
{
    fKeep=1;
}

void AliGenHijing::EvaluateCrossSections()
{
//     Glauber Calculation of geometrical x-section
//
    Float_t xTot=0.;          // barn
    Float_t xTotHard=0.;      // barn 
    Float_t xPart=0.;         // barn
    Float_t xPartHard=0.;     // barn 
    Float_t sigmaHard=0.1;    // mbarn
    Float_t bMin=0.;
    Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
    const Float_t kdib=0.2;
    Int_t   kMax=Int_t((bMax-bMin)/kdib)+1;


    printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
    printf("\n Target     Radius (fm): %f \n",fHijing->GetHIPR1(35));    
    Int_t i;
    Float_t oldvalue=0.;
    
    for (i=0; i<kMax; i++)
    {
	Float_t xb=bMin+i*kdib;
	Float_t ov;
	ov=fHijing->Profile(xb);
	Float_t gb =  2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*(1.-TMath::Exp(-fHijing->GetHINT1(12)*ov));
	Float_t gbh = 2.*0.01*fHijing->GetHIPR1(40)*kdib*xb*sigmaHard*ov;
	xTot+=gb;
	xTotHard+=gbh;
	if (xb > fMinImpactParam && xb < fMaxImpactParam)
	{
	    xPart+=gb;
	    xPartHard+=gbh;
	}
	if (xTot == 0.) continue;
	
	if ((xTot-oldvalue)/xTot<0.0001) break;
	oldvalue=xTot;
	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
	printf("\n Hard  cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
    }
    printf("\n Total cross section (barn): %f \n",xTot);
    printf("\n Hard  cross section (barn): %f \n \n",xTotHard);
    printf("\n Partial       cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
    printf("\n Partial  hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
}

Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
{
//
// Looks recursively if one of the daughters has been selected
//
//    printf("\n Consider daughters %d:",iparticle->GetPdgCode());
    Int_t imin=-1;
    Int_t imax=-1;
    Int_t i;
    Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
    Bool_t selected=kFALSE;
    if (hasDaughters) {
	imin=iparticle->GetFirstDaughter();
	imax=iparticle->GetLastDaughter();       
	for (i=imin; i<= imax; i++){
	    TParticle *  jparticle       = (TParticle *) particles->At(i);	
	    Int_t ip=jparticle->GetPdgCode();
	    if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
	    if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	}
    } else {
	return kFALSE;
    }

    return selected;
}


Bool_t AliGenHijing::SelectFlavor(Int_t pid)
{
// Select flavor of particle
// 0: all
// 4: charm and beauty
// 5: beauty
    if (fFlavor == 0) return kTRUE;
    
    Int_t ifl=TMath::Abs(pid/100);
    if (ifl > 10) ifl/=10;
    return ((fFlavor==4 && (ifl==4 || ifl==5))  || 
	    (fFlavor==5 &&  ifl==5));

}

void AliGenHijing::MakeHeader()
{
// Builds the event header, to be called after each event
    AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
//    header->SetDate(date);
//    header->SetRunNumber(run);
//    header->SetEventNumber(event);
    header->SetNProduced(fHijing->GetNATT());
    header->SetImpactParameter(fHijing->GetHINT1(19));
    header->SetTotalEnergy(fHijing->GetEATT());
    header->SetHardScatters(fHijing->GetJATT());
    header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
    header->SetCollisions(fHijing->GetN0(),
			  fHijing->GetN01(),
			  fHijing->GetN10(),
			  fHijing->GetN11());
}

AliGenHijing& AliGenHijing::operator=(const  AliGenHijing& rhs)
{
// Assignment operator
    return *this;
}
Commit	Line	Data
110410f9	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
	16	/*
	17	$Log$
73bc6a56	18	Revision 1.7 2000/10/02 21:28:06 fca
	19	Removal of useless dependecies via forward declarations
	20
94de3818	21	Revision 1.6 2000/09/11 13:23:37 morsch
	22	Write last seed to file (fortran lun 50) and reed back from same lun using calls to
	23	luget_hijing and luset_hijing.
	24
b78ff279	25	Revision 1.5 2000/09/07 16:55:40 morsch
	26	fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
	27
ba244bac	28	Revision 1.4 2000/07/11 18:24:56 fca
	29	Coding convention corrections + few minor bug fixes
	30
aee8290b	31	Revision 1.3 2000/06/30 12:08:36 morsch
	32	In member data: char* replaced by TString, Init takes care of resizing the strings to
	33	8 characters required by Hijing.
	34
3b8bd63f	35	Revision 1.2 2000/06/15 14:15:05 morsch
	36	Add possibility for heavy flavor selection: charm and beauty.
	37
eb342d2d	38	Revision 1.1 2000/06/09 20:47:27 morsch
	39	AliGenerator interface class to HIJING using THijing (test version)
	40
110410f9	41	*/
	42
	43	#include "AliGenHijing.h"
3b8bd63f	44	#include "AliGenHijingEventHeader.h"
110410f9	45	#include "AliRun.h"
94de3818	46	#include "AliMC.h"
110410f9	47
	48	#include <TArrayI.h>
	49	#include <TParticle.h>
	50	#include <THijing.h>
	51
	52	ClassImp(AliGenHijing)
	53
	54	AliGenHijing::AliGenHijing()
	55	:AliGenerator()
	56	{
	57	// Constructor
	58	}
	59
	60	AliGenHijing::AliGenHijing(Int_t npart)
	61	:AliGenerator(npart)
	62	{
	63	// Default PbPb collisions at 5. 5 TeV
	64	//
	65	SetEnergyCMS();
	66	SetImpactParameterRange();
	67	SetTarget();
	68	SetProjectile();
	69	fKeep=0;
	70	fQuench=1;
	71	fShadowing=1;
	72	fTrigger=0;
	73	fDecaysOff=1;
	74	fEvaluate=0;
	75	fSelectAll=0;
eb342d2d	76	fFlavor=0;
110410f9	77	}
	78
	79	AliGenHijing::AliGenHijing(const AliGenHijing & Hijing)
	80	{
	81	// copy constructor
	82	}
	83
	84
	85	AliGenHijing::~AliGenHijing()
	86	{
	87	// Destructor
	88	}
	89
	90	void AliGenHijing::Init()
	91	{
	92	// Initialisation
3b8bd63f	93	fFrame.Resize(8);
	94	fTarget.Resize(8);
	95	fProjectile.Resize(8);
	96
110410f9	97	SetMC(new THijing(fEnergyCMS, fFrame, fProjectile, fTarget,
	98	fAProjectile, fZProjectile, fATarget, fZTarget,
	99	fMinImpactParam, fMaxImpactParam));
	100
	101	fHijing=(THijing*) fgMCEvGen;
ba244bac	102
110410f9	103	fHijing->SetIHPR2(3, fTrigger);
	104	fHijing->SetIHPR2(4, fQuench);
	105	fHijing->SetIHPR2(6, fShadowing);
	106	fHijing->SetIHPR2(12, fDecaysOff);
	107	fHijing->SetIHPR2(21, fKeep);
b78ff279	108	fHijing->Rluset(50,0);
ba244bac	109	fHijing->Initialize();
b78ff279	110
b78ff279	111
110410f9	112	//
	113	if (fEvaluate) EvaluateCrossSections();
	114	}
	115
	116	void AliGenHijing::Generate()
	117	{
	118	// Generate one event
	119
	120	Float_t polar[3] = {0,0,0};
	121	Float_t origin[3]= {0,0,0};
	122	Float_t origin0[3]= {0,0,0};
	123	Float_t p[3], random[6];
	124	Float_t tof;
	125
	126	static TClonesArray *particles;
	127	// converts from mm/c to s
	128	const Float_t kconv=0.001/2.999792458e8;
	129	//
	130	Int_t nt=0;
	131	Int_t jev=0;
	132	Int_t j, kf, ks, imo;
eb342d2d	133	kf=0;
eb342d2d	134
110410f9	135	if(!particles) particles=new TClonesArray("TParticle",10000);
	136
	137	fTrials=0;
	138	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
aee8290b	139	if(fVertexSmear==kPerEvent) {
110410f9	140	gMC->Rndm(random,6);
	141	for (j=0;j<3;j++) {
	142	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	143	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	144	// fHijing->SetMSTP(151,0);
	145	}
aee8290b	146	} else if (fVertexSmear==kPerTrack) {
110410f9	147	// fHijing->SetMSTP(151,0);
	148	for (j=0;j<3;j++) {
	149	// fHijing->SetPARP(151+j, fOsigma[j]*10.);
	150	}
	151	}
	152	while(1)
	153	{
	154
	155	fHijing->GenerateEvent();
	156	fTrials++;
	157	fHijing->ImportParticles(particles,"All");
	158	Int_t np = particles->GetEntriesFast();
	159	printf("\n **************************************************%d\n",np);
	160	Int_t nc=0;
	161	if (np == 0 ) continue;
	162	Int_t i;
	163	Int_t * newPos = new Int_t[np];
eb342d2d	164	for (i = 0; i<np; i++) *(newPos+i)=i;
110410f9	165
eb342d2d	166	for (i = 0; i<np; i++) {
110410f9	167	TParticle * iparticle = (TParticle *) particles->At(i);
	168
	169	Bool_t hasMother = (iparticle->GetFirstMother() >=0);
	170	Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
	171	Bool_t selected = kTRUE;
	172	Bool_t hasSelectedDaughters = kFALSE;
	173
110410f9	174
eb342d2d	175	kf = iparticle->GetPdgCode();
eb342d2d	176	if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
110410f9	177	if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
	178	//
	179	// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
	180	//
	181
	182	if (selected \|\| hasSelectedDaughters) {
	183	nc++;
	184	ks = iparticle->GetStatusCode();
	185	p[0]=iparticle->Px();
	186	p[1]=iparticle->Py();
	187	p[2]=iparticle->Pz();
	188	origin[0]=origin0[0]+iparticle->Vx()/10;
	189	origin[1]=origin0[1]+iparticle->Vy()/10;
	190	origin[2]=origin0[2]+iparticle->Vz()/10;
	191	tof=kconv*iparticle->T();
	192	imo=-1;
	193	if (hasMother) {
	194	imo=iparticle->GetFirstMother();
	195	imo=*(newPos+imo);
	196	}
	197
	198	// printf("\n selected iparent %d %d %d \n",i, kf, imo);
	199	if (hasDaughter) {
	200	gAlice->SetTrack(0,imo,kf,p,origin,polar,
	201	tof,"Primary",nt);
	202	} else {
	203	gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
	204	tof,"Secondary",nt);
	205	}
	206	*(newPos+i)=nt;
	207	} // selected
	208	} // particle loop
	209	delete newPos;
	210	printf("\n I've put %i particles on the stack \n",nc);
	211	if (nc > 0) {
	212	jev+=nc;
	213	if (jev >= fNpart \|\| fNpart == -1) {
	214	fKineBias=Float_t(fNpart)/Float_t(fTrials);
	215	printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
	216	break;
	217	}
	218	}
	219	} // event loop
b78ff279	220	fHijing->Rluget(50,-1);
110410f9	221	}
	222
	223	Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
	224	{
	225	// Perform kinematic selection
	226	Float_t px=particle->Px();
	227	Float_t py=particle->Py();
	228	Float_t pz=particle->Pz();
	229	Float_t e=particle->Energy();
	230
	231	//
	232	// transverse momentum cut
	233	Float_t pt=TMath::Sqrt(pxpx+pypy);
	234	if (pt > fPtMax \|\| pt < fPtMin)
	235	{
	236	// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
	237	return kFALSE;
	238	}
	239	//
	240	// momentum cut
	241	Float_t p=TMath::Sqrt(pxpx+pypy+pz*pz);
	242	if (p > fPMax \|\| p < fPMin)
	243	{
	244	// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
	245	return kFALSE;
	246	}
	247
	248	//
	249	// theta cut
	250	Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
	251	if (theta > fThetaMax \|\| theta < fThetaMin)
	252	{
	253
	254	// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
	255	return kFALSE;
	256	}
	257
	258	//
	259	// rapidity cut
	260	Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
	261	if (y > fYMax \|\| y < fYMin)
	262	{
	263	// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
	264	return kFALSE;
	265	}
	266
	267	//
	268	// phi cut
	269	Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
	270	if (phi > fPhiMax \|\| phi < fPhiMin)
	271	{
	272	// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
	273	return kFALSE;
	274	}
	275
	276	return kTRUE;
	277	}
	278
	279	void AliGenHijing::KeepFullEvent()
	280	{
	281	fKeep=1;
	282	}
	283
	284	void AliGenHijing::EvaluateCrossSections()
285	{
286	// Glauber Calculation of geometrical x-section
287	//
288	Float_t xTot=0.; // barn
289	Float_t xTotHard=0.; // barn
290	Float_t xPart=0.; // barn
291	Float_t xPartHard=0.; // barn
292	Float_t sigmaHard=0.1; // mbarn
293	Float_t bMin=0.;
294	Float_t bMax=fHijing->GetHIPR1(34)+fHijing->GetHIPR1(35);
295	const Float_t kdib=0.2;
296	Int_t kMax=Int_t((bMax-bMin)/kdib)+1;
297
298
299	printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
300	printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
301	Int_t i;
302	Float_t oldvalue=0.;
303
304	for (i=0; i<kMax; i++)
305	{
306	Float_t xb=bMin+i*kdib;
307	Float_t ov;
308	ov=fHijing->Profile(xb);
309	Float_t gb = 2.0.01fHijing->GetHIPR1(40)kdibxb(1.-TMath::Exp(-fHijing->GetHINT1(12)ov));
310	Float_t gbh = 2.0.01fHijing->GetHIPR1(40)kdibxbsigmaHardov;
311	xTot+=gb;
312	xTotHard+=gbh;
313	if (xb > fMinImpactParam && xb < fMaxImpactParam)
314	{
315	xPart+=gb;
316	xPartHard+=gbh;
317	}
73bc6a56	318	if (xTot == 0.) continue;
110410f9	319
73bc6a56	320	if ((xTot-oldvalue)/xTot<0.0001) break;
110410f9	321	oldvalue=xTot;
	322	printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
	323	printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
	324	}
	325	printf("\n Total cross section (barn): %f \n",xTot);
	326	printf("\n Hard cross section (barn): %f \n \n",xTotHard);
	327	printf("\n Partial cross section (barn): %f %f \n",xPart, xPart/xTot*100.);
	328	printf("\n Partial hard cross section (barn): %f %f \n",xPartHard, xPartHard/xTotHard*100.);
	329	}
	330
	331	Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
	332	{
	333	//
	334	// Looks recursively if one of the daughters has been selected
	335	//
	336	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	337	Int_t imin=-1;
	338	Int_t imax=-1;
	339	Int_t i;
	340	Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
	341	Bool_t selected=kFALSE;
	342	if (hasDaughters) {
	343	imin=iparticle->GetFirstDaughter();
	344	imax=iparticle->GetLastDaughter();
110410f9	345	for (i=imin; i<= imax; i++){
110410f9	346	TParticle * jparticle = (TParticle *) particles->At(i);
eb342d2d	347	Int_t ip=jparticle->GetPdgCode();
eb342d2d	348	if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
110410f9	349	if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	350	}
	351	} else {
	352	return kFALSE;
	353	}
	354
	355	return selected;
	356	}
	357
	358
eb342d2d	359	Bool_t AliGenHijing::SelectFlavor(Int_t pid)
	360	{
	361	// Select flavor of particle
	362	// 0: all
	363	// 4: charm and beauty
	364	// 5: beauty
	365	if (fFlavor == 0) return kTRUE;
	366
	367	Int_t ifl=TMath::Abs(pid/100);
	368	if (ifl > 10) ifl/=10;
	369	return ((fFlavor==4 && (ifl==4 \|\| ifl==5)) \|\|
	370	(fFlavor==5 && ifl==5));
	371
	372	}
110410f9	373
3b8bd63f	374	void AliGenHijing::MakeHeader()
	375	{
	376	// Builds the event header, to be called after each event
	377	AliGenHijingEventHeader* header = new AliGenHijingEventHeader("Hijing");
	378	// header->SetDate(date);
	379	// header->SetRunNumber(run);
	380	// header->SetEventNumber(event);
	381	header->SetNProduced(fHijing->GetNATT());
	382	header->SetImpactParameter(fHijing->GetHINT1(19));
	383	header->SetTotalEnergy(fHijing->GetEATT());
	384	header->SetHardScatters(fHijing->GetJATT());
	385	header->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
	386	header->SetCollisions(fHijing->GetN0(),
	387	fHijing->GetN01(),
	388	fHijing->GetN10(),
	389	fHijing->GetN11());
	390	}
	391
110410f9	392	AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
	393	{
	394	// Assignment operator
	395	return *this;
	396	}