// Andreas Morsch (andreas.morsch@cern.ch)
//
+#include <TClonesArray.h>
#include <TGraph.h>
#include <THijing.h>
#include <TLorentzVector.h>
#include "AliGenHijing.h"
#include "AliGenHijingEventHeader.h"
-#include "AliRun.h"
#include "AliHijingRndm.h"
+#include "AliLog.h"
+#include "AliRun.h"
ClassImp(AliGenHijing)
fKeep(0),
fQuench(1),
fShadowing(1),
- fDecaysOff(1),
+ fDecaysOff(3),
fTrigger(0),
fEvaluate(0),
fSelectAll(0),
fFlavor(0),
- fEnergyCMS(5500.),
fKineBias(0.),
fTrials(0),
fXsection(0.),
fHijing(0),
- fPtHardMin(0.),
- fPtHardMax(1.e4),
+ fPtHardMin(2.0),
+ fPtHardMax(-1),
fSpectators(1),
fDsigmaDb(0),
fDnDb(0),
fTargetSpecp(0),
fLHC(kFALSE),
fRandomPz(kFALSE),
- fNoHeavyQuarks(kFALSE)
+ fNoHeavyQuarks(kFALSE),
+ fHeader(AliGenHijingEventHeader("Hijing")),
+ fSigmaNN(-1),
+ fNoElas(0)
{
-// Constructor
- fParticles = 0;
- AliHijingRndm::SetHijingRandom(GetRandom());
+ // Constructor
+ fEnergyCMS = 5500.;
+ AliHijingRndm::SetHijingRandom(GetRandom());
}
AliGenHijing::AliGenHijing(Int_t npart)
fKeep(0),
fQuench(1),
fShadowing(1),
- fDecaysOff(1),
+ fDecaysOff(3),
fTrigger(0),
fEvaluate(0),
fSelectAll(0),
fFlavor(0),
- fEnergyCMS(5500.),
fKineBias(0.),
fTrials(0),
fXsection(0.),
fHijing(0),
- fPtHardMin(0.),
- fPtHardMax(1.e4),
+ fPtHardMin(2.0),
+ fPtHardMax(-1),
fSpectators(1),
fDsigmaDb(0),
fDnDb(0),
fTargetSpecp(0),
fLHC(kFALSE),
fRandomPz(kFALSE),
- fNoHeavyQuarks(kFALSE)
+ fNoHeavyQuarks(kFALSE),
+ fHeader(AliGenHijingEventHeader("Hijing")),
+ fSigmaNN(-1),
+ fNoElas(0)
{
// Default PbPb collisions at 5. 5 TeV
//
+ fEnergyCMS = 5500.;
fName = "Hijing";
fTitle= "Particle Generator using HIJING";
//
- fParticles = new TClonesArray("TParticle",10000);
//
// Set random number generator
AliHijingRndm::SetHijingRandom(GetRandom());
+
}
-AliGenHijing::AliGenHijing(const AliGenHijing & hijing):
- AliGenMC(hijing),
- fFrame("CMS"),
- fMinImpactParam(0.),
- fMaxImpactParam(5.),
- fKeep(0),
- fQuench(1),
- fShadowing(1),
- fDecaysOff(1),
- fTrigger(0),
- fEvaluate(0),
- fSelectAll(0),
- fFlavor(0),
- fEnergyCMS(5500.),
- fKineBias(0.),
- fTrials(0),
- fXsection(0.),
- fHijing(0),
- fPtHardMin(0.),
- fPtHardMax(1.e4),
- fSpectators(1),
- fDsigmaDb(0),
- fDnDb(0),
- fPtMinJet(-2.5),
- fEtaMinJet(-20.),
- fEtaMaxJet(+20.),
- fPhiMinJet(0.),
- fPhiMaxJet(2. * TMath::Pi()),
- fRadiation(3),
- fSimpleJet(kFALSE),
- fNoGammas(kFALSE),
- fProjectileSpecn(0),
- fProjectileSpecp(0),
- fTargetSpecn(0),
- fTargetSpecp(0),
- fLHC(kFALSE),
- fRandomPz(kFALSE),
- fNoHeavyQuarks(kFALSE)
-{
-// copy constructor
-}
-
-
AliGenHijing::~AliGenHijing()
{
// Destructor
if ( fDsigmaDb) delete fDsigmaDb;
if ( fDnDb) delete fDnDb;
- delete fParticles;
}
void AliGenHijing::Init()
fHijing->SetIHPR2(6, fShadowing);
fHijing->SetIHPR2(12, fDecaysOff);
fHijing->SetIHPR2(21, fKeep);
- fHijing->SetHIPR1(10, fPtMinJet);
+ fHijing->SetHIPR1(8, fPtHardMin);
+ fHijing->SetHIPR1(9, fPtHardMax);
+ fHijing->SetHIPR1(10, fPtMinJet);
+ if (fSigmaNN>0)
+ fHijing->SetHIPR1(31, fSigmaNN/2.);
fHijing->SetHIPR1(50, fSimpleJet);
+ //
+ // Switching off elastic scattering
+ if (fNoElas)
+ fHijing->SetIHPR2(14, 0);
//
// Quenching
//
//
}
+void AliGenHijing::SetSeed(UInt_t seed)
+{
+ AliHijingRndm::GetHijingRandom()->SetSeed(seed);
+}
+
void AliGenHijing::Generate()
{
// Generate one event
Float_t polar[3] = {0,0,0};
Float_t origin[3] = {0,0,0};
Float_t origin0[3] = {0,0,0};
+ Float_t time0 = 0.;
Float_t p[3];
Float_t tof;
// converts from mm/c to s
- const Float_t kconv = 0.001/2.999792458e8;
+ const Float_t kconv = 0.001/2.99792458e8;
//
Int_t nt = 0;
Int_t jev = 0;
fTrials = 0;
+
for (j = 0;j < 3; j++) origin0[j] = fOrigin[j];
+ time0 = fTimeOrigin;
+
if(fVertexSmear == kPerEvent) {
Vertex();
for (j=0; j < 3; j++) origin0[j] = fVertex[j];
+ time0 = fTime;
}
Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.;
+
while(1)
{
// Generate one event
// --------------------------------------------------------------------------
fHijing->GenerateEvent();
fTrials++;
- fHijing->ImportParticles(fParticles,"All");
+ fNprimaries = 0;
+ fHijing->ImportParticles(&fParticles,"All");
if (fTrigger != kNoTrigger) {
if (!CheckTrigger()) continue;
}
if (fLHC) Boost();
- Int_t np = fParticles->GetEntriesFast();
- printf("\n **************************************************%d\n",np);
+ Int_t np = fParticles.GetEntriesFast();
Int_t nc = 0;
if (np == 0 ) continue;
Int_t i;
// Get event vertex
//
- TParticle * iparticle = (TParticle *) fParticles->At(0);
fVertex[0] = origin0[0];
fVertex[1] = origin0[1];
fVertex[2] = origin0[2];
-
+ fTime = time0;
//
// First select parent particles
//
-
+ TParticle * iparticle = 0;
for (i = 0; i < np; i++) {
- iparticle = (TParticle *) fParticles->At(i);
+ iparticle = (TParticle *) fParticles.At(i);
// Is this a parent particle ?
if (Stable(iparticle)) continue;
//
for (i = 0; i<np; i++) {
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ iparticle = (TParticle *) fParticles.At(i);
// Is this a final state particle ?
if (!Stable(iparticle)) continue;
} // selected
} // particle loop final state
-//
-// Time of the interactions
- Float_t tInt = 0.;
- if (fPileUpTimeWindow > 0.) tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.);
-
//
// Write particles to stack
for (i = 0; i<np; i++) {
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ iparticle = (TParticle *) fParticles.At(i);
Bool_t hasMother = (iparticle->GetFirstMother() >=0);
Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0);
if (pSelected[i]) {
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() + sign * origin0[2] / 3.e10;
- if (fPileUpTimeWindow > 0.) tof += tInt;
-
+ tof = time0+kconv * iparticle->T();
+
imo = -1;
TParticle* mother = 0;
if (hasMother) {
imo = iparticle->GetFirstMother();
- mother = (TParticle *) fParticles->At(imo);
+ mother = (TParticle *) fParticles.At(imo);
imo = (mother->GetPdgCode() != 92) ? newPos[imo] : -1;
} // if has mother
Bool_t tFlag = (fTrackIt && !hasDaughter);
PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks);
-
-
+ fNprimaries++;
KeepTrack(nt);
newPos[i] = nt;
} // if selected
delete[] newPos;
delete[] pSelected;
- printf("\n I've put %i particles on the stack \n",nc);
+ AliInfo(Form("\n I've put %i particles on the stack \n",nc));
if (nc > 0) {
jev += nc;
if (jev >= fNpart || fNpart == -1) {
fKineBias = Float_t(fNpart)/Float_t(fTrials);
- printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+ AliInfo(Form("\n Trials: %i %i %i\n",fTrials, fNpart, jev));
break;
}
}
} // event loop
+
MakeHeader();
SetHighWaterMark(nt);
}
printf("\n Projectile Radius (fm): %f \n",fHijing->GetHIPR1(34));
printf("\n Target Radius (fm): %f \n",fHijing->GetHIPR1(35));
+ printf("\n Inelastic and total cross section (mb) %f %f \n",fHijing->GetHINT1(12), fHijing->GetHINT1(13));
Int_t i;
Float_t oldvalue= 0.;
Float_t* b = new Float_t[kMax];
Float_t* si1 = new Float_t[kMax];
Float_t* si2 = new Float_t[kMax];
-
+ for (i = 0; i < kMax; i++){
+ b[i] = 0.;
+ si1[i] = 0.;
+ si2[i] = 0.;
+ }
+
for (i = 0; i < kMax; i++)
{
Float_t xb = bMin+i*kdib;
fDnDb = new TGraph(i, b, si2);
}
-Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle)
+Bool_t AliGenHijing::DaughtersSelection(const TParticle* iparticle)
{
//
// Looks recursively if one of the daughters has been selected
imin = iparticle->GetFirstDaughter();
imax = iparticle->GetLastDaughter();
for (i = imin; i <= imax; i++){
- TParticle * jparticle = (TParticle *) fParticles->At(i);
+ TParticle * jparticle = (TParticle *) fParticles.At(i);
Int_t ip = jparticle->GetPdgCode();
if (KinematicSelection(jparticle,0)&&SelectFlavor(ip)) {
selected=kTRUE; break;
return res;
}
-Bool_t AliGenHijing::Stable(TParticle* particle) const
+Bool_t AliGenHijing::Stable(const TParticle* particle) const
{
// Return true for a stable particle
//
void AliGenHijing::MakeHeader()
{
// Builds the event header, to be called after each event
- AliGenEventHeader* header = new AliGenHijingEventHeader("Hijing");
- ((AliGenHijingEventHeader*) header)->SetNProduced(fHijing->GetNATT());
- ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19));
- ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT());
- ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT());
- ((AliGenHijingEventHeader*) header)->SetParticipants(fHijing->GetNP(), fHijing->GetNT());
- ((AliGenHijingEventHeader*) header)->SetCollisions(fHijing->GetN0(),
- fHijing->GetN01(),
- fHijing->GetN10(),
- fHijing->GetN11());
- ((AliGenHijingEventHeader*) header)->SetSpectators(fProjectileSpecn, fProjectileSpecp,
- fTargetSpecn,fTargetSpecp);
- ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(fHijing->GetHINT1(20));
-
-
+ fHeader.SetNProduced(fNprimaries);
+ fHeader.SetImpactParameter(fHijing->GetHINT1(19));
+ fHeader.SetTotalEnergy(fHijing->GetEATT());
+ fHeader.SetHardScatters(fHijing->GetJATT());
+ fHeader.SetParticipants(fHijing->GetNP(), fHijing->GetNT());
+ fHeader.SetCollisions(fHijing->GetN0(),
+ fHijing->GetN01(),
+ fHijing->GetN10(),
+ fHijing->GetN11());
+ fHeader.SetSpectators(fProjectileSpecn, fProjectileSpecp,
+ fTargetSpecn,fTargetSpecp);
+ fHeader.SetReactionPlaneAngle(fHijing->GetHINT1(20));
+ fHeader.SetTrueNPart(fHijing->GetNPART());
// 4-momentum vectors of the triggered jets.
//
fHijing->GetHINT1(37),
fHijing->GetHINT1(38),
fHijing->GetHINT1(39));
- ((AliGenHijingEventHeader*) header)->SetJets(jet1, jet2, jet3, jet4);
+ fHeader.SetJets(jet1, jet2, jet3, jet4);
// Bookkeeping for kinematic bias
- ((AliGenHijingEventHeader*) header)->SetTrials(fTrials);
+ fHeader.SetTrials(fTrials);
// Event Vertex
- header->SetPrimaryVertex(fVertex);
- AddHeader(header);
- fCollisionGeometry = (AliGenHijingEventHeader*) header;
-}
-
-void AliGenHijing::AddHeader(AliGenEventHeader* header)
-{
- // Passes header either to the container or to gAlice
- if (fContainer) {
- fContainer->AddHeader(header);
- } else {
- gAlice->SetGenEventHeader(header);
+ fHeader.SetPrimaryVertex(fVertex);
+ fHeader.SetInteractionTime(fTime);
+
+ Int_t nsd1 = 0,nsd2 = 0,ndd = 0;
+ Int_t nT = fHijing->GetNT();
+ Int_t nP = fHijing->GetNP();
+ for (Int_t i = 1; i <= nP; ++i) {
+ for (Int_t j = 1; j <= nT; ++j) {
+ Int_t tp = fHijing->GetNFP(i, 5);
+ Int_t tt = fHijing->GetNFT(j, 5);
+ if (tp == 2)
+ nsd1++;
+ if (tt == 2)
+ nsd2++;
+ if (tp == 2 && tt == 2)
+ ndd++;
+ }
}
+ fHeader.SetNDiffractive(nsd1, nsd2, ndd);
+ AddHeader(&fHeader);
+ fCollisionGeometry = &fHeader;
}
} else if (fTrigger == 2) {
// Gamma Jet
//
- Int_t np = fParticles->GetEntriesFast();
+ Int_t np = fParticles.GetEntriesFast();
for (Int_t i = 0; i < np; i++) {
- TParticle* part = (TParticle*) fParticles->At(i);
+ TParticle* part = (TParticle*) fParticles.At(i);
Int_t kf = part->GetPdgCode();
Int_t ksp = part->GetUniqueID();
if (kf == 22 && ksp == 40) {
} // fTrigger == 2
return triggered;
}
-
-
-void AliGenHijing::Copy(TObject &) const
-{
- Fatal("Copy","Not implemented!\n");
-}
-
-AliGenHijing& AliGenHijing::operator=(const AliGenHijing& rhs)
-{
- rhs.Copy(*this);
- return (*this);
-}
-