Pythia Histogram Problem

Dear Rooters,

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.

Any suggestions or ideas would be appreciated.

Thank you.
pythiatest2.C (1.81 KB)

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.

Rene
pythiatest2.C (1.89 KB)

Dear Rene,

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.

Kind regards.

well! it is not always 23. In 10 events you have two cases at 25.
It depends on the selected process.

Rene

Dear Rene,

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.

You cannot print charge outside the loop. The variable is only defined inside the loop.

Rene

But if I write

printf(“charge %d\n”,charge);

after the line

if (finalstate != 1) continue;

Macro does not give any charge value. But I can see the result histogram.

My aim is to plot final state charged particle histogram. Do you think my tree is fill with these type particles? or should I change the location of

part->Fill();

Thank you very much.

see me previous comments above and the code in attachment

Rene
pythiatest2.C (2.03 KB)

Dear Rene,

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

part->Fill();

line to above

}

Could you please help me to solve this problem.

Thank you.

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).

Cheers,
Philippe.[/code]

Could you please tell me to get only charged particles should I move

part->Fill();

line inside of the upper paranthesis like that,
part->Fill();

}
}
part->Write();

Thank you.

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 }

Cheers,
Philippe.

Dear Philippe,

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.

Thank you very much.

Hi,

What did you try? You may want to re-read the user’s guide chapter on collection.

Cheers,
Philippe.

Dear Philippe,

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?

I think something is wrong, but I don’t know what :frowning:

Thank you.

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.

Cheers,
Philippe.

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);

 Int_t fKF = 0, Int_t fKS = 0 ;
 fKF=prt->GetKF();
 fKS=prt->GetKS();
 TParticlePDG *pdg = TDatabasePDG::Instance()->GetParticle(fKS);
 
 Int_t charge;
 charge  = TDatabasePDG::Instance()->GetParticle(fKF)->Charge();
 
 if (charge == 0) continue; 
if(fKS!=1) 
 {prt->TClonesArray::RemoveAt(j)}
printf("charge=%d\n, fKS=%d\n",charge,fKS);       
                                                                      }

part->Fill();
}

At this stage I see right charge values but fKS values are not always 1. Could you please tell me what should I do?

Thank you.

Could you please let me know how can I remove an element from the array with using TClonesArray::RemoveAt()

if (charge==0)
{
// remove the element from the array
}

I try to use this my code but it does not work. (Please see my above message)

I think something is wrong but I could not find what it is?

Thank you very much,

Kind regards.

Hi,

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);

for (Int_t i = 0; i < nEvents; i++) {

  pythia->GenerateEvent();

  Int_t nsel = 0;
  selected->Clear();
  for (Int_t j=0; j<particles->GetEntries(); j++) {
     prt=(TMCParticle*)particles->At(j);

     Int_t fKF = 0, Int_t fKS = 0 ;
     fKF=prt->GetKF();
     fKS=prt->GetKS();
     TParticlePDG *pdg = TDatabasePDG::Instance()->GetParticle(fKS);

     Int_t charge;
     charge = TDatabasePDG::Instance()->GetParticle(fKF)->Charge();

     if (charge != 0 && fKS==1) {
        new ((*selected)[nsel++]) TMCParticle( *prt );
        printf("accepted charge=%d\n, fKS=%d\n",charge,fKS);
     } else {
         printf("rejected charge=%d\n, fKS=%d\n",charge,fKS);
     }
  }

part->Fill();
}
[/code]

Cheers,
Philippe

*** Note that I added a missing call to Clear in my code fragment