[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.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 heavy mesons (J/Psi, upsilon or phi) in the
// the desired theta, phi and momentum windows; Gaussian smearing 
// on the vertex is done if selected


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

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