I have write a file in which i am making collision of particles and observing the number of particles produced. In all the collision events, the total produced particles are about 1 million. I have assigned the value of mass using gRandom->Uniform(0-2000). But I want to relate the mass with probability, like I want that for lower mass value the number of particles produced is more and for higher mass particle the production probability is less and if we fill it in histogram then it must follow something like 1/sqrt(x) type relation. How can i obtain such plot for mass.
Hi @Sachin_Rana,
welcome to the ROOT forum. Could you maybe share the code you wrote already?
Cheers,
Marta
void lorentz_vector()
{
TCanvas *c1 =new TCanvas();
c1->Divide(2,2);
TH1F *histm = new TH1F("histm" , "Mass" ,500, 0,2000);
TH1F *histpx = new TH1F("histpx" , "Px" ,500, 0,12);
TH1F *histpy = new TH1F("histpy" , "Py" ,500, 0,12 );
TH1F *histpz = new TH1F("histpz" , "Pz" ,500, 0,12 );
TRandom2 *random= new TRandom2(1);
double_t px, py, pz,p, trans, e, m ;
for (int i=1; i<=10000; i++){
m = gRandom->Uniform(0,2000); // use this for random gaussian value
// cout<<m<<endl;
TF1 *f1 = new TF1("f1","[0]*TMath::Power(([1]/[2]),(x/[2]))*(TMath::Exp(-([1]/[2])))/TMath::Gamma((x/[2])+1)",.02,12);// use this for user defined random poisson value ;
f1->SetParameters(1,1,1);
px = f1->GetRandom();
py = f1->GetRandom();
pz = f1->GetRandom();
trans = sqrt(px*px+py*py);
p=sqrt(px*px+py*py+pz*pz);
e = sqrt(px*px+py*py+pz*pz+m*m);
histm->Fill(m);
histpx->Fill(px);
histpy->Fill(py);
histpz->Fill(pz);
// cout<<" Mass: "<< m<<" Px: "<<px<<" Py: " << py<<" Pz: " << pz<<" Pt: " << trans <<" Energy: " << e <<endl;
}
TLorentzVector v;
v.SetPxPyPzE(px,py,pz,e);
cout<<"X-component of momentum using lorentz vector= "<<v.Px()<<endl;
cout<< "Y-component of momentum using lorentz vector= " <<v.Py()<< endl;
cout<<"Z-component of momentum using lorentz vector= " <<v.Pz()<<endl;
cout << "Energy using lorentz_vector= " <<v.E()<<endl;
double_t mag, theta, phi, pt, eta;
mag= v.Rho();
theta=v.Theta()*180/TMath::Pi(); // we have multiply this with 180/pi to convert it to degree
phi = v.Phi();
pt = v.Perp();
eta = v.PseudoRapidity();
cout <<" Magnitude of selected lorentz vector: "<< mag << " THETA: " << theta << " PHI:" << phi << " Pt By Lorentz: "<< pt << " Pseudorapidity: " << eta<< endl;
double pseudo,energy;
pseudo=0.5*log((p+pz)/(p-pz));
energy=sqrt(p*p+m*m);
cout<<"pseudorapidity using p and pz is="<<pseudo<<endl;
cout<<"energy using formula is="<<energy<<endl;
cout<<"transverse momentum using formula is="<<trans<<endl;
double pseudorapi=-2.303*0.5*log(tan(theta/2));
cout<<"psedorapidity using theta ="<<pseudorapi<<endl;
double pxx=trans*cos(phi);
cout<<"pxx=pt*cos(phi)="<<pxx<<endl;
//py=pt*sin(phi)
//pz=pt*sinh(eta)
c1->cd(1);
histm->Draw();
c1->cd(2);
histpx->Draw();
c1->cd(3);
histpy->Draw();
c1->cd(4);
histpz->Draw();
}
I am unable to share the link so i am pasting the code here.
Hi @Sachin_Rana,
thank you for the code. Iām still not sure if I fully understand what you want to do, but basically instead of m = gRandom->Uniform(0,2000); you want to have the mass distributed as 1/sqrt(m)? In that case you can first define your custom function (in same way as you do it for Poisson, which is also an already defined function in ROOT so you can just use that instead) and then randomly sample from it using the GetRandom?
And here is more information on TRandom class: ROOT: TRandom Class Reference
Cheers,
Marta
Like what i am trying to say is that lets say total number of produced particles are 1 million and in these produced particles there are 50 types of particles(for example lets say out of 1 million, 100k are electrons, 100k are positrons, 5k are muon and so on), so we assign particle id to each type of particle. Now production of particles is in such a way that for any two types of particle the probability of their production will follow the following relation.
Lets say i am relating probability for muon and pion:
P(muon)=x*P(pion) : Where x =(mass of pion)/(mass of muon)
So in this way i want to fill histogram for mass.
Hi,
this is a small program doing what you are asking for
void ParticleGen_Mass()
{
const int n_particles=4;
double mass[n_particles]={105.7,139.6,938.3,1115.7};//mu+ mass, pi+ mass, proton, lambda mass
double mass_ratio[n_particles];
TH1D *h =new TH1D("h","mass generated",2300,50,1200.);
TRandom3 *r = new TRandom3();
for(int i=0;i<n_particles;i++)
mass_ratio[i]=mass[0]/mass[i];
int gen_particle=0;
int idx_particle;
while(gen_particle<1e4)
{
idx_particle = floor(r->Uniform(n_particles));
if( r->Uniform(1) < mass_ratio[idx_particle])
{
h->Fill(mass[idx_particle]);
gen_particle++;
}
}
h->Draw();
}
I filled the histogram of mass, you generate the particles randomly with the probability distribution you asked for.
The only important thing is that the particle with the smallest mass is the first element of the mass vector.
Best,
Stefano
1 Like
This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.