// 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)
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),
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),
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),
- fFragPhotonInCalo(kFALSE),
- fPi0InCalo(kFALSE) ,
- fCheckEMCAL(kFALSE),
- fCheckPHOS(kFALSE),
- fFragPhotonOrPi0MinPt(0),
- fPHOSMinPhi(219.),
- fPHOSMaxPhi(321.),
- fPHOSEta(0.13),
- fEMCALMinPhi(79.),
- fEMCALMaxPhi(191.),
- fEMCALEta(0.71)
-{
-// 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:
//
Int_t i;
-
+ 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);
polar[0], polar[1], polar[2],
kPPrimary, nt, 1., ks);
pParent[i] = nt;
- KeepTrack(nt);
+ KeepTrack(nt);
+ fNprimaries++;
} // PushTrack loop
}
} else {
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 && iparticle->GetFirstMother() != 11))//No photon from hadron decay
+ (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
return 0;
}
+
+ // 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++) {
}
//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;
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
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 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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff