[u/mrichter/AliRoot.git] / THerwig / AliGenHerwig.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 Herwig as an external generator
// The main Herwig options are accessable for the user through this interface.
// Uses the THerwig implementation of TGenerator.

#include "AliGenHerwig.h"
#include "AliHerwigRndm.h"
#include "AliRun.h"

#include <TParticle.h>
#include "THerwig6.h"

#include "Riostream.h"
#include "AliMC.h"

#include "driver.h"

ClassImp(AliGenHerwig)


  AliGenHerwig::AliGenHerwig() :
    AliGenMC(),
    fAutPDF("GRV"),
    fModPDF(5),
    fStrucFunc(kCTEQ5L),
    fKeep(0),
    fDecaysOff(1),
    fTrigger(0),
    fSelectAll(0),
    fFlavor(0),
    fEnergyCMS(14000),
    fMomentum1(7000),
    fMomentum2(7000),
    fKineBias(1),
    fTrials(0),
    fXsection(0),
    fHerwig(0x0),
    fProcess(0),
    fPtHardMin(0),
    fPtRMS(0),
    fMaxPr(10),
    fMaxErrors(1000),
    fEnSoft(1),
    fEv1Pr(0),
    fEv2Pr(0)
{
// Constructor
}

AliGenHerwig::AliGenHerwig(Int_t npart)
    :AliGenMC(npart)
{
    SetBeamMomenta();
    SetTarget();
    SetProjectile();
    SetStrucFunc(kCTEQ5L);
    fKeep=0;
    fTrigger=0;
    fDecaysOff=1;
    fSelectAll=0;
    fFlavor=0;
    fPtHardMin=0.;
    fPtRMS=0.0;
    fEnSoft=1.0;
    fMaxPr=10;
    fMaxErrors=1000;
    fEv1Pr=0;
    fEv2Pr=0;
    // Set random number generator   
    AliHerwigRndm::SetHerwigRandom(GetRandom());
}

AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
    :AliGenMC(Herwig)
{
// Copy constructor
    Herwig.Copy(*this);
}


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

void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
{
  fEv1Pr = ++eventFirst;
  fEv2Pr = ++eventLast;
  if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
}

void AliGenHerwig::Init()
{
// Initialisation
  fTarget.Resize(8);
  fProjectile.Resize(8);
  SetMC(new THerwig6());
  fHerwig=(THerwig6*) fMCEvGen;
  // initialize common blocks
  fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
  // reset parameters according to user needs
  InitPDF();
  fHerwig->SetPTMIN(fPtHardMin);
  fHerwig->SetPTRMS(fPtRMS);
  fHerwig->SetMAXPR(fMaxPr);
  fHerwig->SetMAXER(fMaxErrors);
  fHerwig->SetENSOF(fEnSoft);
  
  fHerwig->SetEV1PR(fEv1Pr);
  fHerwig->SetEV2PR(fEv2Pr);

// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
//       RMASS(1)=0.32
//       RMASS(2)=0.32
//       RMASS(3)=0.5
//       RMASS(4)=1.55
//       RMASS(5)=4.95
//       RMASS(6)=174.3
//       RMASS(13)=0.75

  fHerwig->SetRMASS(4,1.2);
  fHerwig->SetRMASS(5,4.75);
  
  if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);

  fHerwig->Hwusta("PI0     ");

  // compute parameter dependent constants
  fHerwig->PrepareRun();
}

void AliGenHerwig::InitPDF()
{
  switch(fStrucFunc)
    {
    case kGRVLO98:
      fModPDF=12;
      fAutPDF="GRV";
      break;
    case kCTEQ4M:
      fModPDF=34;
      fAutPDF="CTEQ";
      break;
    case kCTEQ5L:
      fModPDF=46;
      fAutPDF="CTEQ";
      break;
    case kCTEQ5M:
      fModPDF=48;
      fAutPDF="CTEQ";
      break;
    default:
      cerr << "This structure function is not inplemented " << fStrucFunc << endl;
      break;
    }
  fAutPDF.Resize(20);      
  fHerwig->SetMODPDF(1,fModPDF);
  fHerwig->SetMODPDF(2,fModPDF);
  fHerwig->SetAUTPDF(1,fAutPDF);
  fHerwig->SetAUTPDF(2,fAutPDF);
}

void AliGenHerwig::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[4], random[6];
  
  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) {
    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]));
    }
  }
  
  while(1)
    {
	fHerwig->GenerateEvent();
	fTrials++;
	fHerwig->ImportParticles(particles,"All");
	Int_t np = particles->GetEntriesFast()-1;
	if (np == 0 ) continue;
	
	Int_t nc=0;
	
	Int_t * newPos = new Int_t[np];
	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
	
	for (Int_t i = 0; i<np; i++) {
	    TParticle *  iparticle       = (TParticle *) particles->At(i);
	    imo = iparticle->GetFirstMother();
	    kf        = iparticle->GetPdgCode();
	    ks        = iparticle->GetStatusCode();
	    if (ks != 3 && 
		KinematicSelection(iparticle,0)) 
	    {
		nc++;
		p[0]=iparticle->Px();
		p[1]=iparticle->Py();
		p[2]=iparticle->Pz();
		p[3]=iparticle->Energy();
		origin[0]=origin0[0]+iparticle->Vx()/10;
		origin[1]=origin0[1]+iparticle->Vy()/10;
		origin[2]=origin0[2]+iparticle->Vz()/10;
		Float_t tof = kconv*iparticle->T();
		Int_t   iparent = (imo > -1) ? newPos[imo] : -1;
		Int_t   trackIt = (ks == 1) && fTrackIt;
		PushTrack(trackIt, iparent, kf,
			  p[0], p[1], p[2], p[3],
			  origin[0], origin[1], origin[2], 
			  tof,
			  polar[0], polar[1], polar[2],
			  kPPrimary, nt, fHerwig->GetEVWGT(), ks);
		KeepTrack(nt);
		newPos[i]=nt;
	    } // end of if: selection of particle
	} // end of for: particle loop
	if (newPos) delete[] newPos;
	//      MakeHeader();
	if (nc > 0) {
	    jev+=nc;
	    if (jev >= fNpart || fNpart == -1) {
		fKineBias=Float_t(fNpart)/Float_t(fTrials);
		break;
	    }
	}
    }
  SetHighWaterMark(nt);
//  adjust weight due to kinematic selection
  AdjustWeights();
//  get cross-section
  fXsection=fHerwig->GetAVWGT();
}

void AliGenHerwig::AdjustWeights()
{
// Adjust the weights after generation of all events
    TParticle *part;
    Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
    for (Int_t i=0; i<ntrack; i++) {
        part= gAlice->GetMCApp()->Particle(i);
        part->SetWeight(part->GetWeight()*fKineBias);
    }
}
 

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

Bool_t AliGenHerwig::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,0)&&SelectFlavor(ip)) {
		selected=kTRUE; break;
	    }
	    if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	}
    } else {
	return kFALSE;
    }

    return selected;
}


Bool_t AliGenHerwig::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 == ifl);
}

Bool_t AliGenHerwig::Stable(TParticle*  particle)
{
// Return true for a stable particle
//
    Int_t kf = TMath::Abs(particle->GetPdgCode());
    
    if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
	 
    {
	return kTRUE;
    } else {
	return kFALSE;
    }
}

void AliGenHerwig::FinishRun()
{
  fHerwig->Hwefin();
}


AliGenHerwig& AliGenHerwig::operator=(const  AliGenHerwig& rhs)
{
// Assignment operator
    rhs.Copy(*this);
    return (*this);
}
Commit	Line	Data
618e1dc0	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$ */
618e1dc0	17
	18
	19	// Generator using Herwig as an external generator
	20	// The main Herwig options are accessable for the user through this interface.
	21	// Uses the THerwig implementation of TGenerator.
	22
	23	#include "AliGenHerwig.h"
93831bcf	24	#include "AliHerwigRndm.h"
618e1dc0	25	#include "AliRun.h"
	26
	27	#include <TParticle.h>
	28	#include "THerwig6.h"
	29
	30	#include "Riostream.h"
