/*
$Log$
+Revision 1.11 2000/10/20 13:22:26 morsch
+- skip particle type 92 (string)
+- Charmed and beauty baryions (5122, 4122) are considered as stable consistent with
+ mesons.
+
+Revision 1.10 2000/10/17 15:10:20 morsch
+Write first all the parent particles to the stack and then the final state particles.
+
+Revision 1.9 2000/10/17 13:38:59 morsch
+Protection against division by zero in EvaluateCrossSection() and KinematicSelection(..) (FCA)
+
+Revision 1.8 2000/10/17 12:46:31 morsch
+Protect EvaluateCrossSections() against division by zero.
+
+Revision 1.7 2000/10/02 21:28:06 fca
+Removal of useless dependecies via forward declarations
+
+Revision 1.6 2000/09/11 13:23:37 morsch
+Write last seed to file (fortran lun 50) and reed back from same lun using calls to
+luget_hijing and luset_hijing.
+
+Revision 1.5 2000/09/07 16:55:40 morsch
+fHijing->Initialize(); after change of parameters. (Dmitri Yurevitch Peressounko)
+
+Revision 1.4 2000/07/11 18:24:56 fca
+Coding convention corrections + few minor bug fixes
+
Revision 1.3 2000/06/30 12:08:36 morsch
In member data: char* replaced by TString, Init takes care of resizing the strings to
8 characters required by Hijing.
#include "AliGenHijing.h"
#include "AliGenHijingEventHeader.h"
#include "AliRun.h"
+#include "AliMC.h"
#include <TArrayI.h>
#include <TParticle.h>
fMinImpactParam, fMaxImpactParam));
fHijing=(THijing*) fgMCEvGen;
- fHijing->Initialize();
+
fHijing->SetIHPR2(3, fTrigger);
fHijing->SetIHPR2(4, fQuench);
fHijing->SetIHPR2(6, fShadowing);
fHijing->SetIHPR2(12, fDecaysOff);
fHijing->SetIHPR2(21, fKeep);
+ fHijing->Rluset(50,0);
+ fHijing->Initialize();
+
+
//
if (fEvaluate) EvaluateCrossSections();
}
if (np == 0 ) continue;
Int_t i;
Int_t * newPos = new Int_t[np];
+
for (i = 0; i<np; i++) *(newPos+i)=i;
-
+//
+// First write parent particles
+//
+
for (i = 0; i<np; i++) {
TParticle * iparticle = (TParticle *) particles->At(i);
-
+// Is this a parent particle ?
+ if (Stable(iparticle)) continue;
+//
Bool_t hasMother = (iparticle->GetFirstMother() >=0);
- Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
Bool_t selected = kTRUE;
Bool_t hasSelectedDaughters = kFALSE;
-
-
+
+
kf = iparticle->GetPdgCode();
+ ks = iparticle->GetStatusCode();
+ if (kf == 92) continue;
+
if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
- if (hasDaughter && !selected) hasSelectedDaughters = DaughtersSelection(iparticle, particles);
+ hasSelectedDaughters = DaughtersSelection(iparticle, particles);
//
// Put particle on the stack if it is either selected or it is the mother of at least one seleted particle
//
-
if (selected || hasSelectedDaughters) {
nc++;
- ks = iparticle->GetStatusCode();
p[0]=iparticle->Px();
p[1]=iparticle->Py();
p[2]=iparticle->Pz();
imo=-1;
if (hasMother) {
imo=iparticle->GetFirstMother();
- imo=*(newPos+imo);
+ TParticle* mother= (TParticle *) particles->At(imo);
+ imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
}
-
-// printf("\n selected iparent %d %d %d \n",i, kf, imo);
- if (hasDaughter) {
- gAlice->SetTrack(0,imo,kf,p,origin,polar,
- tof,"Primary",nt);
- } else {
- gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
- tof,"Secondary",nt);
+// Put particle on the stack ...
+// printf("\n set track mother: %d %d %d %d %d %d ",i,imo, kf, nt+1, selected, hasSelectedDaughters);
+
+ gAlice->SetTrack(0,imo,kf,p,origin,polar,
+ tof,"Primary",nt);
+// ... and keep it there
+ gAlice->KeepTrack(nt);
+//
+ *(newPos+i)=nt;
+ } // selected
+ } // particle loop parents
+//
+// Now write the final state particles
+//
+
+ for (i = 0; i<np; i++) {
+ TParticle * iparticle = (TParticle *) particles->At(i);
+// Is this a final state particle ?
+ if (!Stable(iparticle)) continue;
+//
+ Bool_t hasMother = (iparticle->GetFirstMother() >=0);
+ Bool_t selected = kTRUE;
+ kf = iparticle->GetPdgCode();
+ if (!fSelectAll) selected = KinematicSelection(iparticle)&&SelectFlavor(kf);
+//
+// Put particle on the stack if selected
+//
+ if (selected) {
+ nc++;
+ ks = iparticle->GetStatusCode();
+ p[0]=iparticle->Px();
+ p[1]=iparticle->Py();
+ p[2]=iparticle->Pz();
+ origin[0]=origin0[0]+iparticle->Vx()/10;
+ origin[1]=origin0[1]+iparticle->Vy()/10;
+ origin[2]=origin0[2]+iparticle->Vz()/10;
+ tof=kconv*iparticle->T();
+ imo=-1;
+
+ if (hasMother) {
+ imo=iparticle->GetFirstMother();
+ TParticle* mother= (TParticle *) particles->At(imo);
+ imo = (mother->GetPdgCode() != 92) ? imo=*(newPos+imo) : -1;
}
+// Put particle on the stack
+ gAlice->SetTrack(fTrackIt,imo,kf,p,origin,polar,
+ tof,"Secondary",nt);
+
+// printf("\n set track final: %d %d %d",imo, kf, nt);
+ gAlice->KeepTrack(nt);
*(newPos+i)=nt;
} // selected
- } // particle loop
+ } // particle loop final state
+
delete newPos;
+
printf("\n I've put %i particles on the stack \n",nc);
if (nc > 0) {
jev+=nc;
}
}
} // event loop
+ fHijing->Rluget(50,-1);
}
Bool_t AliGenHijing::KinematicSelection(TParticle *particle)
{
// Perform kinematic selection
- Float_t px=particle->Px();
- Float_t py=particle->Py();
- Float_t pz=particle->Pz();
- Float_t e=particle->Energy();
+ Double_t px=particle->Px();
+ Double_t py=particle->Py();
+ Double_t pz=particle->Pz();
+ Double_t e=particle->Energy();
//
// transverse momentum cut
- Float_t pt=TMath::Sqrt(px*px+py*py);
+ Double_t pt=TMath::Sqrt(px*px+py*py);
if (pt > fPtMax || pt < fPtMin)
{
// printf("\n failed pt cut %f %f %f \n",pt,fPtMin,fPtMax);
}
//
// momentum cut
- Float_t p=TMath::Sqrt(px*px+py*py+pz*pz);
+ Double_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);
//
// theta cut
- Float_t theta = Float_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
+ Double_t theta = Double_t(TMath::ATan2(Double_t(pt),Double_t(pz)));
if (theta > fThetaMax || theta < fThetaMin)
{
//
// rapidity cut
- Float_t y = 0.5*TMath::Log((e+pz)/(e-pz));
+ Double_t y;
+ if(e<=pz) y = 99;
+ else if (e<=-pz) y = -99;
+ else y = 0.5*TMath::Log((e+pz)/(e-pz));
if (y > fYMax || y < fYMin)
{
// printf("\n failed y cut %f %f %f \n",y,fYMin,fYMax);
//
// phi cut
- Float_t phi=Float_t(TMath::ATan2(Double_t(py),Double_t(px)));
+ Double_t phi=Double_t(TMath::ATan2(Double_t(py),Double_t(px)));
if (phi > fPhiMax || phi < fPhiMin)
{
// printf("\n failed phi cut %f %f %f \n",phi,fPhiMin,fPhiMax);
xPartHard+=gbh;
}
- if ((xTot-oldvalue)/oldvalue<0.0001) break;
+ if(oldvalue) if ((xTot-oldvalue)/oldvalue<0.0001) break;
oldvalue=xTot;
printf("\n Total cross section (barn): %d %f %f \n",i, xb, xTot);
printf("\n Hard cross section (barn): %d %f %f \n\n",i, xb, xTotHard);
for (i=imin; i<= imax; i++){
TParticle * jparticle = (TParticle *) particles->At(i);
Int_t ip=jparticle->GetPdgCode();
- if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {selected=kTRUE; break;}
+ if (KinematicSelection(jparticle)&&SelectFlavor(ip)) {
+ selected=kTRUE; break;
+ }
if (DaughtersSelection(jparticle, particles)) {selected=kTRUE; break; }
}
} else {
Int_t ifl=TMath::Abs(pid/100);
if (ifl > 10) ifl/=10;
- return ((fFlavor==4 && (ifl==4 || ifl==5)) ||
- (fFlavor==5 && ifl==5));
+ return (fFlavor == ifl);
+}
+Bool_t AliGenHijing::Stable(TParticle* particle)
+{
+ Int_t kf = TMath::Abs(particle->GetPdgCode());
+
+ if ( (particle->GetFirstDaughter() < 0 ) || (kf == 1000*fFlavor+122))
+
+ {
+ return kTRUE;
+ } else {
+ return kFALSE;
+ }
}
void AliGenHijing::MakeHeader()
// Assignment operator
return *this;
}
+
+
+
+
+
+
+
+
+
+
+
+
+
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff