[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"

ClassImp(AliGenHerwig)


  AliGenHerwig::AliGenHerwig() :
    AliGenMC(),
    fAutPDF("GRV"),
    fModPDF(5),
    fStrucFunc(kGRVHO),
    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)
{
// Constructor
}

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

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


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

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);
  // compute parameter dependent constants
  fHerwig->PrepareRun();
}

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