Using TGenPhaseSpace I decay 4 particles. However, particle 3 (line 235) produces a triangular curve as opposed to one that is more Gaussian in nature. It should not do this as this creates two sets of two identical particles- therefore there should be two sets of two identical results. The spiking in the graphs is the case for whichever particle is called in the third position (and the order they are referenced does not change the physics). I have found one thing in particular that corrects this- changing the second value on line 43 from 4 to 11. I do not know why this works as this is the number of decay products, so should work when set to 4. I also found listing particles with like mass together reduced this issue slightly, but did not solve it. Does anyone know why changing these values fix es the issue, or more importantly, why they become spiked at all?
–After further investigation, the final particle called, whether there are 3 or 4 particles and regardless of which particle it is, always comes out as a triangular, rather than Gaussian shape.
Thanks in advance!
This is the code producing ‘spiked peaks’:
#include “TGenPhaseSpace.h”
#include “TLorentzVector.h”
#include “TROOT.h”
#include “TStyle.h”
#include “TCanvas.h”
#include “TSystem.h”
#include “TH2.h”
#include “TH1.h”
#include
using namespace std;
void test7_1(UInt_t NrOfEvents = 1000000) {
if (!gROOT->GetClass(“TGenPhaseSpace”)) gSystem->Load(“libPhysics”);
TLorentzVector target(0.0, 0.0, 0.0, 1.876);
TLorentzVector beam(0.0, 0.0, 2.2, 2.2);
TLorentzVector W = beam + target;
// cout << “W=”<< W << endl;
//(Momentum, Energy units are Gev/C, GeV)
Double_t masses[4] = {0.938, 0.938, 0.493, 0.493} ;
/* int c=0;
TRandom3 p(-1);
for (UInt_t iEvent = 0; iEvent < NrOfEvents; ++iEvent) {
Double_t num = p.Uniform(0,360);
if (num > 0.05){
c=c+1;
}
}
cout << “c=”<< c << endl;*/
double phiMin = 10.0; // degrees
gStyle->SetOptStat(0);
//int detectn = 10000; //detect n will eventually be an output based on detector efficiency and angle p
TGenPhaseSpace event;
event.SetDecay(W, 4, masses);
TCanvas *c1 = new TCanvas(“c1”,“c1”, 50, 50, 2000, 1200);
c1->Divide(4,5);
TCanvas *c2 = new TCanvas(“c2”,“c2”, 50, 50, 2000, 1200);
c2->Divide(2,2);
TCanvas *c3 = new TCanvas(“c3”,“c3”, 50, 50, 2000, 1200);
c3->Divide(4,5);
TH2F *dist = new TH2F(“dist”,“dist”, 100,0,6, 100,0,6);
dist->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
dist->GetYaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h1 = new TH1F(“h1”,“kminus mass^2”, 100,0,3);
h1->SetFillColor(kRed);
h1->GetXaxis()->SetTitle(“Mass^2 [GeV/C^2]”);
h1->GetYaxis()->SetTitle(“Counts”);
TH1F *h3 = new TH1F(“h3”,“Proton mass^2”, 100,0,3);
h3->SetFillColor(kRed);
h3->GetXaxis()->SetTitle(“Mass^2 [GeV/C^2]”);
TH1F *h4 = new TH1F(“h4”,“Neutron mass^2”, 100,0,3);
h4->SetFillColor(kRed);
h4->GetXaxis()->SetTitle(“Mass^2 [GeV/C^2]”);
TH1F *h5 = new TH1F(“h5”,“Kplus Mass^2”, 100,0,3);
h5->SetFillColor(kRed);
h5->GetXaxis()->SetTitle(“Mass^2 [GeV/C^2]”);
TH1F *h6 = new TH1F(“h6”,“Rho Km”,100,0,3);
h6->GetXaxis()->SetTitle(“Radial Distance, Rho [rads]”);
TH1F *h7 = new TH1F(“h7”,“Rho Kp”,100,0,3);
h7->GetXaxis()->SetTitle(“Radial Distance, Rho [rads]”);
TH1F *h8 = new TH1F(“h8”,“Rho Proton”,100,0,3);
h8->GetXaxis()->SetTitle(“Radial Distance, Rho [rads]”);
TH1F *h9 = new TH1F(“h9”,“Rho Neutron”,100,0,3);
h9->GetXaxis()->SetTitle(“Radial Distance, Rho [rads]”);
TH1F *h10 = new TH1F(“h10”,“Phi Km”,100,-180,180);
h10->GetXaxis()->SetTitle(“Azimuthal Angle, Phi [degs]”);
TH1F *h11 = new TH1F(“h11”,“Phi Kp”,100,-180,180);
h11->GetXaxis()->SetTitle(“Azimuthal Angle, Phi [degs]”);
TH1F *h12 = new TH1F(“h12”, “Phi P”,100,-180,180);
h12->GetXaxis()->SetTitle(“Azimuthal Angle, Phi [degs]”);
TH1F *h13 = new TH1F(“h13”, “Phi N”,100,-180,180);
h13->GetXaxis()->SetTitle(“Azimuthal Angle, Phi [degs]”);
TH1F *h14 = new TH1F(“h14”, “Km Theta”,100,-10,180);
h14->SetFillColor(kRed);
h14->GetXaxis()->SetTitle(“Polar Angle, Theta [degs]”);
TH1F *h15 = new TH1F(“h15”, “Kp Theta”,100,-10,180);
h15->SetFillColor(kRed);
h15->GetXaxis()->SetTitle(“Polar Angle, Theta [degs]”);
TH1F *h16 = new TH1F(“h16”, “Proton Theta”,100,-10,180);
h16->SetFillColor(kRed);
h16->GetXaxis()->SetTitle(“Polar Angle, Theta [degs]”);
TH1F *h17 = new TH1F(“h17”, “Neutron Theta”,100,-10,180);
h17->SetFillColor(kRed);
h17->GetXaxis()->SetTitle(“Polar Angle, Theta [degs]”);
TH2F *h18 = new TH2F(“h18”, “Kp + Km m^2”,100,0.2,0.3, 100,0.2,0.3);
h18->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
h18->GetYaxis()->SetTitle(“Momentum[GeV/c]”);
TH2F *h19 = new TH2F(“h19”, “P + N m^2”,100,0.8,1, 100,0.8,1);
h19->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
h19->GetYaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h20 = new TH1F(“h20”, “Km Px”,100,-2,2);
h20->SetFillColor(kRed);
h20->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h21 = new TH1F(“h21”, “Kp Px”,100,-2,2);
h21->SetFillColor(kRed);
h21->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h22 = new TH1F(“h22”, “P Px”,100,-2,2);
h22->SetFillColor(kRed);
h22->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h23 = new TH1F(“h23”, “N Px”,100,-2,2);
h23->SetFillColor(kRed);
h23->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h24 = new TH1F(“h24”, “Km Py”,100,-2,2);
h24->SetFillColor(kRed);
h24->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h25 = new TH1F(“h25”, “Kp Py”,100,-2,2);
h25->SetFillColor(kRed);
h25->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h26 = new TH1F(“h26”, “P Py”,100,-2,2);
h26->SetFillColor(kRed);
h26->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h27 = new TH1F(“h27”, “N Py”,100,-2,2);
h27->SetFillColor(kRed);
h27->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h28 = new TH1F(“h28”, “Km Pz”,100,-2,2);
h28->SetFillColor(kRed);
h28->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h29 = new TH1F(“h29”, “Kp Pz”,100,-2,2);
h29->SetFillColor(kRed);
h29->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h30 = new TH1F(“h30”, “P Pz”,100,-2,2);
h30->SetFillColor(kRed);
h30->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *h31 = new TH1F(“h31”, “N Pz”,100,-2,2);
h31->SetFillColor(kRed);
h31->GetXaxis()->SetTitle(“Momentum[GeV/c]”);
TH1F *kmm2_cut = new TH1F(“kmm2_cut”,“K minus mass^2”, 100,0,3);
kmm2_cut->SetFillColor(kGreen);
TH1F *pm2_cut = new TH1F(“pm2_cut”,“Proton mass^2”, 100,0,3);
pm2_cut->SetFillColor(kGreen);
TH1F *nm2_cut = new TH1F(“nm2_cut”,“Neutron mass^2”, 100,0,3);
nm2_cut->SetFillColor(kGreen);
TH1F *kpm2_cut = new TH1F(“kpm2_cut”,“K plus mass^2”, 100,0,3);
kpm2_cut->SetFillColor(kGreen);
TH1F *kmt_cut = new TH1F(“kmt_cut”,“Km Theta”, 100,-10,180);
kmt_cut->SetFillColor(kGreen);
TH1F *kpt_cut = new TH1F(“kpt_cut”,“Kp Theta”, 100,-10,180);
kpt_cut->SetFillColor(kGreen);
TH1F *pt_cut = new TH1F(“pt_cut”,“Proton Theta”, 100,-10,180);
pt_cut->SetFillColor(kGreen);
TH1F *nt_cut = new TH1F(“nt_cut”,“Neutron Theta”, 100,-10,180);
nt_cut->SetFillColor(kGreen);
TH1F *kmpx_cut = new TH1F(“kmpx_cut”,“Km Px”, 100,-3,3);
kmpx_cut->SetFillColor(kGreen);
TH1F *kppx_cut = new TH1F(“kppx_cut”,“Kp Px”, 100,-3,3);
kppx_cut->SetFillColor(kGreen);
TH1F *ppx_cut = new TH1F(“ppx_cut”,“P Px”, 100,-3,3);
ppx_cut->SetFillColor(kGreen);
TH1F *npx_cut = new TH1F(“npx_cut”,“N Px”, 100,-3,3);
npx_cut->SetFillColor(kGreen);
TH1F *kmpy_cut = new TH1F(“kmpy_cut”,“Km Py”, 100,-2,2);
kmpy_cut->SetFillColor(kGreen);
TH1F *kppy_cut = new TH1F(“kppy_cut”,“Kp Py”, 100,-2,2);
kppy_cut->SetFillColor(kGreen);
TH1F *ppy_cut = new TH1F(“ppy_cut”,“P Py”, 100,-2,2);
ppy_cut->SetFillColor(kGreen);
TH1F *npy_cut = new TH1F(“npy_cut”,“N Py”, 100,-2,2);
npy_cut->SetFillColor(kGreen);
TH1F *kmpz_cut = new TH1F(“kmpz_cut”,“Km Pz”, 100,0,3);
kmpz_cut->SetFillColor(kGreen);
TH1F *kppz_cut = new TH1F(“kppz_cut”,“Kp Pz”, 100,0,3);
kppz_cut->SetFillColor(kGreen);
TH1F *ppz_cut = new TH1F(“ppz_cut”,“P Pz”, 100,0,3);
ppz_cut->SetFillColor(kGreen);
TH1F *npz_cut = new TH1F(“npz_cut”,“N Pz”, 100,0,3);
npz_cut->SetFillColor(kGreen);
int nCut1 = 0;
double neutronpx_phi;
double protonpx_phi;
double kpluspx_phi;
double kminuspx_phi;
double neutronpy_phi;
double protonpy_phi;
double kpluspy_phi;
double kminuspy_phi;
double neutronpz_phi;
double protonpz_phi;
double kpluspz_phi;
double kminuspz_phi;
double neutront_phi;
double protont_phi;
double kminust_phi;
double kplust_phi;
double kminus_phi;
double kpluss_phi;
double protons_phi;
double neutrons_phi;
double kminus_rho;
double kplus_rho;
double proton_rho;
double neutron_rho;
double kminuss_phi;
double kplus_phi; //??
double proton_phi;
double neutron_phi;
double kminus_theta;
double kplus_theta;
double proton_theta;
double neutron_theta;
double kminuss_px;
double kplus_px;
double proton_px;
double neutron_px;
double kminus_py;
double kplus_py;
double proton_py;
double neutron_py;
double kminus_pz;
double kplus_pz;
double proton_pz;
double neutron_pz;
for (Int_t n=0; n<NrOfEvents; n++) { //instead of being 100000 events this will be 100000 runs where each event is ‘tested’ to see if it is detected
if( n % 1000 == 0 )
cout << n << endl;
Double_t weight = event.Generate();
TLorentzVector *pNeutron = event.GetDecay(1);
TLorentzVector *pKp = event.GetDecay(2);
TLorentzVector *pKm = event.GetDecay(3);
TLorentzVector *pProton = event.GetDecay(0);
TLorentzVector pNKp = *pNeutron + *pKp;
TLorentzVector pKmP = *pKm + *pProton;
TLorentzVector pKms = *pKm;
TLorentzVector pNs = *pNeutron;
TLorentzVector pPs = *pProton;
TLorentzVector pKps = *pKp;
dist->Fill(pNKp.M2() ,pKmP.M2() ,weight); //histogram of momentum for the two pairs of particles
h1->Fill(pKms.M2());//1D histogram of individual particle mass^2
h3->Fill(pNs.M2());
h4->Fill(pPs.M2());
h5->Fill(pKps.M2());
kminus_phi = pKms.Phi()*TMath::RadToDeg();
// cout << kminus_phi << endl;
if (kminus_phi < 0.0 || kminus_phi > phiMin) {
nCut1++;
kmm2_cut->Fill(pKms.M2());
}
protons_phi = pPs.Phi()*TMath::RadToDeg();
// cout << protons_phi << endl;
if (protons_phi < 0.0 || protons_phi > phiMin) {
nCut1++;
pm2_cut->Fill(pPs.M2());
}
neutrons_phi = pNs.Phi()*TMath::RadToDeg();
// cout << neutrons_phi << endl;
if (neutrons_phi < 0.0 || neutrons_phi > phiMin) {
nCut1++;
nm2_cut->Fill(pNs.M2());
}
kpluss_phi = pKps.Phi()*TMath::RadToDeg();
// cout << kpluss_phi << endl;
if (kpluss_phi < 0.0 || kpluss_phi > phiMin) {
nCut1++;
kpm2_cut->Fill(pKps.M2());
}
kminus_rho = pKms.Rho();
h6->Fill(pKms.Rho());
kplus_rho = pKps.Rho();
h7->Fill(pKps.Rho());
proton_rho = pPs.Rho();
h8->Fill(pPs.Rho());
neutron_rho = pNs.Rho();
h9->Fill(pNs.Rho());
kminuss_phi = pKms.Phi();
h10->Fill(pKms.Phi()*TMath::RadToDeg());
kplus_phi = pKps.Phi();
h11->Fill(pKps.Phi()*TMath::RadToDeg());
proton_phi = pPs.Phi();
h12->Fill(pPs.Phi()*TMath::RadToDeg());
neutron_phi = pNs.Phi();
h13->Fill(pNs.Phi()*TMath::RadToDeg());
kminus_theta = pKms.Theta();
h14->Fill(pKms.Theta()*TMath::RadToDeg());
kplus_theta = pKps.Theta();
h15->Fill(pKps.Theta()*TMath::RadToDeg());
proton_theta = pPs.Theta();
h16->Fill(pPs.Theta()*TMath::RadToDeg());
neutron_theta = pNs.Theta();
h17->Fill(pNs.Theta()*TMath::RadToDeg());
h18->Fill(pKps.M2() ,pKms.M2() ,weight);
h19->Fill(pNs.M2() ,pPs.M2() ,weight);
kminuss_px = pKms.Px();
h20->Fill(pKms.Px());
kplus_px = pKps.Px();
h21->Fill(pKps.Px());
proton_px = pPs.Px();
h22->Fill(pPs.Px());
neutron_px = pNs.Px();
h23->Fill(pNs.Px());
kminus_py = pKms.Py();
h24->Fill(pKms.Py());
kplus_py = pKps.Py();
h25->Fill(pKps.Py());
proton_py = pPs.Py();
h26->Fill(pPs.Py());
neutron_py = pNs.Py();
h27->Fill(pNs.Py());
kminus_pz = pKms.Pz();
h28->Fill(pKms.Pz());
kplus_pz = pKps.Pz();
h29->Fill(pKps.Pz());
proton_pz = pPs.Pz();
h30->Fill(pPs.Pz());
neutron_pz = pNs.Pz();
h31->Fill(pNs.Pz());
kminust_phi = pKms.Phi()*TMath::RadToDeg();
// cout << kminust_phi << endl;
if (kminust_phi < 0.0 || kminust_phi > phiMin) {
nCut1++;
kmt_cut->Fill(pKms.Theta()*TMath::RadToDeg());
}
kplust_phi = pKps.Phi()*TMath::RadToDeg();
// cout << kplust_phi << endl;
if (kplust_phi < 0.0 || kplust_phi > phiMin) {
nCut1++;
kpt_cut->Fill(pKps.Theta()*TMath::RadToDeg());
}
protont_phi = pPs.Phi()*TMath::RadToDeg();
// cout << protont_phi << endl;
if (protont_phi < 0.0 || protont_phi > phiMin) {
nCut1++;
pt_cut->Fill(pPs.Theta()*TMath::RadToDeg());
}
neutront_phi = pNs.Phi()*TMath::RadToDeg();
// cout << neutront_phi << endl;
if (neutront_phi < 0.0 || neutront_phi > phiMin) {
nCut1++;
nt_cut->Fill(pNs.Theta()*TMath::RadToDeg());
}
kminuspx_phi = pKms.Phi();
// cout << kminuspx_phi << endl;
if (kminuspx_phi < 0.0 || kminuspx_phi > phiMin) {
nCut1++;
kmpx_cut->Fill(pKms.Px());
}
kpluspx_phi = pKps.Phi();
// cout << kpluspx_phi << endl;
if (kpluspx_phi < 0.0 || kpluspx_phi > phiMin) {
nCut1++;
kppx_cut->Fill(pKps.Px());
}
protonpx_phi = pPs.Phi();
// cout << protonpx_phi << endl;
if (protonpx_phi < 0.0 || protonpx_phi > phiMin) {
nCut1++;
ppx_cut->Fill(pPs.Px());
}
neutronpx_phi = pNs.Phi();
// cout << neutronpx_phi << endl;
if (neutronpx_phi < 0.0 || neutronpx_phi > phiMin) {
nCut1++;
npx_cut->Fill(pNs.Px());
}
kminuspy_phi = pKms.Phi();
// cout << kminuspy_phi << endl;
if (kminuspy_phi < 0.0 || kminuspy_phi > phiMin) {
nCut1++;
kmpy_cut->Fill(pKms.Py());
}
kpluspy_phi = pKps.Phi();
// cout << kpluspy_phi << endl;
if (kpluspy_phi < 0.0 || kpluspy_phi > phiMin) {
nCut1++;
kppy_cut->Fill(pKps.Py());
}
protonpy_phi = pPs.Phi();
// cout << protonpy_phi << endl;
if (protonpy_phi < 0.0 || protonpy_phi > phiMin) {
nCut1++;
ppy_cut->Fill(pPs.Py());
}
neutronpy_phi = pNs.Phi();
// cout << neutronpy_phi << endl;
if (neutronpy_phi < 0.0 || neutronpy_phi > phiMin) {
nCut1++;
npy_cut->Fill(pNs.Py());
}
kminuspz_phi = pKms.Phi();
// cout << kminuspz_phi << endl;
if (kminuspz_phi < 0.0 || kminuspz_phi > phiMin) {
nCut1++;
kmpz_cut->Fill(pKms.Pz());
}
kpluspz_phi = pKps.Phi();
// cout << kpluspz_phi << endl;
if (kpluspz_phi < 0.0 || kpluspz_phi > phiMin) {
nCut1++;
kppz_cut->Fill(pKps.Pz());
}
protonpz_phi = pPs.Phi();
// cout << protonpz_phi << endl;
if (protonpz_phi < 0.0 || protonpz_phi > phiMin) {
nCut1++;
ppz_cut->Fill(pPs.Pz());
}
neutronpz_phi = pNs.Phi();
// cout << neutronpz_phi << endl;
if (neutronpz_phi < 0.0 || neutronpz_phi > phiMin) {
nCut1++;
npz_cut->Fill(pNs.Pz());
}
}
//cout << nCut1 << endl;
c2->cd(1);
dist->Draw(“col”);
c1->cd(1);
h1->Draw();
kmm2_cut->Draw(“same”);
c1->cd(2);
h3->Draw();
pm2_cut->Draw(“same”);
c1->cd(3);
h4->Draw();
nm2_cut->Draw(“same”);
c1->cd(4);
h5->Draw();
kpm2_cut->Draw(“same”);
c1->cd(5);
h6->Draw();
c1->cd(6);
h7->Draw();
c1->cd(7);
h8->Draw();
c1->cd(8);
h9->Draw();
c1->cd(9);
h10->Draw();
c1->cd(10);
h11->Draw();
c1->cd(11);
h12->Draw();
c1->cd(12);
h13->Draw();
c1->cd(13);
h14->Draw();
kmt_cut->Draw(“same”);
c1->cd(14);
h15->Draw();
kpt_cut->Draw(“same”);
c1->cd(15);
h16->Draw();
pt_cut->Draw(“same”);
c1->cd(16);
h17->Draw();
nt_cut->Draw(“same”);
c2->cd(2);
h18->Draw(“col”);
c2->cd(3);
h19->Draw(“col”);
c3->cd(1);
h20->Draw();
kmpx_cut->Draw(“same”);
c3->cd(2);
h21->Draw();
kppx_cut->Draw(“same”);
c3->cd(3);
h22->Draw();
ppx_cut->Draw(“same”);
c3->cd(4);
h23->Draw();
npx_cut->Draw(“same”);
c3->cd(5);
h24->Draw();
kmpy_cut->Draw(“same”);
c3->cd(6);
h25->Draw();
kppy_cut->Draw(“same”);
c3->cd(7);
h26->Draw();
ppy_cut->Draw(“same”);
c3->cd(8);
h27->Draw();
npy_cut->Draw(“same”);
c3->cd(9);
h28->Draw();
kmpz_cut->Draw(“same”);
c3->cd(10);
h29->Draw();
kppz_cut->Draw(“same”);
c3->cd(11);
h30->Draw();
ppz_cut->Draw(“same”);
c3->cd(12);
h31->Draw();
npz_cut->Draw(“same”);
}