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