From cc363573d939674632d1b0c00e457d0df0caf063 Mon Sep 17 00:00:00 2001 From: morsch Date: Fri, 15 Nov 2002 00:43:06 +0000 Subject: [PATCH] Changes for kPyJets - initial and final state g-radiation + pt-kick default - trigger based on parton clusters (using pyclus) - trigger jets are stored in header. --- EVGEN/AliGenPythia.cxx | 150 ++++++++++++++++++++++++++++++++++++----- EVGEN/AliGenPythia.h | 21 +++++- 2 files changed, 150 insertions(+), 21 deletions(-) diff --git a/EVGEN/AliGenPythia.cxx b/EVGEN/AliGenPythia.cxx index 0820dcb61eb..e603c18dd4d 100644 --- a/EVGEN/AliGenPythia.cxx +++ b/EVGEN/AliGenPythia.cxx @@ -15,6 +15,9 @@ /* $Log$ +Revision 1.63 2002/10/14 14:55:35 hristov +Merging the VirtualMC branch to the main development branch (HEAD) + Revision 1.52.4.4 2002/10/10 16:40:08 hristov Updating VirtualMC to v3-09-02 @@ -203,8 +206,11 @@ Introduction of the Copyright and cvs Log #include "AliRun.h" #include "AliPythia.h" #include "AliPDG.h" +#include "AliConst.h" #include #include +#include +#include ClassImp(AliGenPythia) @@ -218,6 +224,7 @@ AliGenPythia::AliGenPythia() SetEventListRange(); SetJetPhiRange(); SetJetEtaRange(); + SetJetEtRange(); SetGammaPhiRange(); SetGammaEtaRange(); } @@ -250,6 +257,7 @@ AliGenPythia::AliGenPythia(Int_t npart) SetEventListRange(); SetJetPhiRange(); SetJetEtaRange(); + SetJetEtRange(); SetGammaPhiRange(); SetGammaEtaRange(); // Options determining what to keep in the stack (Heavy flavour generation) @@ -284,7 +292,8 @@ void AliGenPythia::SetEventListRange(Int_t eventFirst, Int_t eventLast) void AliGenPythia::Init() { // Initialisation - SetMC(AliPythia::Instance()); + + SetMC(AliPythia::Instance()); fPythia=(AliPythia*) fgMCEvGen; // fParentWeight=1./Float_t(fNpart); @@ -309,8 +318,15 @@ void AliGenPythia::Init() fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc); - // fPythia->Pylist(0); - // fPythia->Pystat(2); +// initial state radiation + fPythia->SetMSTP(61,fGinit); +// final state radiation + fPythia->SetMSTP(71,fGfinal); +// + fPythia->SetMSTP(91,1); + fPythia->SetMSTU(16,1); + fPythia->SetMSTJ(1,1); +// // Parent and Children Selection switch (fProcess) { @@ -705,18 +721,36 @@ void AliGenPythia::SetNuclei(Int_t a1, Int_t a2) } -void AliGenPythia::MakeHeader() const +void AliGenPythia::MakeHeader() { // Builds the event header, to be called after each event AliGenEventHeader* header = new AliGenPythiaEventHeader("Pythia"); +// +// Event type ((AliGenPythiaEventHeader*) header)->SetProcessType(fPythia->GetMSTI(1)); +// +// Number of trials ((AliGenPythiaEventHeader*) header)->SetTrials(fTrials); +// +// Event Vertex header->SetPrimaryVertex(fEventVertex); +// +// Jets that have triggered + if (fProcess == kPyJets) + { + Int_t ntrig, njet; + Float_t jets[4][10]; + GetJets(5.0, 1, njet, ntrig, jets); + for (Int_t i = 0; i < ntrig; i++) { + ((AliGenPythiaEventHeader*) header)->AddJet(jets[0][i], jets[1][i], jets[2][i], + jets[3][i]); + } + } gAlice->SetGenEventHeader(header); } -Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) const +Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) { // Check the kinematic trigger condition // @@ -732,19 +766,15 @@ Bool_t AliGenPythia::CheckTrigger(TParticle* jet1, TParticle* jet2) const Bool_t triggered = kFALSE; if (fProcess == kPyJets) { - //Check eta range first... - if ( - ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet) && - (phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet)) - - || - - ((eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet) && - (phi[1] < fPhiMaxJet && phi[1] > fPhiMinJet)) - ) - { - triggered = kTRUE; - } + Int_t njets = 0; + Int_t ntrig = 0; + Float_t jets[4][10]; +// +// Use Pythia clustering on parton level to determine jet axis +// + GetJets(5.0, 1, njets, ntrig, jets); + if (ntrig) triggered = kTRUE; +// } else { Int_t ij = 0; Int_t ig = 1; @@ -773,6 +803,90 @@ AliGenPythia& AliGenPythia::operator=(const AliGenPythia& rhs) return *this; } +void AliGenPythia::LoadEvent() +{ +// +// Load event into Pythia Common Block +// + + + Int_t npart = (Int_t) (gAlice->TreeK())->GetEntries(); + (fPythia->GetPyjets())->N = npart; + + for (Int_t part = 0; part < npart; part++) { + TParticle *MPart = gAlice->Particle(part); + Int_t kf = MPart->GetPdgCode(); + Int_t ks = MPart->GetStatusCode(); + Float_t px = MPart->Px(); + Float_t py = MPart->Py(); + Float_t pz = MPart->Pz(); + Float_t e = MPart->Energy(); + Float_t p = TMath::Sqrt(px * px + py * py + pz * pz); + Float_t m = TMath::Sqrt(e * e - p * p); + + + (fPythia->GetPyjets())->P[0][part] = px; + (fPythia->GetPyjets())->P[1][part] = py; + (fPythia->GetPyjets())->P[2][part] = pz; + (fPythia->GetPyjets())->P[3][part] = e; + (fPythia->GetPyjets())->P[4][part] = m; + + (fPythia->GetPyjets())->K[1][part] = kf; + (fPythia->GetPyjets())->K[0][part] = ks; + } +} + + +void AliGenPythia::GetJets(Float_t dist, Int_t part, Int_t& nJets, Int_t& nJetsTrig, Float_t jets[4][10]) +{ +// +// Calls the Pythia clustering algorithm to find jets in the current event +// + Int_t n = fPythia->GetN(); + nJets = 0; + nJetsTrig = 0; +// +// Configure cluster algorithm +// + fPythia->SetPARU(43, dist); + fPythia->SetMSTU(41, part); +// +// Call cluster algorithm +// + fPythia->Pyclus(nJets); +// +// Loading jets from common block +// + Int_t i; + for (i = 0; i < nJets; i++) { + Float_t px = (fPythia->GetPyjets())->P[0][n+i]; + Float_t py = (fPythia->GetPyjets())->P[1][n+i]; + Float_t pz = (fPythia->GetPyjets())->P[2][n+i]; + Float_t e = (fPythia->GetPyjets())->P[3][n+i]; + Float_t pt = TMath::Sqrt(px * px + py * py); + Float_t phi = TMath::ATan2(py,px); + Float_t theta = TMath::ATan2(pt,pz); + Float_t et = e * TMath::Sin(theta); + Float_t eta = -TMath::Log(TMath::Tan(theta / 2.)); + + if ( + eta > fEtaMinJet && eta < fEtaMaxJet && + phi > fPhiMinJet && eta < fPhiMaxJet && + et > fEtMinJet && et < fEtMaxJet + ) + { +// printf("\n*Trigger-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg); + jets[0][nJetsTrig] = px; + jets[1][nJetsTrig] = py; + jets[2][nJetsTrig] = pz; + jets[3][nJetsTrig] = e; + nJetsTrig++; + + } else { +// printf("\n........-Jet #%d: %10.3f %10.3f %10.3f %10.3f \n", i, pt, et, eta, phi * kRaddeg); + } + } +} #ifdef never diff --git a/EVGEN/AliGenPythia.h b/EVGEN/AliGenPythia.h index b717d3c5b78..07ad58392de 100644 --- a/EVGEN/AliGenPythia.h +++ b/EVGEN/AliGenPythia.h @@ -46,10 +46,17 @@ class AliGenPythia : public AliGenMC {fPtHardMin = ptmin; fPtHardMax = ptmax; } virtual void SetYHard(Float_t ymin = -1.e10, Float_t ymax = 1.e10) {fYHardMin = ymin; fYHardMax = ymax; } + // Set initial and final state gluon radiation + virtual void SetGluonRadiation(Int_t iIn, Int_t iFin) + {fGinit = iIn; fGfinal = iFin;} + virtual void SetPtKick(Float_t kt) + {fPtKick = kt;} // set centre of mass energy virtual void SetEnergyCMS(Float_t energy = 5500) {fEnergyCMS = energy;} // treat protons as inside nuclei virtual void SetNuclei(Int_t a1, Int_t a2); + virtual void SetJetEtRange(Float_t etmin = 0., Float_t etmax = 1.e4) + {fEtMinJet = etmin; fEtMaxJet = etmax;} virtual void SetJetEtaRange(Float_t etamin = -20., Float_t etamax = 20.) {fEtaMinJet = etamin; fEtaMaxJet = etamax;} virtual void SetJetPhiRange(Float_t phimin = -180., Float_t phimax = 180.) @@ -78,6 +85,9 @@ class AliGenPythia : public AliGenMC // get cross section of process virtual Float_t GetXsection() const {return fXsection;} + // get triggered jets + void GetJets(Float_t dist, Int_t part, Int_t& njets, Int_t& ntrig, Float_t[4][10]); + void LoadEvent(); // Getters virtual Process_t GetProcess() {return fProcess;} virtual StrucFunc_t GetStrucFunc() {return fStrucFunc;} @@ -86,7 +96,7 @@ class AliGenPythia : public AliGenMC virtual Float_t GetEnergyCMS() {return fEnergyCMS;} virtual void GetNuclei(Int_t& a1, Int_t& a2) {a1 = fNucA1; a2 = fNucA2;} - virtual void GetJetEtaRange(Float_t& etamin, Float_t& etamax) + virtual void GetJetEtRange(Float_t& etamin, Float_t& etamax) {etamin = fEtaMinJet; etamax = fEtaMaxJet;} virtual void GetJetPhiRange(Float_t& phimin, Float_t& phimax) {phimin = fPhiMinJet*180./TMath::Pi(); phimax = fPhiMaxJet*180/TMath::Pi();} @@ -96,7 +106,7 @@ class AliGenPythia : public AliGenMC {phimin = fPhiMinGamma*180./TMath::Pi(); phimax = fPhiMaxGamma*180./TMath::Pi();} // virtual void FinishRun(); - Bool_t CheckTrigger(TParticle* jet1, TParticle* jet2) const; + Bool_t CheckTrigger(TParticle* jet1, TParticle* jet2); // Assignment Operator AliGenPythia & operator=(const AliGenPythia & rhs); @@ -104,7 +114,7 @@ class AliGenPythia : public AliGenMC // adjust the weight from kinematic cuts void AdjustWeights(); Int_t GenerateMB(); - void MakeHeader() const; + void MakeHeader(); TClonesArray* fParticles; //Particle List @@ -127,10 +137,15 @@ class AliGenPythia : public AliGenMC Float_t fYHardMax; //higher y-hard cut Int_t fNucA1; //mass number nucleus side 1 Int_t fNucA2; //mass number nucleus side 2 + Int_t fGinit; //initial state gluon radiation + Int_t fGfinal; //final state gluon radiation + Float_t fPtKick; //Transverse momentum kick Bool_t fFullEvent; //!Write Full event if true AliDecayer *fDecayer; //!Pointer to the decayer instance Int_t fDebugEventFirst; //!First event to debug Int_t fDebugEventLast; //!Last event to debug + Float_t fEtMinJet; //Minimum et of triggered Jet + Float_t fEtMaxJet; //Maximum et of triggered Jet Float_t fEtaMinJet; //Minimum eta of triggered Jet Float_t fEtaMaxJet; //Maximum eta of triggered Jet Float_t fPhiMinJet; //Minimum phi of triggered Jet -- 2.31.1