Hi,
I’m trying to plot transverse momentum in pythia6 with ROOT, but I have a problem, please someone could help me. This is the macro that I use in root
[code]#ifndef CINT
#include “TApplication.h”
#include “TPythia6.h”
#include “TFile.h”
#include “TError.h”
#include “TTree.h”
#include “TClonesArray.h”
#include “TH1.h”
#include “TF1.h”
#include “TStyle.h”
#include “TLatex.h”
#include “TCanvas.h”
#include “Riostream.h”
#include
using namespace std;
#endif
// This funtion just loads the needed libraries if we’re executing from
// an interactive session.
void loadLibraries()
{
#ifdef CINT
// Load the Event Generator abstraction library, Pythia 6
// library, and the Pythia 6 interface library.
gSystem->Load(“libEG”);
gSystem->Load(“usr/local/lib/root/libPythia6”); //change to your setup
gSystem->Load(“libEGPythia6”);
#endif
}
float m;
int ks, kf;
//int makeEventSample(Int_t nEvents)
// This is the main part
//
void pythia_E()
{
// Load needed libraries
loadLibraries();
// Create an instance of the Pythia event generator …
TPythia6* mt = new TPythia6;
// Select production
mt->SetMSEL(0);
mt->SetMSUB(271,1); //process 271
mt->SetMSUB(272,1); //process 272
mt->SetIMSS(1,1); // MSSM simulation
mt->SetMSTP(81,0); //no multiples interacctions
// … and initialise it to run p+p at 10000 GeV in CMS
mt->Initialize(“cms”, “p”, “p”, 14000.);
// (re)creating histogram file
TFile* file = new TFile(“pythia_E.root”, “recreate”);
// booking histogram
TH1F* energy = new TH1F(“momentum”, "Transverse Momentum for Squarks “,100,0.,1000.);
// Now we make some events
for (Int_t i = 0; i < 1000; i++)
{
// Make an event.
mt->GenerateEvent();
// Print the first event
if(i==0) mt->Pylist(1);
// loop over all particles in the event
for (Int_t part=1; part<=mt->GetN(); part++)
{
// select only final state particles
if (ks = mt->GetK(part,1) == 21)
{
// flavor code (~u_L= 1000002)
if(kf = mt->GetK(part,2) == 1000002)
{
// transverse momentum
m = mt->GetP(part,10);
// fill histogram
energy->Fill(m,1.);
}
}
}
}
// print statistics
mt->Pystat(1);
// open graphics window
c1 = new TCanvas(“pythia_E”,” ",800,600);
// plot energy distribution - histogram
energy->Draw();
// save histogram file
energy->Write();
}
[/code]
I guess I have problem withGetK or GetP.
Anybody who knows how to solve this problem?
Thanks.
LINA[/i]