// andreas.morsch@cern.ch
//
+#include <TClonesArray.h>
#include <TDatabasePDG.h>
#include <TParticle.h>
#include <TPDGCode.h>
#include "AliStack.h"
#include "AliRunLoader.h"
#include "AliMC.h"
-#include "pyquenCommon.h"
+#include "PyquenCommon.h"
ClassImp(AliGenPythia)
AliGenMC(),
fProcess(kPyCharm),
fStrucFunc(kCTEQ5L),
- fEnergyCMS(5500.),
fKineBias(0.),
fTrials(0),
fTrialsRun(0),
fHFoff(kFALSE),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerMultiplicity(0),
+ fTriggerMultiplicityEta(0),
fCountMode(kCountAll),
fHeader(0),
fRL(0),
- fFileName(0)
+ fFileName(0),
+ fFragPhotonInCalo(kFALSE),
+ fPi0InCalo(kFALSE) ,
+ fPhotonInCalo(kFALSE),
+ fCheckEMCAL(kFALSE),
+ fCheckPHOS(kFALSE),
+ fCheckPHOSeta(kFALSE),
+ fFragPhotonOrPi0MinPt(0),
+ fPhotonMinPt(0),
+ fPHOSMinPhi(219.),
+ fPHOSMaxPhi(321.),
+ fPHOSEta(0.13),
+ fEMCALMinPhi(79.),
+ fEMCALMaxPhi(191.),
+ fEMCALEta(0.71)
+
{
// Default Constructor
+ fEnergyCMS = 5500.;
SetNuclei(0,0);
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
:AliGenMC(npart),
fProcess(kPyCharm),
fStrucFunc(kCTEQ5L),
- fEnergyCMS(5500.),
fKineBias(0.),
fTrials(0),
fTrialsRun(0),
fHFoff(kFALSE),
fTriggerParticle(0),
fTriggerEta(0.9),
+ fTriggerMultiplicity(0),
+ fTriggerMultiplicityEta(0),
fCountMode(kCountAll),
fHeader(0),
fRL(0),
- fFileName(0)
+ fFileName(0),
+ fFragPhotonInCalo(kFALSE),
+ fPi0InCalo(kFALSE) ,
+ fPhotonInCalo(kFALSE),
+ fCheckEMCAL(kFALSE),
+ fCheckPHOS(kFALSE),
+ fCheckPHOSeta(kFALSE),
+ fFragPhotonOrPi0MinPt(0),
+ fPhotonMinPt(0),
+ fPHOSMinPhi(219.),
+ fPHOSMaxPhi(321.),
+ fPHOSEta(0.13),
+ fEMCALMinPhi(79.),
+ fEMCALMaxPhi(191.),
+ fEMCALEta(0.71)
{
// default charm production at 5. 5 TeV
// semimuonic decay
// structure function GRVHO
//
+ fEnergyCMS = 5500.;
fName = "Pythia";
fTitle= "Particle Generator using PYTHIA";
SetForceDecay();
// Set random number generator
if (!AliPythiaRndm::GetPythiaRandom())
AliPythiaRndm::SetPythiaRandom(GetRandom());
- fParticles = new TClonesArray("TParticle",1000);
SetNuclei(0,0);
}
-AliGenPythia::AliGenPythia(const AliGenPythia & Pythia)
- :AliGenMC(Pythia),
- fProcess(kPyCharm),
- fStrucFunc(kCTEQ5L),
- fEnergyCMS(5500.),
- fKineBias(0.),
- fTrials(0),
- fTrialsRun(0),
- fQ(0.),
- fX1(0.),
- fX2(0.),
- fEventTime(0.),
- fInteractionRate(0.),
- fTimeWindow(0.),
- fCurSubEvent(0),
- fEventsTime(0),
- fNev(0),
- fFlavorSelect(0),
- fXsection(0.),
- fPythia(0),
- fPtHardMin(0.),
- fPtHardMax(1.e4),
- fYHardMin(-1.e10),
- fYHardMax(1.e10),
- fGinit(kTRUE),
- fGfinal(kTRUE),
- fHadronisation(kTRUE),
- fNpartons(0),
- fReadFromFile(kFALSE),
- fQuench(kFALSE),
- fPtKick(1.),
- fFullEvent(kTRUE),
- fDecayer(new AliDecayerPythia()),
- fDebugEventFirst(-1),
- fDebugEventLast(-1),
- fEtMinJet(0.),
- fEtMaxJet(1.e4),
- fEtaMinJet(-20.),
- fEtaMaxJet(20.),
- fPhiMinJet(0.),
- fPhiMaxJet(2.* TMath::Pi()),
- fJetReconstruction(kCell),
- fEtaMinGamma(-20.),
- fEtaMaxGamma(20.),
- fPhiMinGamma(0.),
- fPhiMaxGamma(2. * TMath::Pi()),
- fPycellEtaMax(2.),
- fPycellNEta(274),
- fPycellNPhi(432),
- fPycellThreshold(0.),
- fPycellEtSeed(4.),
- fPycellMinEtJet(10.),
- fPycellMaxRadius(1.),
- fStackFillOpt(kFlavorSelection),
- fFeedDownOpt(kTRUE),
- fFragmentation(kTRUE),
- fSetNuclei(kFALSE),
- fNewMIS(kFALSE),
- fHFoff(kFALSE),
- fTriggerParticle(0),
- fTriggerEta(0.9),
- fCountMode(kCountAll),
- fHeader(0),
- fRL(0),
- fFileName(0)
-{
-// copy constructor
- Pythia.Copy(*this);
-}
-
AliGenPythia::~AliGenPythia()
{
// Destructor
//
fParentWeight=1./Float_t(fNpart);
//
-// Forward Paramters to the AliPythia object
- fDecayer->SetForceDecay(fForceDecay);
- fDecayer->Init();
fPythia->SetCKIN(3,fPtHardMin);
// pt - kick
if (fPtKick > 0.) {
fPythia->SetMSTP(91,1);
- fPythia->SetPARP(91,fPtKick);
+ fPythia->SetPARP(91,fPtKick);
+ fPythia->SetPARP(93, 4. * fPtKick);
} else {
fPythia->SetMSTP(91,0);
}
} else {
fRL = 0x0;
}
+ //
+ fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
+ // Forward Paramters to the AliPythia object
+ fDecayer->SetForceDecay(fForceDecay);
// Switch off Heavy Flavors on request
if (fHFoff) {
+ // Maximum number of quark flavours used in pdf
fPythia->SetMSTP(58, 3);
+ // Maximum number of flavors that can be used in showers
fPythia->SetMSTJ(45, 3);
- for (Int_t i = 156; i <= 160; i++) fPythia->SetMDME(i, 1, 0);
+ // Switch off g->QQbar splitting in decay table
+ ((AliDecayerPythia*) fDecayer)->HeavyFlavourOff();
}
- //
- fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
+
+ fDecayer->Init();
+
// Parent and Children Selection
switch (fProcess)
case kPyJpsi:
fParentSelect[0] = 443;
break;
+ case kPyMbDefault:
case kPyMb:
case kPyMbNonDiffr:
case kPyMbMSEL1:
// printf("Calling hadronisation %d\n", fPythia->GetN());
fPythia->Pyexec();
fTrials++;
- fPythia->ImportParticles(fParticles,"All");
+ fPythia->ImportParticles(&fParticles,"All");
Boost();
//
//
//
Int_t i;
-
- Int_t np = fParticles->GetEntriesFast();
+ fNprimaries = 0;
+ Int_t np = fParticles.GetEntriesFast();
if (np == 0) continue;
//
Int_t nc = 0; // Total n. of selected particles
Int_t nParents = 0; // Selected parents
Int_t nTkbles = 0; // Trackable particles
- if (fProcess != kPyMb && fProcess != kPyJets &&
+ if (fProcess != kPyMbDefault &&
+ fProcess != kPyMb &&
+ fProcess != kPyJets &&
fProcess != kPyDirectGamma &&
fProcess != kPyMbNonDiffr &&
fProcess != kPyMbMSEL1 &&
- fProcess != kPyW && fProcess != kPyZ &&
- fProcess != kPyCharmppMNRwmi && fProcess != kPyBeautyppMNRwmi) {
+ fProcess != kPyW &&
+ fProcess != kPyZ &&
+ fProcess != kPyCharmppMNRwmi &&
+ fProcess != kPyBeautyppMNRwmi) {
for (i = 0; i < np; i++) {
- TParticle* iparticle = (TParticle *) fParticles->At(i);
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
Int_t ks = iparticle->GetStatusCode();
kf = CheckPDGCode(iparticle->GetPdgCode());
// No initial state partons
if (kf >= fFlavorSelect && kf <= 6) {
Int_t idau = iparticle->GetFirstDaughter() - 1;
if (idau > -1) {
- TParticle* daughter = (TParticle *) fParticles->At(idau);
+ TParticle* daughter = (TParticle *) fParticles.At(idau);
Int_t pdgD = daughter->GetPdgCode();
if (pdgD == 91 || pdgD == 92) {
Int_t jmin = daughter->GetFirstDaughter() - 1;
Int_t jmax = daughter->GetLastDaughter() - 1;
- for (Int_t j = jmin; j <= jmax; j++)
- ((TParticle *) fParticles->At(j))->SetFirstMother(i+1);
+ for (Int_t jp = jmin; jp <= jmax; jp++)
+ ((TParticle *) fParticles.At(jp))->SetFirstMother(i+1);
} // is string or cluster
} // has daughter
} // heavy quark
if (ipa > -1) {
- TParticle * mother = (TParticle *) fParticles->At(ipa);
+ TParticle * mother = (TParticle *) fParticles.At(ipa);
kfMo = TMath::Abs(mother->GetPdgCode());
}
if (nc > 0) {
for (i = 0; i<np; i++) {
if (!pSelected[i]) continue;
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
kf = CheckPDGCode(iparticle->GetPdgCode());
Int_t ks = iparticle->GetStatusCode();
p[0] = iparticle->Px();
polar[0], polar[1], polar[2],
kPPrimary, nt, 1., ks);
pParent[i] = nt;
- KeepTrack(nt);
+ KeepTrack(nt);
+ fNprimaries++;
} // PushTrack loop
}
} else {
- Int_t np = (fHadronisation) ? fParticles->GetEntriesFast() : fNpartons;
+ Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
Int_t* pParent = new Int_t[np];
for (i=0; i< np; i++) pParent[i] = -1;
if (fProcess == kPyJets || fProcess == kPyDirectGamma) {
- TParticle* jet1 = (TParticle *) fParticles->At(6);
- TParticle* jet2 = (TParticle *) fParticles->At(7);
+ TParticle* jet1 = (TParticle *) fParticles.At(6);
+ TParticle* jet2 = (TParticle *) fParticles.At(7);
if (!CheckTrigger(jet1, jet2)) {
delete [] pParent;
return 0;
}
}
+ // Select jets with fragmentation photon or pi0 going to PHOS or EMCAL
+ if (fProcess == kPyJets && (fFragPhotonInCalo || fPi0InCalo) ) {
+
+ Bool_t ok = kFALSE;
+
+ Int_t pdg = 0;
+ if (fFragPhotonInCalo) pdg = 22 ; // Photon
+ else if (fPi0InCalo) pdg = 111 ; // Pi0
+
+ for (i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
+ if(iparticle->GetStatusCode()==1 && iparticle->GetPdgCode()==pdg &&
+ iparticle->Pt() > fFragPhotonOrPi0MinPt){
+ Int_t imother = iparticle->GetFirstMother() - 1;
+ TParticle* pmother = (TParticle *) fParticles.At(imother);
+ if(pdg == 111 ||
+ (pdg == 22 && pmother->GetStatusCode() != 11))//No photon from hadron decay
+ {
+ Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
+ if((fCheckEMCAL && IsInEMCAL(phi,eta)) ||
+ (fCheckPHOS && IsInPHOS(phi,eta)) )
+ ok =kTRUE;
+ }
+ }
+ }
+ if(!ok)
+ return 0;
+ }
+
+ // Check for minimum multiplicity
+ if (fTriggerMultiplicity > 0) {
+ Int_t multiplicity = 0;
+ for (i = 0; i < np; i++) {
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
+
+ Int_t statusCode = iparticle->GetStatusCode();
+
+ // Initial state particle
+ if (statusCode > 20)
+ continue;
+
+ // skip quarks and gluons
+ Int_t pdgCode = TMath::Abs(iparticle->GetPdgCode());
+ if (pdgCode <= 10 || pdgCode == 21)
+ continue;
+
+ if (fTriggerMultiplicityEta > 0 && TMath::Abs(iparticle->Eta()) > fTriggerMultiplicityEta)
+ continue;
+
+ TParticlePDG* pdgPart = iparticle->GetPDG();
+ if (pdgPart && pdgPart->Charge() == 0)
+ continue;
+
+ ++multiplicity;
+ }
+
+ if (multiplicity < fTriggerMultiplicity) {
+ delete [] pParent;
+ return 0;
+ }
+
+ Printf("Triggered on event with multiplicity of %d > %d", multiplicity, fTriggerMultiplicity);
+ }
+
+ // Select events with a photon pt > min pt going to PHOS eta acceptance or exactly PHOS eta phi
+ if ((fProcess == kPyJets || fProcess == kPyDirectGamma) && fPhotonInCalo && (fCheckPHOSeta || fCheckPHOS)){
+
+ Bool_t okd = kFALSE;
+
+ Int_t pdg = 22;
+ Int_t iphcand = -1;
+ for (i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
+ Float_t phi = iparticle->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(iparticle->Eta());//in calos etamin=-etamax
+
+ if(iparticle->GetStatusCode() == 1
+ && iparticle->GetPdgCode() == pdg
+ && iparticle->Pt() > fPhotonMinPt
+ && eta < fPHOSEta){
+
+ // first check if the photon is in PHOS phi
+ if(IsInPHOS(phi,eta)){
+ okd = kTRUE;
+ break;
+ }
+ if(fCheckPHOSeta) iphcand = i; // candiate photon to rotate in phi
+
+ }
+ }
+
+ if(!okd && iphcand != -1) // execute rotation in phi
+ RotatePhi(iphcand,okd);
+
+ if(!okd)
+ return 0;
+ }
+
if (fTriggerParticle) {
Bool_t triggered = kFALSE;
for (i = 0; i < np; i++) {
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
kf = CheckPDGCode(iparticle->GetPdgCode());
if (kf != fTriggerParticle) continue;
if (iparticle->Pt() == 0.) continue;
Float_t yQ;
Int_t pdgQ;
for(i=0; i<np; i++) {
- hvq = (TParticle*)fParticles->At(i);
+ hvq = (TParticle*)fParticles.At(i);
pdgQ = hvq->GetPdgCode();
if(TMath::Abs(pdgQ) != fFlavorSelect) continue;
if(pdgQ>0) { theQ=kTRUE; } else { theQbar=kTRUE; }
}
//Introducing child cuts in case kPyW, kPyZ, kPyMb, and kPyMbNonDiff
- if ( (fProcess == kPyW || fProcess == kPyZ || fProcess == kPyMb || fProcess == kPyMbNonDiffr)
+ if ( (fProcess == kPyW ||
+ fProcess == kPyZ ||
+ fProcess == kPyMbDefault ||
+ fProcess == kPyMb ||
+ fProcess == kPyMbNonDiffr)
&& (fCutOnChild == 1) ) {
if ( !CheckKinematicsOnChild() ) {
delete[] pParent;
for (i = 0; i < np; i++) {
Int_t trackIt = 0;
- TParticle * iparticle = (TParticle *) fParticles->At(i);
+ TParticle * iparticle = (TParticle *) fParticles.At(i);
kf = CheckPDGCode(iparticle->GetPdgCode());
Int_t ks = iparticle->GetStatusCode();
Int_t km = iparticle->GetFirstMother();
origin[0], origin[1], origin[2], tof,
polar[0], polar[1], polar[2],
kPPrimary, nt, 1., ks);
+ fNprimaries++;
//
// Special Treatment to store color-flow
//
//
// Event Vertex
fHeader->SetPrimaryVertex(fVertex);
+
+//
+// Number of primaries
+ fHeader->SetNProduced(fNprimaries);
//
// Jets that have triggered
// Pass header
//
AddHeader(fHeader);
+ fHeader = 0x0;
}
-void AliGenPythia::AddHeader(AliGenEventHeader* header)
-{
- // Add header to container or runloader
- if (fContainer) {
- fContainer->AddHeader(header);
- } else {
- AliRunLoader::GetRunLoader()->GetHeader()->SetGenEventHeader(header);
- }
-}
-
-
Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2)
{
// Check the kinematic trigger condition
Int_t nPartAcc = 0; //number of particles in the acceptance range
Int_t numberOfAcceptedParticles = 1;
if (fNumberOfAcceptedParticles != 0) { numberOfAcceptedParticles = fNumberOfAcceptedParticles; }
- Int_t npart = fParticles->GetEntriesFast();
+ Int_t npart = fParticles.GetEntriesFast();
for (j = 0; j<npart; j++) {
- TParticle * jparticle = (TParticle *) fParticles->At(j);
+ TParticle * jparticle = (TParticle *) fParticles.At(j);
kcode = TMath::Abs( CheckPDGCode(jparticle->GetPdgCode()) );
ks = jparticle->GetStatusCode();
km = jparticle->GetFirstMother();
return checking;
}
-
-AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs)
-{
-// Assignment operator
- rhs.Copy(*this);
- return *this;
-}
-
void AliGenPythia::LoadEvent(AliStack* stack, Int_t flag, Int_t reHadr)
{
//
return;
}
+Bool_t AliGenPythia::IsInEMCAL(Float_t phi, Float_t eta)
+{
+ // Is particle in EMCAL acceptance?
+ // phi in degrees, etamin=-etamax
+ if(phi > fEMCALMinPhi && phi < fEMCALMaxPhi &&
+ eta < fEMCALEta )
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
+Bool_t AliGenPythia::IsInPHOS(Float_t phi, Float_t eta)
+{
+ // Is particle in PHOS acceptance?
+ // Acceptance slightly larger considered.
+ // phi in degrees, etamin=-etamax
+ if(phi > fPHOSMinPhi && phi < fPHOSMaxPhi &&
+ eta < fPHOSEta )
+ return kTRUE;
+ else
+ return kFALSE;
+}
+
+void AliGenPythia::RotatePhi(Int_t iphcand, Bool_t& okdd)
+{
+ //calculate the new position random between fPHOSMinPhi and fPHOSMaxPhi
+ Double_t phiPHOSmin = TMath::Pi()*fPHOSMinPhi/180;
+ Double_t phiPHOSmax = TMath::Pi()*fPHOSMaxPhi/180;
+ Double_t phiPHOS = gRandom->Uniform(phiPHOSmin,phiPHOSmax);
+
+ //calculate deltaphi
+ TParticle* ph = (TParticle *) fParticles.At(iphcand);
+ Double_t phphi = ph->Phi();
+ Double_t deltaphi = phiPHOS - phphi;
+
+
+
+ //loop for all particles and produce the phi rotation
+ Int_t np = (fHadronisation) ? fParticles.GetEntriesFast() : fNpartons;
+ Double_t oldphi, newphi;
+ Double_t newVx, newVy, R, Vz, time;
+ Double_t newPx, newPy, pt, Pz, e;
+ for(Int_t i=0; i< np; i++) {
+ TParticle* iparticle = (TParticle *) fParticles.At(i);
+ oldphi = iparticle->Phi();
+ newphi = oldphi + deltaphi;
+ if(newphi < 0) newphi = 2*TMath::Pi() + newphi; // correct angle
+ if(newphi > 2*TMath::Pi()) newphi = newphi - 2*TMath::Pi(); // correct angle
+
+ R = iparticle->R();
+ newVx = R*TMath::Cos(newphi);
+ newVy = R*TMath::Sin(newphi);
+ Vz = iparticle->Vz(); // don't transform
+ time = iparticle->T(); // don't transform
+
+ pt = iparticle->Pt();
+ newPx = pt*TMath::Cos(newphi);
+ newPy = pt*TMath::Sin(newphi);
+ Pz = iparticle->Pz(); // don't transform
+ e = iparticle->Energy(); // don't transform
+
+ // apply rotation
+ iparticle->SetProductionVertex(newVx, newVy, Vz, time);
+ iparticle->SetMomentum(newPx, newPy, Pz, e);
+
+ } //end particle loop
+
+ // now let's check that we put correctly the candidate photon in PHOS
+ Float_t phi = ph->Phi()*180./TMath::Pi(); //Convert to degrees
+ Float_t eta =TMath::Abs(ph->Eta());//in calos etamin=-etamax
+ if(IsInPHOS(phi,eta))
+ okdd = kTRUE;
+}
+
+
#ifdef never
void AliGenPythia::Streamer(TBuffer &R__b)
{
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff