Re: [ROOT] TPythia6: problem setting parameter

From: Rene Brun (Rene.Brun@cern.ch)
Date: Fri Jul 18 2003 - 16:35:52 MEST


Hi Jennifer,

Could you clarify which version of Pythia you are using?

The recent Root releases, including today release 3.05/06 work only
with Pythia 6.2. You should not use them with Pythia6.1

Rene Brun

Jennifer Klay wrote:
> 
> 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