[u/mrichter/AliRoot.git] / EVGEN / AliGenParam.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.18  2000/06/29 21:08:27  morsch
All paramatrisation libraries derive from the pure virtual base class AliGenLib.
This allows to pass a pointer to a library directly to AliGenParam and avoids the
use of function pointers in Config.C.

Revision 1.17  2000/06/09 20:33:30  morsch
All coding rule violations except RS3 corrected

Revision 1.16  2000/05/02 07:51:31  morsch
- Control precision of pT sampling TF1::SetNpx(..)
- Correct initialisation of child-cuts in all constructors.
- Most coding rule violations corrected.

Revision 1.15  2000/04/03 15:42:12  morsch
Cuts on primary particles are separated from those on the decay products. Methods
SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.

Revision 1.14  1999/11/09 07:38:48  fca
Changes for compatibility with version 2.23 of ROOT

Revision 1.13  1999/11/04 11:30:31  fca
Correct the logics for SetForceDecay

Revision 1.12  1999/11/03 17:43:20  fca
New version from G.Martinez & A.Morsch

Revision 1.11  1999/09/29 09:24:14  fca
Introduction of the Copyright and cvs Log

*/

#include "AliGenParam.h"
#include "AliGenMUONlib.h"
#include "AliRun.h"
#include "AliPythia.h"
#include <TParticle.h>
#include <TF1.h>

ClassImp(AliGenParam)

//------------------------------------------------------------

  //Begin_Html
  /*
    <img src="picts/AliGenParam.gif">
  */
  //End_Html

//____________________________________________________________
//____________________________________________________________
AliGenParam::AliGenParam()
                 :AliGenerator()
{
// Deafault constructor
    fPtPara = 0;
    fYPara  = 0;
    fParam  = jpsi_p;
    fAnalog = analog;
    SetCutOnChild();
    SetChildMomentumRange();
    SetChildPtRange();
    SetChildPhiRange();
    SetChildThetaRange();  
    SetDeltaPt();
}

AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library,  Param_t param, char* tname):AliGenerator(npart)
{
// Constructor using number of particles parameterisation id and library
    
    fPtParaFunc = Library->GetPt(param, tname);
    fYParaFunc  = Library->GetY (param, tname);
    fIpParaFunc = Library->GetIp(param, tname);
    
    fPtPara = 0;
    fYPara  = 0;
    fParam  = param;
    fAnalog = analog;
    fChildSelect.Set(5);
    for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
    SetForceDecay();
    SetCutOnChild();
    SetChildMomentumRange();
    SetChildPtRange();
    SetChildPhiRange();
    SetChildThetaRange(); 
    SetDeltaPt(); 
}

//____________________________________________________________

AliGenParam::AliGenParam(Int_t npart, Param_t param, char* tname):AliGenerator(npart)
{
// Constructor using parameterisation id and number of particles
//      
    AliGenLib* Library = new AliGenMUONlib();
 
    fPtParaFunc = Library->GetPt(param, tname);
    fYParaFunc  = Library->GetY (param, tname);
    fIpParaFunc = Library->GetIp(param, tname);
    
    fPtPara = 0;
    fYPara  = 0;
    fParam  = param;
    fAnalog = analog;
    fChildSelect.Set(5);
    for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
    SetForceDecay();
    SetCutOnChild();
    SetChildMomentumRange();
    SetChildPtRange();
    SetChildPhiRange();
    SetChildThetaRange(); 
    SetDeltaPt(); 
}

AliGenParam::AliGenParam(Int_t npart, Param_t param,
                         Double_t (*PtPara) (Double_t*, Double_t*),
                         Double_t (*YPara ) (Double_t* ,Double_t*),
		         Int_t    (*IpPara) ())			 
    :AliGenerator(npart)
{
// Constructor
// Gines Martinez 1/10/99 
    fPtParaFunc = PtPara; 
    fYParaFunc  = YPara;  
    fIpParaFunc = IpPara;
//  
    fPtPara = 0;
    fYPara  = 0;
    fParam  = param;
    fAnalog = analog;
    fChildSelect.Set(5);
    for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
    SetForceDecay();
    SetCutOnChild();
    SetChildMomentumRange();
    SetChildPtRange();
    SetChildPhiRange();
    SetChildThetaRange();  
    SetDeltaPt();
}