5d12ce38	31	#include "AliMC.h"
618e1dc0	32
e2054d85	33	#include "driver.h"
e2054d85	34
014a9521	35	ClassImp(AliGenHerwig)
618e1dc0	36
7cdba479	37
a1d47d99	38	AliGenHerwig::AliGenHerwig() :
	39	AliGenMC(),
	40	fAutPDF("GRV"),
	41	fModPDF(5),
2ad1e1c2	42	fStrucFunc(kCTEQ5L),
a1d47d99	43	fKeep(0),
	44	fDecaysOff(1),
	45	fTrigger(0),
	46	fSelectAll(0),
	47	fFlavor(0),
	48	fEnergyCMS(14000),
	49	fMomentum1(7000),
	50	fMomentum2(7000),
	51	fKineBias(1),
	52	fTrials(0),
	53	fXsection(0),
	54	fHerwig(0x0),
	55	fProcess(0),
	56	fPtHardMin(0),
	57	fPtRMS(0),
	58	fMaxPr(10),
	59	fMaxErrors(1000),
e2054d85	60	fEnSoft(1),
	61	fEv1Pr(0),
	62	fEv2Pr(0)
618e1dc0	63	{
	64	// Constructor
	65	}
	66
	67	AliGenHerwig::AliGenHerwig(Int_t npart)
	68	:AliGenMC(npart)
	69	{
	70	SetBeamMomenta();
	71	SetTarget();
	72	SetProjectile();
2ad1e1c2	73	SetStrucFunc(kCTEQ5L);
618e1dc0	74	fKeep=0;
	75	fTrigger=0;
	76	fDecaysOff=1;
	77	fSelectAll=0;
	78	fFlavor=0;
93831bcf	79	fPtHardMin=0.;
618e1dc0	80	fPtRMS=0.0;
618e1dc0	81	fEnSoft=1.0;
93831bcf	82	fMaxPr=10;
618e1dc0	83	fMaxErrors=1000;
e2054d85	84	fEv1Pr=0;
e2054d85	85	fEv2Pr=0;
93831bcf	86	// Set random number generator
93831bcf	87	AliHerwigRndm::SetHerwigRandom(GetRandom());
618e1dc0	88	}
	89
	90	AliGenHerwig::AliGenHerwig(const AliGenHerwig & Herwig)
014a9521	91	:AliGenMC(Herwig)
618e1dc0	92	{
014a9521	93	// Copy constructor
014a9521	94	Herwig.Copy(*this);
618e1dc0	95	}
	96
	97
	98	AliGenHerwig::~AliGenHerwig()
	99	{
	100	// Destructor
	101	}
	102
e2054d85	103	void AliGenHerwig::SetEventListRange(Int_t eventFirst, Int_t eventLast)
	104	{
	105	fEv1Pr = ++eventFirst;
	106	fEv2Pr = ++eventLast;
	107	if ( fEv2Pr == -1 ) fEv2Pr = fEv2Pr;
	108	}
	109
618e1dc0	110	void AliGenHerwig::Init()
	111	{
	112	// Initialisation
	113	fTarget.Resize(8);
	114	fProjectile.Resize(8);
	115	SetMC(new THerwig6());
d274cd60	116	fHerwig=(THerwig6*) fMCEvGen;
618e1dc0	117	// initialize common blocks
	118	fHerwig->Initialize(fProjectile, fTarget, fMomentum1, fMomentum2, fProcess);
	119	// reset parameters according to user needs
	120	InitPDF();
	121	fHerwig->SetPTMIN(fPtHardMin);
	122	fHerwig->SetPTRMS(fPtRMS);
	123	fHerwig->SetMAXPR(fMaxPr);
	124	fHerwig->SetMAXER(fMaxErrors);
	125	fHerwig->SetENSOF(fEnSoft);
e2054d85	126
	127	fHerwig->SetEV1PR(fEv1Pr);
	128	fHerwig->SetEV2PR(fEv2Pr);
	129
	130	// C---D,U,S,C,B,T QUARK AND GLUON MASSES (IN THAT ORDER)
	131	// RMASS(1)=0.32
	132	// RMASS(2)=0.32
	133	// RMASS(3)=0.5
	134	// RMASS(4)=1.55
	135	// RMASS(5)=4.95
	136	// RMASS(6)=174.3
	137	// RMASS(13)=0.75
	138
	139	fHerwig->SetRMASS(4,1.2);
	140	fHerwig->SetRMASS(5,4.75);
	141
	142	if ( fProcess < 0 ) strncpy(VVJIN.QQIN,fFileName.Data(),50);
	143
	144	fHerwig->Hwusta("PI0 ");
	145
618e1dc0	146	// compute parameter dependent constants
	147	fHerwig->PrepareRun();
	148	}
	149
	150	void AliGenHerwig::InitPDF()
	151	{
	152	switch(fStrucFunc)
	153	{
618e1dc0	154	case kGRVLO98:
	155	fModPDF=12;
	156	fAutPDF="GRV";
	157	break;
618e1dc0	158	case kCTEQ4M:
	159	fModPDF=34;
	160	fAutPDF="CTEQ";
	161	break;
	162	case kCTEQ5L:
	163	fModPDF=46;
	164	fAutPDF="CTEQ";
	165	break;
e2054d85	166	case kCTEQ5M:
	167	fModPDF=48;
	168	fAutPDF="CTEQ";
	169	break;
618e1dc0	170	default:
	171	cerr << "This structure function is not inplemented " << fStrucFunc << endl;
	172	break;
	173	}
	174	fAutPDF.Resize(20);
	175	fHerwig->SetMODPDF(1,fModPDF);
	176	fHerwig->SetMODPDF(2,fModPDF);
	177	fHerwig->SetAUTPDF(1,fAutPDF);
	178	fHerwig->SetAUTPDF(2,fAutPDF);
	179	}
	180
	181	void AliGenHerwig::Generate()
	182	{
	183	// Generate one event
	184
	185	Float_t polar[3] = {0,0,0};
	186	Float_t origin[3]= {0,0,0};
	187	Float_t origin0[3]= {0,0,0};
	188	Float_t p[4], random[6];
	189
	190	static TClonesArray *particles;
	191	// converts from mm/c to s
	192	const Float_t kconv=0.001/2.999792458e8;
	193	//
	194	Int_t nt=0;
	195	Int_t jev=0;
	196	Int_t j, kf, ks, imo;
	197	kf=0;
	198
	199	if(!particles) particles=new TClonesArray("TParticle",10000);
	200
	201	fTrials=0;
	202	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
	203	if(fVertexSmear==kPerEvent) {
	204	Rndm(random,6);
	205	for (j=0;j<3;j++) {
	206	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	207	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	208	}
	209	}
	210
	211	while(1)
	212	{
	213	fHerwig->GenerateEvent();
	214	fTrials++;
	215	fHerwig->ImportParticles(particles,"All");
	216	Int_t np = particles->GetEntriesFast()-1;
	217	if (np == 0 ) continue;
	218
	219	Int_t nc=0;
	220
	221	Int_t * newPos = new Int_t[np];
	222	for (Int_t i = 0; i<np; i++) *(newPos+i)=-1;
	223
	224	for (Int_t i = 0; i<np; i++) {
	225	TParticle * iparticle = (TParticle *) particles->At(i);
	226	imo = iparticle->GetFirstMother();
	227	kf = iparticle->GetPdgCode();
	228	ks = iparticle->GetStatusCode();
	229	if (ks != 3 &&
	230	KinematicSelection(iparticle,0))
	231	{
	232	nc++;
	233	p[0]=iparticle->Px();
234	p[1]=iparticle->Py();
235	p[2]=iparticle->Pz();
236	p[3]=iparticle->Energy();
237	origin[0]=origin0[0]+iparticle->Vx()/10;
238	origin[1]=origin0[1]+iparticle->Vy()/10;
239	origin[2]=origin0[2]+iparticle->Vz()/10;
240	Float_t tof = kconv*iparticle->T();
241	Int_t iparent = (imo > -1) ? newPos[imo] : -1;
242	Int_t trackIt = (ks == 1) && fTrackIt;
93831bcf	243	PushTrack(trackIt, iparent, kf,
	244	p[0], p[1], p[2], p[3],
	245	origin[0], origin[1], origin[2],
	246	tof,
	247	polar[0], polar[1], polar[2],
e2054d85	248	kPPrimary, nt, fHerwig->GetEVWGT(), ks);
618e1dc0	249	KeepTrack(nt);
	250	newPos[i]=nt;
	251	} // end of if: selection of particle
	252	} // end of for: particle loop
	253	if (newPos) delete[] newPos;
618e1dc0	254	// MakeHeader();
618e1dc0	255	if (nc > 0) {
	256	jev+=nc;
	257	if (jev >= fNpart \|\| fNpart == -1) {
	258	fKineBias=Float_t(fNpart)/Float_t(fTrials);
618e1dc0	259	break;
	260	}
	261	}
	262	}
	263	SetHighWaterMark(nt);
	264	// adjust weight due to kinematic selection
	265	AdjustWeights();
	266	// get cross-section
	267	fXsection=fHerwig->GetAVWGT();
	268	}
	269
	270	void AliGenHerwig::AdjustWeights()
	271	{
	272	// Adjust the weights after generation of all events
	273	TParticle *part;
5d12ce38	274	Int_t ntrack=gAlice->GetMCApp()->GetNtrack();
618e1dc0	275	for (Int_t i=0; i<ntrack; i++) {
5d12ce38	276	part= gAlice->GetMCApp()->Particle(i);
618e1dc0	277	part->SetWeight(part->GetWeight()*fKineBias);
	278	}
	279	}
	280
	281
	282	void AliGenHerwig::KeepFullEvent()
	283	{
	284	fKeep=1;
	285	}
	286
	287	Bool_t AliGenHerwig::DaughtersSelection(TParticle* iparticle, TClonesArray* particles)
	288	{
	289	//
	290	// Looks recursively if one of the daughters has been selected
	291	//
	292	// printf("\n Consider daughters %d:",iparticle->GetPdgCode());
	293	Int_t imin=-1;
	294	Int_t imax=-1;
	295	Int_t i;
	296	Bool_t hasDaughters= (iparticle->GetFirstDaughter() >=0);
	297	Bool_t selected=kFALSE;
	298	if (hasDaughters) {
	299	imin=iparticle->GetFirstDaughter();
	300	imax=iparticle->GetLastDaughter();
	301	for (i=imin; i<= imax; i++){
	302	TParticle * jparticle = (TParticle *) particles->At(i);
	303	Int_t ip=jparticle->GetPdgCode();
	304	if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
	305	selected=kTRUE; break;
	306	}
	307	if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
	308	}
	309	} else {
	310	return kFALSE;
	311	}
	312
	313	return selected;
	314	}
	315
	316
	317	Bool_t AliGenHerwig::SelectFlavor(Int_t pid)
	318	{
	319	// Select flavor of particle
	320	// 0: all
	321	// 4: charm and beauty
	322	// 5: beauty
	323	if (fFlavor == 0) return kTRUE;
	324
	325	Int_t ifl=TMath::Abs(pid/100);
	326	if (ifl > 10) ifl/=10;
	327	return (fFlavor == ifl);
	328	}
	329
	330	Bool_t AliGenHerwig::Stable(TParticle* particle)
	331	{
	332	// Return true for a stable particle
	333	//
	334	Int_t kf = TMath::Abs(particle->GetPdgCode());
	335
	336	if ( (particle->GetFirstDaughter() < 0 ) \|\| (kf == 1000*fFlavor+122))
	337
	338	{
	339	return kTRUE;
	340	} else {
341	return kFALSE;
342	}
343	}
344
345	void AliGenHerwig::FinishRun()
346	{
347	fHerwig->Hwefin();
348	}
349
350
351	AliGenHerwig& AliGenHerwig::operator=(const AliGenHerwig& rhs)
352	{
353	// Assignment operator
014a9521	354	rhs.Copy(*this);
014a9521	355	return (*this);
618e1dc0	356	}
	357
	358