[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.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();
}

//____________________________________________________________

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

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)) 
//

// printf("Generate !!!!!!!!!!!!!\n");

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