AliGenParam::AliGenParam(const AliGenParam & Paramd)
{
// copy constructor
}

//____________________________________________________________
AliGenParam::~AliGenParam()
{
// Destructor
    delete  fPtPara;
    delete  fYPara;
}

//____________________________________________________________
void AliGenParam::Init()
{
// Initialisation
    SetMC(new AliPythia());
    fPythia= (AliPythia*) fgMCEvGen;
    
  //Begin_Html
  /*
    <img src="picts/AliGenParam.gif">
  */
  //End_Html
 
    fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
//  Set representation precision to 10 MeV
    Int_t npx= Int_t((fPtMax-fPtMin)/fDeltaPt);
    
    fPtPara->SetNpx(npx);
    
    fYPara  = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
    TF1* ptPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
    TF1* yPara  = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);

//
// dN/dy| y=0
    Double_t y1=0;
    Double_t y2=0;
    
    fdNdy0=fYParaFunc(&y1,&y2);
//
// Integral over generation region
    Float_t intYS  = yPara ->Integral(fYMin, fYMax);
    Float_t intPt0 = ptPara->Integral(0,15);
    Float_t intPtS = ptPara->Integral(fPtMin,fPtMax);
    Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
    if (fAnalog == analog) {
	fYWgt  = intYS/fdNdy0;
	fPtWgt = intPtS/intPt0;
	fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
    } else {
	fYWgt = intYS/fdNdy0;
	fPtWgt = (fPtMax-fPtMin)/intPt0;
	fParentWeight = fYWgt*fPtWgt*phiWgt/fNpart;
    }
//
// particle decay related initialization
    fPythia->DefineParticles();
// semimuonic decays of charm and beauty
    fPythia->ForceDecay(fForceDecay);
//
    switch (fForceDecay) 
    {
    case semielectronic:
    case dielectron:
    case b_jpsi_dielectron:
    case b_psip_dielectron:
	fChildSelect[0]=11;	
	break;
    case semimuonic:
    case dimuon:
    case b_jpsi_dimuon:
    case b_psip_dimuon:
	fChildSelect[0]=13;
	break;
    case pitomu:
	fChildSelect[0]=13;
	break;
    case katomu:
	fChildSelect[0]=13;
	break;
    case nodecay:
	break;
    case all:
	break;
    }
}

//____________________________________________________________
void AliGenParam::Generate()
{
//
// Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
// Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and 
// antineutrons in the the desired theta, phi and momentum windows; 
// Gaussian smearing on the vertex is done if selected. 
// The decay of heavy mesons is done using lujet, 
//    and the childern particle are tracked by GEANT
// However, light mesons are directly tracked by GEANT 
// setting fForceDecay = nodecay (SetForceDecay(nodecay)) 
//


  Float_t polar[3]= {0,0,0};  // Polarisation of the parent particle (for GEANT tracking)
  Float_t origin0[3];         // Origin of the generated parent particle (for GEANT tracking)
  Float_t pt, pl, ptot;       // Transverse, logitudinal and total momenta of the parent particle
  Float_t phi, theta;         // Phi and theta spherical angles of the parent particle momentum
  Float_t p[3], pc[3], 
          och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
  Float_t ty, xmt;
  Int_t nt, i, j, kfch[10];
  Float_t  wgtp, wgtch;
  Double_t dummy;
  static TClonesArray *particles;
  //
  if(!particles) particles=new TClonesArray("TParticle",1000);
  //
  Float_t random[6];
 
// Calculating vertex position per event
  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]));
      }
  }
  Int_t ipa=0;
