/*
$Log$
+Revision 1.68 2002/12/11 09:16:16 morsch
+Use GetJets to fill header.
+
+Revision 1.67 2002/12/09 15:24:09 morsch
+Same trigger routine can use Pycell or Pyclus.
+
+Revision 1.66 2002/12/09 08:22:56 morsch
+UA1 jet finder (Pycell) for software triggering added.
+
+Revision 1.65 2002/11/19 08:57:10 morsch
+Configuration of pt-kick added.
+
+Revision 1.64 2002/11/15 00:43:06 morsch
+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.
+
+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
+
+Revision 1.62 2002/09/24 10:00:01 morsch
+CheckTrigger() corrected.
+
+Revision 1.61 2002/07/30 09:52:38 morsch
+Call SetGammaPhiRange() and SetGammaEtaRange() in the constructor.
+
+Revision 1.60 2002/07/19 14:49:03 morsch
+Typo corrected.
+
+Revision 1.59 2002/07/19 14:35:36 morsch
+Count total number of trials. Print mean Q, x1, x2.
+
+Revision 1.58 2002/07/17 10:04:09 morsch
+SetYHard method added.
+
Revision 1.57 2002/05/22 13:22:53 morsch
Process kPyMbNonDiffr added.
// andreas.morsch@cern.ch
//
+#include <TDatabasePDG.h>
+#include <TParticle.h>
+#include <TPDGCode.h>
+#include <TSystem.h>
+#include <TTree.h>
+#include "AliConst.h"
+#include "AliDecayerPythia.h"
#include "AliGenPythia.h"
#include "AliGenPythiaEventHeader.h"
-#include "AliDecayerPythia.h"
-#include "AliRun.h"
#include "AliPythia.h"
-#include "AliPDG.h"
-#include <TParticle.h>
-#include <TSystem.h>
+#include "AliRun.h"
ClassImp(AliGenPythia)
SetEventListRange();
SetJetPhiRange();
SetJetEtaRange();
+ SetJetEtRange();
+ SetGammaPhiRange();
+ SetGammaEtaRange();
+ SetPtKick();
}
AliGenPythia::AliGenPythia(Int_t npart)
SetStrucFunc();
SetForceDecay();
SetPtHard();
+ SetYHard();
SetEnergyCMS();
fDecayer = new AliDecayerPythia();
// Set random number generator
SetEventListRange();
SetJetPhiRange();
SetJetEtaRange();
+ SetJetEtRange();
+ SetGammaPhiRange();
+ SetGammaEtaRange();
+ SetJetReconstructionMode();
+ SetPtKick();
// Options determining what to keep in the stack (Heavy flavour generation)
fStackFillOpt = kFlavorSelection; // Keep particle with selected flavor
fFeedDownOpt = kTRUE; // allow feed down from higher family
void AliGenPythia::Init()
{
// Initialisation
- SetMC(AliPythia::Instance());
+
+ SetMC(AliPythia::Instance());
fPythia=(AliPythia*) fgMCEvGen;
//
fParentWeight=1./Float_t(fNpart);
fPythia->SetMSTP(111,0);
}
+
+// initial state radiation
+ fPythia->SetMSTP(61,fGinit);
+// final state radiation
+ fPythia->SetMSTP(71,fGfinal);
+// pt - kick
+ if (fPtKick > 0.) {
+ fPythia->SetMSTP(91,1);
+ fPythia->SetPARP(91,fPtKick);
+ } else {
+ fPythia->SetMSTP(91,0);
+ }
+
+ // fPythia->SetMSTJ(1,2);
+ //
fPythia->ProcInit(fProcess,fEnergyCMS,fStrucFunc);
- // fPythia->Pylist(0);
- // fPythia->Pystat(2);
// Parent and Children Selection
switch (fProcess)
{
case kPyDirectGamma:
break;
}
+//
+// This counts the total number of calls to Pyevnt() per run.
+ fTrialsRun = 0;
+ fQ = 0.;
+ fX1 = 0.;
+ fX2 = 0.;
+ fNev = 0 ;
+//
AliGenMC::Init();
}
pSelected[i] = 0;
trackIt[i] = 0;
}
- // printf("\n **************************************************%d\n",np);
+
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 && fProcess != kPyDirectGamma &&
+ if (fProcess != kPyMb && fProcess != kPyJets &&
+ fProcess != kPyDirectGamma &&
fProcess != kPyMbNonDiffr) {
for (i = 0; i<np; i++) {
TParticle * mother = (TParticle *) fParticles->At(ipa);
kfMo = TMath::Abs(mother->GetPdgCode());
}
-// printf("\n particle (all) %d %d %d", i, pSelected[i], kf);
// What to keep in Stack?
Bool_t flavorOK = kFALSE;
Bool_t selectOK = kFALSE;
if (jev >= fNpart || fNpart == -1) {
fKineBias=Float_t(fNpart)/Float_t(fTrials);
printf("\n Trials: %i %i %i\n",fTrials, fNpart, jev);
+
+ fQ += fPythia->GetVINT(51);
+ fX1 += fPythia->GetVINT(41);
+ fX2 += fPythia->GetVINT(42);
+ fTrialsRun += fTrials;
+ fNev++;
MakeHeader();
break;
}
kf = CheckPDGCode(iparticle->GetPdgCode());
Int_t ks = iparticle->GetStatusCode();
Int_t km = iparticle->GetFirstMother();
-// printf("\n Particle: %d %d %d", i, kf, ks);
-
if ((ks == 1 && kf!=0 && KinematicSelection(iparticle, 0)) ||
(ks != 1) ||
(fProcess == kPyJets && ks == 21 && km == 0 && i>1)) {
{
// Print x-section summary
fPythia->Pystat(1);
+ fQ /= fNev;
+ fX1 /= fNev;
+ fX2 /= fNev;
+ printf("\nTotal number of Pyevnt() calls %d\n", fTrialsRun);
+ printf("\nMean Q, x1, x2: %f %f %f\n", fQ, fX1, fX2);
+
+
}
void AliGenPythia::AdjustWeights()
}
-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(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
//
Bool_t triggered = kFALSE;
if (fProcess == kPyJets) {
- //Check eta range first...
- if ((eta[0] < fEtaMaxJet && eta[0] > fEtaMinJet) ||
- (eta[1] < fEtaMaxJet && eta[1] > fEtaMinJet))
- {
- //Eta is okay, now check phi range
- if ((phi[0] < fPhiMaxJet && phi[0] > fPhiMinJet) ||
- (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(njets, ntrig, jets);
+
+ if (ntrig) triggered = kTRUE;
+//
} else {
Int_t ij = 0;
Int_t ig = 1;
}
}
}
-
-
return triggered;
}
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::RecJetsUA1(Float_t eCellMin, Float_t eCellSeed, Float_t eMin, Float_t rMax,
+ Int_t& njets, Float_t jets [4][50])
+{
+//
+// Calls the Pythia jet finding algorithm to find jets in the current event
+//
+//
+// Configure detector (EMCAL like)
+//
+ fPythia->SetPARU(51,2.);
+ fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
+ fPythia->SetMSTU(52,3 * 144);
+//
+// Configure Jet Finder
+//
+ fPythia->SetPARU(58, eCellMin);
+ fPythia->SetPARU(52, eCellSeed);
+ fPythia->SetPARU(53, eMin);
+ fPythia->SetPARU(54, rMax);
+ fPythia->SetMSTU(54, 2);
+//
+// Save jets
+ Int_t n = fPythia->GetN();
+
+//
+// Run Jet Finder
+ fPythia->Pycell(njets);
+ 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];
+
+ jets[0][i] = px;
+ jets[1][i] = py;
+ jets[2][i] = pz;
+ jets[3][i] = e;
+ }
+}
+
+
+
+void AliGenPythia::GetJets(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;
+ if (fJetReconstruction == kCluster) {
+//
+// Configure cluster algorithm
+//
+ fPythia->SetPARU(43, 2.);
+ fPythia->SetMSTU(41, 1);
+//
+// Call cluster algorithm
+//
+ fPythia->Pyclus(nJets);
+//
+// Loading jets from common block
+//
+ } else {
+//
+// Configure detector (EMCAL like)
+//
+ fPythia->SetPARU(51,2.);
+ fPythia->SetMSTU(51,Int_t(96 * 2./0.7));
+ fPythia->SetMSTU(52,3 * 144);
+//
+// Configure Jet Finder
+//
+ fPythia->SetPARU(58, 0.0);
+ fPythia->SetPARU(52, 4.0);
+ fPythia->SetPARU(53, 10.0);
+ fPythia->SetPARU(54, 1.0);
+ fPythia->SetMSTU(54, 2);
+//
+// Run Jet Finder
+ fPythia->Pycell(nJets);
+ }
+
+ 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
+ )
+ {
+ 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
git://git.uio.no / u / mrichter / AliRoot.git / blobdiff