[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.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 "AliGenPHOSlib.h"
#include "AliRun.h"
#include "AliPythia.h"
#include <TDirectory.h>
#include <TFile.h>
#include <TTree.h>
#include <stdlib.h>
#include <TParticle.h>

ClassImp(AliGenParam)

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

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

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

//____________________________________________________________

AliGenParam::AliGenParam(Int_t npart, Param_t param) :AliGenerator(npart)
{
  //
  //  fName="HMESONpara";
  //  fTitle="Heavy Mesons Parametrisation";
  fPtParaFunc = AliGenMUONlib::GetPt(param);
  fYParaFunc  = AliGenMUONlib::GetY(param);
  fIpParaFunc = AliGenMUONlib::GetIp(param);
  
  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();  
}

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)
{
// 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();
}

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

//____________________________________________________________
void AliGenParam::Init()
{
    SetMC(new AliPythia());
    fPythia= (AliPythia*) fgMCEvGen;

//  End of the test !!!
  //Begin_Html
  /*
    <img src="picts/AliGenParam.gif">
  */
  //End_Html
 
  fPtPara = new TF1("Pt-Parametrization",fPtParaFunc,fPtMin,fPtMax,0);
  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==perEvent) {
      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==perTrack) {
	      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])+TMath::Pi();
			  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;

		  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)
{
    for (Int_t i=0; i<5; i++)
    {
	if (fChildSelect[i]==ip) return kTRUE;
    }
    return kFALSE;
}

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