// Generating fNpart particles
  while (ipa<fNpart) {
      while(1) {
//
// particle type
	  Int_t iPart = fIpParaFunc();
	  fChildWeight=(fPythia->GetBraPart(iPart))*fParentWeight;	  
	  Float_t am=fPythia->GetPMAS(fPythia->Lucomp(iPart),1);
	  gMC->Rndm(random,2);
//
// phi
	  phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
//
// y
	  ty=Float_t(TMath::TanH(fYPara->GetRandom()));
//
// pT
	  if (fAnalog == analog) {
	      pt=fPtPara->GetRandom();
	      wgtp=fParentWeight;
	      wgtch=fChildWeight;
	  } else {
	      pt=fPtMin+random[1]*(fPtMax-fPtMin);
	      Double_t ptd=pt;
	      wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
	      wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
	  }
	  xmt=sqrt(pt*pt+am*am);
	  pl=xmt*ty/sqrt(1.-ty*ty);
	  theta=TMath::ATan2(pt,pl);
// Cut on theta
	  if(theta<fThetaMin || theta>fThetaMax) continue;
	  ptot=TMath::Sqrt(pt*pt+pl*pl);
// Cut on momentum
	  if(ptot<fPMin || ptot>fPMax) continue;
	  p[0]=pt*TMath::Cos(phi);
	  p[1]=pt*TMath::Sin(phi);
	  p[2]=pl;
	  if(fVertexSmear==kPerTrack) {
	      gMC->Rndm(random,6);
	      for (j=0;j<3;j++) {
		  origin0[j]=
		      fOrigin[j]+fOsigma[j]*TMath::Cos(2*random[2*j]*TMath::Pi())*
		      TMath::Sqrt(-2*TMath::Log(random[2*j+1]));
	      }
	  }
	  
// Looking at fForceDecay : 
// if fForceDecay != none Primary particle decays using 
// AliPythia and children are tracked by GEANT
//
// if fForceDecay == none Primary particle is tracked by GEANT 
// (In the latest, make sure that GEANT actually does all the decays you want)	  
//
	  if (fForceDecay != nodecay) {
// Using lujet to decay particle
	      Float_t energy=TMath::Sqrt(ptot*ptot+am*am);
	      fPythia->DecayParticle(iPart,energy,theta,phi);
//
// select decay particles
	      Int_t np=fPythia->ImportParticles(particles,"All");
	      Int_t ncsel=0;
	      for (i = 1; i<np; i++) {
		  TParticle *  iparticle = (TParticle *) particles->At(i);
		  Int_t kf = iparticle->GetPdgCode();
//
// children
		  if (ChildSelected(TMath::Abs(kf)))
		  {
		      pc[0]=iparticle->Px();
		      pc[1]=iparticle->Py();
		      pc[2]=iparticle->Pz();
		      och[0]=origin0[0]+iparticle->Vx()/10;
		      och[1]=origin0[1]+iparticle->Vy()/10;
		      och[2]=origin0[2]+iparticle->Vz()/10;
		      if (fCutOnChild) {
			  Float_t ptChild=TMath::Sqrt(pc[0]*pc[0]+pc[1]*pc[1]);
			  Float_t pChild=TMath::Sqrt(ptChild*ptChild+pc[2]*pc[2]);
			  Float_t thetaChild=TMath::ATan2(ptChild,pc[2]);
			  Float_t phiChild=TMath::ATan2(pc[1],pc[0]);
			  Bool_t childok = 
			      ((ptChild    > fChildPtMin    && ptChild    <fChildPtMax)      &&
			       (pChild     > fChildPMin     && pChild     <fChildPMax)       &&
			       (thetaChild > fChildThetaMin && thetaChild <fChildThetaMax)   &&
			       (phiChild   > fChildPhiMin   && phiChild   <fChildPhiMax));
			  if(childok)
			  {
			      pch[ncsel][0]=pc[0];
			      pch[ncsel][1]=pc[1];
			      pch[ncsel][2]=pc[2];
			      kfch[ncsel]=kf;
			      ncsel++;
			  } else {
			      ncsel=-1;
			      break;
			  } // child kine cuts
		      } else {
			  pch[ncsel][0]=pc[0];
			  pch[ncsel][1]=pc[1];
			  pch[ncsel][2]=pc[2];
			  kfch[ncsel]=kf;
			  ncsel++;
		      } // if child selection
		  } // select muon
	      } // decay particle loop
	      Int_t iparent;
	      if ((fCutOnChild && ncsel >0) || !fCutOnChild){
		  ipa++;
//
// parent
		  gAlice->
		      SetTrack(0,-1,iPart,p,origin0,polar,0,"Primary",nt,wgtp);
		  iparent=nt;
		  gAlice->KeepTrack(nt); 
		  for (i=0; i< ncsel; i++) {
		      gAlice->SetTrack(fTrackIt,iparent,kfch[i],
				       &pch[i][0],och,polar,
				       0,"Decay",nt,wgtch);
		      gAlice->KeepTrack(nt); 
		  }
	      }  // Decays by Lujet
	  } // kinematic selection
	  else  // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
	  {
	    gAlice->
		SetTrack(fTrackIt,-1,iPart,p,origin0,polar,0,"Primary",nt,wgtp);
            ipa++; 
	  }
	  break;
    } // while
  } // event loop
}

Bool_t AliGenParam::ChildSelected(Int_t ip)
{
// True if particle is in list of selected children
    for (Int_t i=0; i<5; i++)
    {
	if (fChildSelect[i]==ip) return kTRUE;
    }
    return kFALSE;
}

Bool_t AliGenParam::KinematicSelection(TParticle *particle)
{
// Perform kinematic cuts
    Float_t px=particle->Px();
    Float_t py=particle->Py();
    Float_t pz=particle->Pz();
//
// 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;
    }
    Float_t pt=TMath::Sqrt(px*px+py*py);
    
//
// theta cut
    Float_t  theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
    if (theta > fThetaMax || theta < fThetaMin) 
    {
//	printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
	return kFALSE;
    }

    return kTRUE;
}


AliGenParam& AliGenParam::operator=(const  AliGenParam& rhs)
{
// Assignment operator
    return *this;
}
Commit	Line	Data
4c039060	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$
aee8290b	18	Revision 1.18 2000/06/29 21:08:27 morsch
	19	All paramatrisation libraries derive from the pure virtual base class AliGenLib.
	20	This allows to pass a pointer to a library directly to AliGenParam and avoids the
	21	use of function pointers in Config.C.
	22
b22ee262	23	Revision 1.17 2000/06/09 20:33:30 morsch
	24	All coding rule violations except RS3 corrected
	25
f87cfe57	26	Revision 1.16 2000/05/02 07:51:31 morsch
	27	- Control precision of pT sampling TF1::SetNpx(..)
	28	- Correct initialisation of child-cuts in all constructors.
	29	- Most coding rule violations corrected.
	30
749070b6	31	Revision 1.15 2000/04/03 15:42:12 morsch
	32	Cuts on primary particles are separated from those on the decay products. Methods
	33	SetChildMomentumRange, SetChildPtRange, SetChildPhiRange, SetChildThetaRange added.
	34
cc5d764c	35	Revision 1.14 1999/11/09 07:38:48 fca
	36	Changes for compatibility with version 2.23 of ROOT
	37
084c1b4a	38	Revision 1.13 1999/11/04 11:30:31 fca
	39	Correct the logics for SetForceDecay
	40
28337bc1	41	Revision 1.12 1999/11/03 17:43:20 fca
	42	New version from G.Martinez & A.Morsch
	43
886b6f73	44	Revision 1.11 1999/09/29 09:24:14 fca
	45	Introduction of the Copyright and cvs Log
	46
4c039060	47	*/
4c039060	48
fe4da5cc	49	#include "AliGenParam.h"
	50	#include "AliGenMUONlib.h"
	51	#include "AliRun.h"
fe4da5cc	52	#include "AliPythia.h"
1578254f	53	#include <TParticle.h>
f87cfe57	54	#include <TF1.h>
fe4da5cc	55
	56	ClassImp(AliGenParam)
	57
	58	//------------------------------------------------------------
	59
	60	//Begin_Html
	61	/*
1439f98e	62	<img src="picts/AliGenParam.gif">
fe4da5cc	63	*/
	64	//End_Html
	65
	66	//____________________________________________________________
	67	//____________________________________________________________
	68	AliGenParam::AliGenParam()
	69	:AliGenerator()
	70	{
b22ee262	71	// Deafault constructor
	72	fPtPara = 0;
	73	fYPara = 0;
	74	fParam = jpsi_p;
	75	fAnalog = analog;
	76	SetCutOnChild();
	77	SetChildMomentumRange();
	78	SetChildPtRange();
	79	SetChildPhiRange();
	80	SetChildThetaRange();
	81	SetDeltaPt();
	82	}
	83
	84	AliGenParam::AliGenParam(Int_t npart, AliGenLib * Library, Param_t param, char* tname):AliGenerator(npart)
	85	{
	86	// Constructor using number of particles parameterisation id and library
	87
	88	fPtParaFunc = Library->GetPt(param, tname);
	89	fYParaFunc = Library->GetY (param, tname);
	90	fIpParaFunc = Library->GetIp(param, tname);
	91
	92	fPtPara = 0;
	93	fYPara = 0;
	94	fParam = param;
	95	fAnalog = analog;
	96	fChildSelect.Set(5);
	97	for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
	98	SetForceDecay();
	99	SetCutOnChild();
	100	SetChildMomentumRange();
	101	SetChildPtRange();
	102	SetChildPhiRange();
	103	SetChildThetaRange();
	104	SetDeltaPt();
fe4da5cc	105	}
	106
	107	//____________________________________________________________
