Hi,
Here is the macro I used to run Pythia. It is essentially an adaptation of
the example available in $ROOTSYS/tutorials.
After I generate the events, I plot the status code for pi-zeros in the
TreeViewer and they all have status 11.
Thanks,
Jennifer
==================================================
#include <cstdlib>
#include <iostream>
using namespace std;
#define FILENAME "pythia2.root"
#define TREENAME "tree"
#define BRANCHNAME "particles"
#define HISTNAME "ptSpectra"
#define PDGNUMBER 211
//----------------------------------------------------------------------
void RunPythia(Int_t n=100) {
generateEvents(n);
showEventSample();
}
//----------------------------------------------------------------------
// nEvents is how many events we want.
int generateEvents(Int_t nEvents)
{
//----------------------
// Load needed libraries
//----------------------
loadLibraries();
//-----------------------------------------------------
// Create an instance of the Pythia event generator ...
//-----------------------------------------------------
TPythia6* pythia = new TPythia6;
//--------------------------------------------------
// Initialise it to run p+p at sqrt(5500) GeV in CMS
//--------------------------------------------------
pythia->Initialize("cms", "p", "p", 5500);
//-----------------------------------------------
// Set up all other special codes/settings here:
//-----------------------------------------------
// pythia->SetCKIN(3,10.); //Minimum Jet Et = 10 GeV
// pythia->SetCKIN(7,-1.0); //Minimum rapidity
// pythia->SetCKIN(8,1.0); //Maximum rapidity
// pythia->SetMSTP(111,1); //Fragmentation on (0 = off)
pythia->SetMDCY(111,1,0); //Don't decay pizeros
//----------------------
// Open an output file
//----------------------
TFile* file = TFile::Open(FILENAME, "RECREATE");
if (!file || !file->IsOpen()) {
Error("generateEvents", "Couldn;t open file %s", FILENAME);
return 1;
}
//--------------------------------
// Make a tree in that file ...
//--------------------------------
TTree* tree = new TTree(TREENAME, "Pythia 6 tree");
//------------------------------------------------------------------------------------
// Register the cache of pythia on a branch (TClonesArray of TMCParticle
objects. )
//------------------------------------------------------------------------------------
TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
tree->Branch(BRANCHNAME, &particles);
//---------------------------
// Now we make some events
//---------------------------
for (Int_t i = 0; i < nEvents; i++) {
// Show how far we got every 100'th event.
if (i % 100 == 0)
cout << "Event # " << i << endl;
//-----------------
// Make one event.
//-----------------
pythia->GenerateEvent();
// Maybe you want to have another branch with global event
// information. In that case, you should process that here.
// You can also filter out particles here if you want.
//---------------------------------------------------------
// Now we're ready to fill the tree, and the event is over.
//---------------------------------------------------------
tree->Fill();
}
//--------------------
// Show tree structure
//--------------------
tree->Print();
//--------------------------
// Flush and close the file
//--------------------------
file->Write();
file->Close();
return 0;
}
//----------------------------------------------------
// Start the tree viewer.
//----------------------------------------------------
int showEventSample()
{
// Load needed libraries
loadLibraries();
// Open the file
TFile* file = TFile::Open(FILENAME, "READ");
if (!file || !file->IsOpen()) {
Error("showEventSample", "Couldn;t open file %s", FILENAME);
return 1;
}
// Get the tree
TTree* tree = (TTree*)file->Get(TREENAME);
if (!tree) {
Error("showEventSample", "couldn't get TTree %s", TREENAME);
return 2;
}
TClonesArray* particles;
tree->SetBranchAddress(BRANCHNAME,&particles);
Int_t nEvt = (Int_t)tree->GetEntries();
cout << "There are " << nEvt << " events" << endl;
// tree->Print();
// Start the viewer.
tree->StartViewer();
return 0;
}
//-------------------------------------------------------------
void loadLibraries()
{
// Load the Event Generator abstraction library, Pythia 6
// library, and the Pythia 6 interface library.
gSystem->Load("libEG");
gSystem->Load("libEGPythia6");
gSystem->Load("libPythia6");
}
==============================================================
==============================
Jennifer L. Klay
Lawrence Postdoctoral Fellow
Relativistic Nuclear Collisions Group
Nuclear Science Division
Lawrence Berkeley National Laboratory
==============================
Rene Brun wrote:
> Hi Jennifer,
>
> I do not see any problem with TPythia6::SetMDCY.
>
> Did you initialize the class?
> Could you send a short example showing the problem?
>
> Rene Brun
>
> Jennifer Klay wrote:
> >
> > Hi,
> > I am using the TPythia6 libraries with ROOT to do some simulations on a
> > linux machine running Redhat linux. I have tried running in several
> > different versions of ROOT: 3.02.07, 3.03.09, 3.05.04.
> >
> > Using the $ROOTSYS/tutorials/pythiaExample.C, I tried to set the
> > parameter that turns OFF neutral pion decay:
> >
> > pythia->SetMDCY(111,1,0);
> >
> > Setting or unsetting this has no effect. Neutral pions are always
> > decayed (have status code 11). I checked the value of the parameter at
> > the end of the macro by using the Get method:
> > cout << pythia->GetMDCY(111,1) << endl;
> >
> > and it verifies that the value set is zero.
> >
> > When running PYTHIA as a standalone fortran program, this flag works, so
> > I suspect that the value is not being properly passed through the
> > wrapper code to the fortran implementation.
> >
> > Are there any workarounds to set this parameter and still use the ROOT
> > classes? I want to generate TTrees of the event record so I can run
> > some analysis after generating the events.
> >
> > Thanks,
> > Jennifer
> >
> > ==============================
> > Jennifer L. Klay
> > Lawrence Postdoctoral Fellow
> > Relativistic Nuclear Collisions Group
> > Nuclear Science Division
> > Lawrence Berkeley National Laboratory
> > ==============================
This archive was generated by hypermail 2b29 : Thu Jan 01 2004 - 17:50:13 MET