I would like to plot d_Sigma/d_Pt histogram for final state charged particles with using PYTHIA. I used a small macro which I found in this forum. My problem is when I run the macro I get an empty histogram. I have attached the macro for you. Could you please help me to solve this problem.
We provide a standard tutorial (see $ROOTSYS/tutorials/pythia/pythiaExample.C)
Your code was apparently extracted from this example, but with many mistakes.
See correct code in attachment.
Thank you for your help. I am working on this problem for two days.Finally I find the solution.
I would like to ask you what is the reason of prt always goes to 23. Is it normal? I think it should be number of entries in the particles branch. So, it should not change sequencely.
I have one more problem. If I write
printf(" charge %d\n", charge);
between the lines
}
and
part->Fill();
I could not get only 0 charges. Do you think my tree is fill with only charged final state particles? Could you please let me know what is the reason of this different charge values.
Thank you for this new macro. Yes, I can see the charge values now.
When I write
}
part->Fill();
}
line at the end of the first loop. I think that (but I am not sure) whole particles that fill in a tree have some 0 charges. I think tree is filled with all particles. But I just want charged particles.I am not sure should I move
if (charge == 0) continue;The keyword ‘continue’ request that the code terminates the current iteration of the loop and proceeds directly to the next. In the case of a for loop it jumps to its increment expression. (Aka if charge==0, the code after the continue will not be executed).
Nope nothing like that.
You have an array of particle that you associated with the TTree. The array of particle will always be fully stored if you call part->Fill.
Instead of ‘if (charge==0) continue;’ you should do if (charge==0)
{
// remove the element from the array
}
Thank you for your reply. I have some hesitation on this code working properly. Now it is clear. But I am a new root user, so I don’t know how can I remove an element from array. Could you please give me an example. Or do you know any example macro to do this.
I just want to plot final state charged particle histogram. So , you wrote me that I should remove unwanted element from my array and then fill the tree. I have added
if (charge==0){prt->Delete()}
line in my code. But I have noticed that if I remove
pythia->SetupTest();
line , my charge is showing 0 values. Even if I when I turn off
if (charge==0){prt->Delete()}
line I still get values hich does not have 0 (I get -3 and +3 values) . It is so strange for me. Could you explain this?
if (charge==0){prt->Delete()} This only destroy the content of the object but does NOT remove it from the collection. Try TClonesArray::RemoveAt instead.
Dear Philippe,
Thank you for the hint. My final code is like that,
int makeEventSample(Int_t nEvents) {
loadLibraries();
TPythia6* pythia=TPythia6::Instance();
pythia->Initialize(“CMS”,“p”,“p”, 900);
const int nEvents=1000;
TTree* part = new TTree(“part”,“A tree of particles”);
TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
part->Branch(“particles”, &particles);
for (Int_t i = 0; i < nEvents; i++) {
pythia->GenerateEvent();
for (Int_t j=0; jGetEntries(); j++) {
prt=(TMCParticle*)particles->At(j);
I think you meant to write something like: if (charge == 0 || fKS!=1)
{
prt->TClonesArray::RemoveAt(j)}
printf("charge=%d\n, fKS=%d\n",charge,fKS);
}
Anyway rather than remove element, it is better for TClonesArray to create a new one with only the value you need:[code]int makeEventSample(Int_t nEvents) {
loadLibraries();
TPythia6* pythia=TPythia6::Instance();
pythia->Initialize(“CMS”,“p”,“p”, 900);
const int nEvents=1000;
TTree* part = new TTree(“part”,“A tree of particles”);
TClonesArray* particles = (TClonesArray*)pythia->GetListOfParticles();
TClonesArray* selected = new TClonesArray(“TMCParticle”);
part->Branch(“particles”, &particles);