886b6f73	108
b22ee262	109	AliGenParam::AliGenParam(Int_t npart, Param_t param, char* tname):AliGenerator(npart)
fe4da5cc	110	{
b22ee262	111	// Constructor using parameterisation id and number of particles
	112	//
	113	AliGenLib* Library = new AliGenMUONlib();
	114
	115	fPtParaFunc = Library->GetPt(param, tname);
	116	fYParaFunc = Library->GetY (param, tname);
	117	fIpParaFunc = Library->GetIp(param, tname);
	118
	119	fPtPara = 0;
	120	fYPara = 0;
	121	fParam = param;
	122	fAnalog = analog;
	123	fChildSelect.Set(5);
	124	for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
	125	SetForceDecay();
	126	SetCutOnChild();
	127	SetChildMomentumRange();
	128	SetChildPtRange();
	129	SetChildPhiRange();
	130	SetChildThetaRange();
	131	SetDeltaPt();
fe4da5cc	132	}
fe4da5cc	133
886b6f73	134	AliGenParam::AliGenParam(Int_t npart, Param_t param,
	135	Double_t (PtPara) (Double_t, Double_t*),
	136	Double_t (YPara ) (Double_t ,Double_t*),
	137	Int_t (*IpPara) ())
	138	:AliGenerator(npart)
	139	{
749070b6	140	// Constructor
28337bc1	141	// Gines Martinez 1/10/99
886b6f73	142	fPtParaFunc = PtPara;
	143	fYParaFunc = YPara;
	144	fIpParaFunc = IpPara;
28337bc1	145	//
886b6f73	146	fPtPara = 0;
	147	fYPara = 0;
	148	fParam = param;
	149	fAnalog = analog;
	150	fChildSelect.Set(5);
	151	for (Int_t i=0; i<5; i++) fChildSelect[i]=0;
	152	SetForceDecay();
	153	SetCutOnChild();
749070b6	154	SetChildMomentumRange();
	155	SetChildPtRange();
	156	SetChildPhiRange();
	157	SetChildThetaRange();
	158	SetDeltaPt();
886b6f73	159	}
886b6f73	160
f87cfe57	161
	162	AliGenParam::AliGenParam(const AliGenParam & Paramd)
	163	{
	164	// copy constructor
	165	}
	166
