X-Git-Url: http://git.uio.no/git/?a=blobdiff_plain;f=THijing%2FAliGenHijing.cxx;h=645b05da9ed20883250541ff1719bd33d688d7cc;hb=fcbe44bff5d201f27d03627ac71e429d0abc3b1f;hp=b1167c5c0a1d29b380ee7a27654f57a9e6a5dcab;hpb=39192f0596cda8cf10b41283c032de8a8a83d919;p=u%2Fmrichter%2FAliRoot.git diff --git a/THijing/AliGenHijing.cxx b/THijing/AliGenHijing.cxx index b1167c5c0a1..645b05da9ed 100644 --- a/THijing/AliGenHijing.cxx +++ b/THijing/AliGenHijing.cxx @@ -18,10 +18,11 @@ // Generator using HIJING as an external generator // The main HIJING options are accessable for the user through this interface. // Uses the THijing implementation of TGenerator. +// Author: +// Andreas Morsch (andreas.morsch@cern.ch) // -// andreas.morsch@cern.ch -#include +#include #include #include #include @@ -30,74 +31,109 @@ #include "AliGenHijing.h" #include "AliGenHijingEventHeader.h" -#include "AliRun.h" #include "AliHijingRndm.h" +#include "AliLog.h" +#include "AliRun.h" ClassImp(AliGenHijing) AliGenHijing::AliGenHijing() - :AliGenMC() + :AliGenMC(), + fFrame("CMS"), + fMinImpactParam(0.), + fMaxImpactParam(5.), + fKeep(0), + fQuench(1), + fShadowing(1), + fDecaysOff(1), + fTrigger(0), + fEvaluate(0), + fSelectAll(0), + fFlavor(0), + 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) { -// Constructor - fParticles = 0; - fHijing = 0; - fDsigmaDb = 0; - fDnDb = 0; - AliHijingRndm::SetHijingRandom(GetRandom()); + // Constructor + fEnergyCMS = 5500.; + AliHijingRndm::SetHijingRandom(GetRandom()); } AliGenHijing::AliGenHijing(Int_t npart) - :AliGenMC(npart) + :AliGenMC(npart), + fFrame("CMS"), + fMinImpactParam(0.), + fMaxImpactParam(5.), + fKeep(0), + fQuench(1), + fShadowing(1), + fDecaysOff(1), + fTrigger(0), + fEvaluate(0), + fSelectAll(0), + fFlavor(0), + 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) { // Default PbPb collisions at 5. 5 TeV // + fEnergyCMS = 5500.; fName = "Hijing"; fTitle= "Particle Generator using HIJING"; - - SetEnergyCMS(); - SetImpactParameterRange(); - SetBoostLHC(); - SetJetEtaRange(); - SetJetPhiRange(); - - fKeep = 0; - fQuench = 1; - fShadowing = 1; - fTrigger = 0; - fDecaysOff = 1; - fEvaluate = 0; - fSelectAll = 0; - fFlavor = 0; - fSpectators = 1; - fDsigmaDb = 0; - fDnDb = 0; - fPtMinJet = -2.5; - fRadiation = 3; - // - SetSimpleJets(); - SetNoGammas(); -// - fParticles = new TClonesArray("TParticle",10000); +// // // Set random number generator AliHijingRndm::SetHijingRandom(GetRandom()); - fHijing = 0; - -} - -AliGenHijing::AliGenHijing(const AliGenHijing & hijing): - AliGenMC(hijing) -{ -// copy constructor } - AliGenHijing::~AliGenHijing() { // Destructor if ( fDsigmaDb) delete fDsigmaDb; if ( fDnDb) delete fDnDb; - delete fParticles; } void AliGenHijing::Init() @@ -151,6 +187,15 @@ void AliGenHijing::Init() fHijing->SetHIPR1(11, 2.5); } +// +// Heavy quarks +// + if (fNoHeavyQuarks) { + fHijing->SetIHPR2(49, 1); + } else { + fHijing->SetIHPR2(49, 0); + } + AliGenMC::Init(); @@ -174,40 +219,45 @@ void AliGenHijing::Generate() 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; - Int_t j, kf, ks, imo; + Int_t j, kf, ks, ksp, imo; kf = 0; fTrials = 0; + for (j = 0;j < 3; j++) origin0[j] = fOrigin[j]; if(fVertexSmear == kPerEvent) { Vertex(); for (j=0; j < 3; j++) origin0[j] = fVertex[j]; } + + Float_t sign = (fRandomPz && (Rndm() < 0.5))? -1. : 1.; while(1) { // Generate one event // -------------------------------------------------------------------------- - fSpecn = 0; - fSpecp = 0; + fProjectileSpecn = 0; + fProjectileSpecp = 0; + fTargetSpecn = 0; + fTargetSpecp = 0; // -------------------------------------------------------------------------- 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; @@ -221,7 +271,7 @@ void AliGenHijing::Generate() // Get event vertex // - TParticle * iparticle = (TParticle *) fParticles->At(0); + TParticle * iparticle = (TParticle *) fParticles.At(0); fVertex[0] = origin0[0]; fVertex[1] = origin0[1]; fVertex[2] = origin0[2]; @@ -231,7 +281,7 @@ void AliGenHijing::Generate() // 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; @@ -261,26 +311,31 @@ void AliGenHijing::Generate() // for (i = 0; iAt(i); + iparticle = (TParticle *) fParticles.At(i); // Is this a final state particle ? if (!Stable(iparticle)) continue; Bool_t selected = kTRUE; kf = iparticle->GetPdgCode(); ks = iparticle->GetStatusCode(); + ksp = iparticle->GetUniqueID(); // -------------------------------------------------------------------------- // Count spectator neutrons and protons - if(ks == 0 || ks == 1 || ks == 10 || ks == 11){ - if(kf == kNeutron) fSpecn += 1; - if(kf == kProton) fSpecp += 1; + if(ksp == 0 || ksp == 1){ + if(kf == kNeutron) fProjectileSpecn += 1; + if(kf == kProton) fProjectileSpecp += 1; + } + else if(ksp == 10 || ksp == 11){ + if(kf == kNeutron) fTargetSpecn += 1; + if(kf == kProton) fTargetSpecp += 1; } // -------------------------------------------------------------------------- // if (!fSelectAll) { selected = KinematicSelection(iparticle,0)&&SelectFlavor(kf); - if (!fSpectators && selected) selected = (ks != 0 && ks != 1 && ks != 10 - && ks != 11); + if (!fSpectators && selected) selected = (ksp != 0 && ksp != 1 && ksp != 10 + && ksp != 11); } // // Put particle on the stack if selected @@ -290,34 +345,42 @@ void AliGenHijing::Generate() pSelected[i] = 1; } // selected } // particle loop final state + // -// Write particles to stack +// Time of the interactions + Float_t tInt = 0.; + if (fPileUpTimeWindow > 0.) tInt = fPileUpTimeWindow * (2. * gRandom->Rndm() - 1.); + // +// Write particles to stack + for (i = 0; iAt(i); + iparticle = (TParticle *) fParticles.At(i); Bool_t hasMother = (iparticle->GetFirstMother() >=0); Bool_t hasDaughter = (iparticle->GetFirstDaughter() >=0); - if (pSelected[i]) { kf = iparticle->GetPdgCode(); ks = iparticle->GetStatusCode(); p[0] = iparticle->Px(); p[1] = iparticle->Py(); - p[2] = iparticle->Pz(); + p[2] = iparticle->Pz() * sign; 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(); +// tof = kconv * iparticle->T() + sign * origin0[2] / 3.e10; + tof = kconv * iparticle->T(); + if (fPileUpTimeWindow > 0.) tof += tInt; + imo = -1; TParticle* mother = 0; if (hasMother) { imo = iparticle->GetFirstMother(); - mother = (TParticle *) fParticles->At(imo); - imo = (mother->GetPdgCode() != 92) ? imo = newPos[imo] : -1; + 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); + PushTrack(tFlag,imo,kf,p,origin,polar,tof,kPNoProcess,nt, 1., ks); + fNprimaries++; KeepTrack(nt); newPos[i] = nt; } // if selected @@ -325,12 +388,12 @@ void AliGenHijing::Generate() 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; } } @@ -425,7 +488,7 @@ Bool_t AliGenHijing::DaughtersSelection(TParticle* iparticle) 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; @@ -462,7 +525,7 @@ Bool_t AliGenHijing::SelectFlavor(Int_t pid) return res; } -Bool_t AliGenHijing::Stable(TParticle* particle) +Bool_t AliGenHijing::Stable(TParticle* particle) const { // Return true for a stable particle // @@ -481,7 +544,7 @@ 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)->SetNProduced(fNprimaries); ((AliGenHijingEventHeader*) header)->SetImpactParameter(fHijing->GetHINT1(19)); ((AliGenHijingEventHeader*) header)->SetTotalEnergy(fHijing->GetEATT()); ((AliGenHijingEventHeader*) header)->SetHardScatters(fHijing->GetJATT()); @@ -490,7 +553,11 @@ void AliGenHijing::MakeHeader() fHijing->GetN01(), fHijing->GetN10(), fHijing->GetN11()); - ((AliGenHijingEventHeader*) header)->SetSpectators(fSpecn, fSpecp); + ((AliGenHijingEventHeader*) header)->SetSpectators(fProjectileSpecn, fProjectileSpecp, + fTargetSpecn,fTargetSpecp); + ((AliGenHijingEventHeader*) header)->SetReactionPlaneAngle(fHijing->GetHINT1(20)); + + // 4-momentum vectors of the triggered jets. // @@ -519,10 +586,11 @@ void AliGenHijing::MakeHeader() ((AliGenHijingEventHeader*) header)->SetTrials(fTrials); // Event Vertex header->SetPrimaryVertex(fVertex); - gAlice->SetGenEventHeader(header); + AddHeader(header); fCollisionGeometry = (AliGenHijingEventHeader*) header; } + Bool_t AliGenHijing::CheckTrigger() { // Check the kinematic trigger condition @@ -559,12 +627,12 @@ Bool_t AliGenHijing::CheckTrigger() } 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 ks = part->GetStatusCode(); - if (kf == 22 && ks == 40) { + Int_t ksp = part->GetUniqueID(); + if (kf == 22 && ksp == 40) { Float_t phi = part->Phi(); Float_t eta = part->Eta(); if (eta < fEtaMaxJet && @@ -579,16 +647,3 @@ Bool_t AliGenHijing::CheckTrigger() } // 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); -} -