fe4da5cc	167	//____________________________________________________________
	168	AliGenParam::~AliGenParam()
	169	{
749070b6	170	// Destructor
fe4da5cc	171	delete fPtPara;
	172	delete fYPara;
	173	}
	174
	175	//____________________________________________________________
	176	void AliGenParam::Init()
	177	{
749070b6	178	// Initialisation
fe4da5cc	179	SetMC(new AliPythia());
fe4da5cc	180	fPythia= (AliPythia*) fgMCEvGen;
749070b6	181
fe4da5cc	182	//Begin_Html
fe4da5cc	183	/*
1439f98e	184	<img src="picts/AliGenParam.gif">
fe4da5cc	185	*/
	186	//End_Html
	187
749070b6	188	fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
	189	// Set representation precision to 10 MeV
	190	Int_t npx= Int_t((fPtMax-fPtMin)/fDeltaPt);
	191
	192	fPtPara->SetNpx(npx);
	193
	194	fYPara = new TF1("Y -Parametrization",fYParaFunc,fYMin,fYMax,0);
	195	TF1* ptPara = new TF1("Pt-Parametrization",fPtParaFunc,0,15,0);
	196	TF1* yPara = new TF1("Y -Parametrization",fYParaFunc,-6,6,0);
b7601ac4	197
fe4da5cc	198	//
fe4da5cc	199	// dN/dy\| y=0
749070b6	200	Double_t y1=0;
	201	Double_t y2=0;
	202
	203	fdNdy0=fYParaFunc(&y1,&y2);
fe4da5cc	204	//
fe4da5cc	205	// Integral over generation region
749070b6	206	Float_t intYS = yPara ->Integral(fYMin, fYMax);
	207	Float_t intPt0 = ptPara->Integral(0,15);
	208	Float_t intPtS = ptPara->Integral(fPtMin,fPtMax);
	209	Float_t phiWgt=(fPhiMax-fPhiMin)/2./TMath::Pi();
	210	if (fAnalog == analog) {
	211	fYWgt = intYS/fdNdy0;
	212	fPtWgt = intPtS/intPt0;
	213	fParentWeight = fYWgtfPtWgtphiWgt/fNpart;
	214	} else {
	215	fYWgt = intYS/fdNdy0;
	216	fPtWgt = (fPtMax-fPtMin)/intPt0;
	217	fParentWeight = fYWgtfPtWgtphiWgt/fNpart;
	218	}
fe4da5cc	219	//
fe4da5cc	220	// particle decay related initialization
749070b6	221	fPythia->DefineParticles();
fe4da5cc	222	// semimuonic decays of charm and beauty
749070b6	223	fPythia->ForceDecay(fForceDecay);
fe4da5cc	224	//
	225	switch (fForceDecay)
	226	{
	227	case semielectronic:
	228	case dielectron:
	229	case b_jpsi_dielectron:
	230	case b_psip_dielectron:
	231	fChildSelect[0]=11;
	232	break;
	233	case semimuonic:
	234	case dimuon:
	235	case b_jpsi_dimuon:
	236	case b_psip_dimuon:
	237	fChildSelect[0]=13;
	238	break;
b7601ac4	239	case pitomu:
	240	fChildSelect[0]=13;
	241	break;
	242	case katomu:
	243	fChildSelect[0]=13;
	244	break;
886b6f73	245	case nodecay:
	246	break;
	247	case all:
	248	break;
fe4da5cc	249	}
fe4da5cc	250	}
	251
	252	//____________________________________________________________
	253	void AliGenParam::Generate()
	254	{
28337bc1	255	//
	256	// Generate 'npart' of light and heavy mesons (J/Psi, upsilon or phi, Pion,
	257	// Kaons, Etas, Omegas) and Baryons (proton, antiprotons, neutrons and
	258	// antineutrons in the the desired theta, phi and momentum windows;
	259	// Gaussian smearing on the vertex is done if selected.
cc5d764c	260	// The decay of heavy mesons is done using lujet,
	261	// and the childern particle are tracked by GEANT
	262	// However, light mesons are directly tracked by GEANT
	263	// setting fForceDecay = nodecay (SetForceDecay(nodecay))
28337bc1	264	//
fe4da5cc	265
21aaa175	266
28337bc1	267	Float_t polar[3]= {0,0,0}; // Polarisation of the parent particle (for GEANT tracking)
	268	Float_t origin0[3]; // Origin of the generated parent particle (for GEANT tracking)
	269	Float_t pt, pl, ptot; // Transverse, logitudinal and total momenta of the parent particle
	270	Float_t phi, theta; // Phi and theta spherical angles of the parent particle momentum
	271	Float_t p[3], pc[3],
	272	och[3], pch[10][3]; // Momentum, polarisation and origin of the children particles from lujet
fe4da5cc	273	Float_t ty, xmt;
21aaa175	274	Int_t nt, i, j, kfch[10];
fe4da5cc	275	Float_t wgtp, wgtch;
fe4da5cc	276	Double_t dummy;
09fd3ea2	277	static TClonesArray *particles;
fe4da5cc	278	//
09fd3ea2	279	if(!particles) particles=new TClonesArray("TParticle",1000);
fe4da5cc	280	//
fe4da5cc	281	Float_t random[6];
28337bc1	282
28337bc1	283	// Calculating vertex position per event
fe4da5cc	284	for (j=0;j<3;j++) origin0[j]=fOrigin[j];
aee8290b	285	if(fVertexSmear==kPerEvent) {
cfce8870	286	gMC->Rndm(random,6);
fe4da5cc	287	for (j=0;j<3;j++) {
	288	origin0[j]+=fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	289	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	290	}
	291	}
21aaa175	292	Int_t ipa=0;
28337bc1	293	// Generating fNpart particles
21aaa175	294	while (ipa<fNpart) {
cc5d764c	295	while(1) {
fe4da5cc	296	//
fe4da5cc	297	// particle type
749070b6	298	Int_t iPart = fIpParaFunc();
	299	fChildWeight=(fPythia->GetBraPart(iPart))*fParentWeight;
	300	Float_t am=fPythia->GetPMAS(fPythia->Lucomp(iPart),1);
cfce8870	301	gMC->Rndm(random,2);
fe4da5cc	302	//
	303	// phi
	304	phi=fPhiMin+random[0]*(fPhiMax-fPhiMin);
	305	//
	306	// y
	307	ty=Float_t(TMath::TanH(fYPara->GetRandom()));
	308	//
	309	// pT
82db6e43	310	if (fAnalog == analog) {
fe4da5cc	311	pt=fPtPara->GetRandom();
	312	wgtp=fParentWeight;
	313	wgtch=fChildWeight;
	314	} else {
	315	pt=fPtMin+random[1]*(fPtMax-fPtMin);
	316	Double_t ptd=pt;
	317	wgtp=fParentWeight*fPtParaFunc(& ptd, &dummy);
	318	wgtch=fChildWeight*fPtParaFunc(& ptd, &dummy);
	319	}
	320	xmt=sqrt(ptpt+amam);
	321	pl=xmtty/sqrt(1.-tyty);
	322	theta=TMath::ATan2(pt,pl);
cc5d764c	323	// Cut on theta
fe4da5cc	324	if(theta<fThetaMin \|\| theta>fThetaMax) continue;
fe4da5cc	325	ptot=TMath::Sqrt(ptpt+plpl);
cc5d764c	326	// Cut on momentum
fe4da5cc	327	if(ptot<fPMin \|\| ptot>fPMax) continue;
	328	p[0]=pt*TMath::Cos(phi);
	329	p[1]=pt*TMath::Sin(phi);
	330	p[2]=pl;
aee8290b	331	if(fVertexSmear==kPerTrack) {
cfce8870	332	gMC->Rndm(random,6);
fe4da5cc	333	for (j=0;j<3;j++) {
	334	origin0[j]=
	335	fOrigin[j]+fOsigma[j]TMath::Cos(2random[2j]TMath::Pi())*
	336	TMath::Sqrt(-2TMath::Log(random[2j+1]));
	337	}
	338	}
28337bc1	339
28337bc1	340	// Looking at fForceDecay :
cc5d764c	341	// if fForceDecay != none Primary particle decays using
cc5d764c	342	// AliPythia and children are tracked by GEANT
28337bc1	343	//
cc5d764c	344	// if fForceDecay == none Primary particle is tracked by GEANT
cc5d764c	345	// (In the latest, make sure that GEANT actually does all the decays you want)
fe4da5cc	346	//
886b6f73	347	if (fForceDecay != nodecay) {
28337bc1	348	// Using lujet to decay particle
886b6f73	349	Float_t energy=TMath::Sqrt(ptotptot+amam);
749070b6	350	fPythia->DecayParticle(iPart,energy,theta,phi);
fe4da5cc	351	//
cc5d764c	352	// select decay particles
886b6f73	353	Int_t np=fPythia->ImportParticles(particles,"All");
886b6f73	354	Int_t ncsel=0;
	355	for (i = 1; i<np; i++) {
	356	TParticle * iparticle = (TParticle *) particles->At(i);
	357	Int_t kf = iparticle->GetPdgCode();
fe4da5cc	358	//
fe4da5cc	359	// children
886b6f73	360	if (ChildSelected(TMath::Abs(kf)))
	361	{
	362	pc[0]=iparticle->Px();
	363	pc[1]=iparticle->Py();
	364	pc[2]=iparticle->Pz();
	365	och[0]=origin0[0]+iparticle->Vx()/10;
	366	och[1]=origin0[1]+iparticle->Vy()/10;
	367	och[2]=origin0[2]+iparticle->Vz()/10;
	368	if (fCutOnChild) {
749070b6	369	Float_t ptChild=TMath::Sqrt(pc[0]pc[0]+pc[1]pc[1]);
	370	Float_t pChild=TMath::Sqrt(ptChildptChild+pc[2]pc[2]);
	371	Float_t thetaChild=TMath::ATan2(ptChild,pc[2]);
	372	Float_t phiChild=TMath::ATan2(pc[1],pc[0]);
886b6f73	373	Bool_t childok =
749070b6	374	((ptChild > fChildPtMin && ptChild <fChildPtMax) &&
	375	(pChild > fChildPMin && pChild <fChildPMax) &&
	376	(thetaChild > fChildThetaMin && thetaChild <fChildThetaMax) &&
	377	(phiChild > fChildPhiMin && phiChild <fChildPhiMax));
886b6f73	378	if(childok)
	379	{
	380	pch[ncsel][0]=pc[0];
	381	pch[ncsel][1]=pc[1];
	382	pch[ncsel][2]=pc[2];
	383	kfch[ncsel]=kf;
	384	ncsel++;
	385	} else {
	386	ncsel=-1;
	387	break;
	388	} // child kine cuts
21aaa175	389	} else {
886b6f73	390	pch[ncsel][0]=pc[0];
	391	pch[ncsel][1]=pc[1];
	392	pch[ncsel][2]=pc[2];
	393	kfch[ncsel]=kf;
	394	ncsel++;
	395	} // if child selection
	396	} // select muon
	397	} // decay particle loop
	398	Int_t iparent;
	399	if ((fCutOnChild && ncsel >0) \|\| !fCutOnChild){
	400	ipa++;
21aaa175	401	//
21aaa175	402	// parent
886b6f73	403	gAlice->
749070b6	404	SetTrack(0,-1,iPart,p,origin0,polar,0,"Primary",nt,wgtp);
886b6f73	405	iparent=nt;
f87cfe57	406	gAlice->KeepTrack(nt);
886b6f73	407	for (i=0; i< ncsel; i++) {
	408	gAlice->SetTrack(fTrackIt,iparent,kfch[i],
	409	&pch[i][0],och,polar,
	410	0,"Decay",nt,wgtch);
	411	gAlice->KeepTrack(nt);
	412	}
28337bc1	413	} // Decays by Lujet
21aaa175	414	} // kinematic selection
28337bc1	415	else // nodecay option, so parent will be tracked by GEANT (pions, kaons, eta, omegas, baryons)
	416	{
	417	gAlice->
749070b6	418	SetTrack(fTrackIt,-1,iPart,p,origin0,polar,0,"Primary",nt,wgtp);
28337bc1	419	ipa++;
28337bc1	420	}
fe4da5cc	421	break;
21aaa175	422	} // while
fe4da5cc	423	} // event loop
	424	}
	425
	426	Bool_t AliGenParam::ChildSelected(Int_t ip)
	427	{
749070b6	428	// True if particle is in list of selected children
fe4da5cc	429	for (Int_t i=0; i<5; i++)
	430	{
	431	if (fChildSelect[i]==ip) return kTRUE;
	432	}
	433	return kFALSE;
	434	}
	435
1578254f	436	Bool_t AliGenParam::KinematicSelection(TParticle *particle)
fe4da5cc	437	{
749070b6	438	// Perform kinematic cuts
1578254f	439	Float_t px=particle->Px();
	440	Float_t py=particle->Py();
	441	Float_t pz=particle->Pz();
fe4da5cc	442	//
	443	// momentum cut
	444	Float_t p=TMath::Sqrt(pxpx+pypy+pz*pz);
	445	if (p > fPMax \|\| p < fPMin)
	446	{
	447	// printf("\n failed p cut %f %f %f \n",p,fPMin,fPMax);
	448	return kFALSE;
	449	}
	450	Float_t pt=TMath::Sqrt(pxpx+pypy);
	451
	452	//
	453	// theta cut
	454	Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(p)));
	455	if (theta > fThetaMax \|\| theta < fThetaMin)
	456	{
	457	// printf("\n failed theta cut %f %f %f \n",theta,fThetaMin,fThetaMax);
	458	return kFALSE;
	459	}
	460
	461	return kTRUE;
	462	}
	463
	464
f87cfe57	465	AliGenParam& AliGenParam::operator=(const AliGenParam& rhs)
	466	{
	467	// Assignment operator
	468	return *this;
	469	}
	470
fe4da5cc	471
fe4da5cc	472