Hi All,
This code doesn’t get any errors or warnings, but won’t finish running unless I remove the last few graphs. It creates the canvasses but doesn’t print anything onto them. I can’t see anything obvious that strikes me as wring- any suggestions would be much appreciated!
This is the code:
void test7_1(UInt_t NrOfEvents = 1000000) {
if (!gROOT->GetClass("TGenPhaseSpace")) gSystem->Load("libPhysics");
TLorentzVector target(0.0, 0.0, 0.0, 1.879);
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.493, 0.493, 0.938} ;
/* 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);
TH1F *h1 = new TH1F("h1","kminus mass^2", 100,0,3);
h1->SetFillColor(kRed);
TH1F *h3 = new TH1F("h3","Proton mass^2", 100,0,3);
h3->SetFillColor(kRed);
TH1F *h4 = new TH1F("h4","Neutron mass^2", 100,0,3);
h4->SetFillColor(kRed);
TH1F *h5 = new TH1F("h5","Kplus mass^2", 100,0,3);
h5->SetFillColor(kRed);
TH1F *h6 = new TH1F("h6","Rho Km",100,0,3);
TH1F *h7 = new TH1F("h7","Rho Kp",100,0,3);
TH1F *h8 = new TH1F("h8","Rho Proton",100,0,3);
TH1F *h9 = new TH1F("h9","Rho Neutron",100,0,3);
TH1F *h10 = new TH1F("h10","Phi Km",100,-180,180);
TH1F *h11 = new TH1F("h11","Phi Kp",100,-180,180);
TH1F *h12 = new TH1F("h12", "Phi P",100,-180,180);
TH1F *h13 = new TH1F("h13", "Phi N",100,-180,180);
TH1F *h14 = new TH1F("h14", "Km Theta",100,-10,180);
h14->SetFillColor(kRed);
TH1F *h15 = new TH1F("h15", "Kp Theta",100,-10,180);
h15->SetFillColor(kRed);
TH1F *h16 = new TH1F("h16", "Proton Theta",100,-10,180);
h16->SetFillColor(kRed);
TH1F *h17 = new TH1F("h17", "Neutron Theta",100,-10,180);
h17->SetFillColor(kRed);
TH2F *h18 = new TH2F("h18", "Kp + Km m^2",100,0.2,0.3, 100,0.2,0.3);
TH2F *h19 = new TH2F("h19", "P + N m^2",100,0.8,1, 100,0.8,1);
TH1F *h20 = new TH1F("h20", "Km Px",100,0,3);
h20->SetFillColor(kRed);
TH1F *h21 = new TH1F("h21", "Kp Px",100,0,3);
h21->SetFillColor(kRed);
TH1F *h22 = new TH1F("h22", "P Px",100,0,3);
h22->SetFillColor(kRed);
TH1F *h23 = new TH1F("h23", "N Px",100,0,3);
h23->SetFillColor(kRed);
TH1F *h24 = new TH1F("h24", "Km Py",100,0,3);
h24->SetFillColor(kRed);
TH1F *h25 = new TH1F("h25", "Kp Py",100,0,3);
h25->SetFillColor(kRed);
TH1F *h26 = new TH1F("h26", "P Py",100,0,3);
h26->SetFillColor(kRed);
TH1F *h27 = new TH1F("h27", "N Py",100,0,3);
h27->SetFillColor(kRed);
TH1F *h28 = new TH1F("h28", "Km Pz",100,0,3);
h28->SetFillColor(kRed);
TH1F *h29 = new TH1F("h29", "Kp Pz",100,0,3);
h29->SetFillColor(kRed);
TH1F *h30 = new TH1F("h30", "P Pz",100,0,3);
h30->SetFillColor(kRed);
TH1F *h31 = new TH1F("h31", "N Pz",100,0,3);
h31->SetFillColor(kRed);
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,0,3);
kmpx_cut->SetFillColor(kGreen);
TH1F *kppx_cut = new TH1F("kppx_cut","Kp Px", 100,0,3);
kppx_cut->SetFillColor(kGreen);
TH1F *ppx_cut = new TH1F("ppx_cut","P Px", 100,0,3);
ppx_cut->SetFillColor(kGreen);
TH1F *npx_cut = new TH1F("npx_cut","N Px", 100,0,3);
npx_cut->SetFillColor(kGreen);
TH1F *kmpy_cut = new TH1F("kmpy_cut","Km Py", 100,0,3);
kmpy_cut->SetFillColor(kGreen);
TH1F *kppy_cut = new TH1F("kppy_cut","Kp Py", 100,0,3);
kppy_cut->SetFillColor(kGreen);
TH1F *ppy_cut = new TH1F("ppy_cut","P Py", 100,0,3);
ppy_cut->SetFillColor(kGreen);
TH1F *npy_cut = new TH1F("npy_cut","N Py", 100,0,3);
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 kpluss_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
Double_t weight = event.Generate();
TLorentzVector *pNeutron = event.GetDecay(0);
TLorentzVector *pKp = event.GetDecay(1);
TLorentzVector *pKm = event.GetDecay(2);
TLorentzVector *pProton = event.GetDecay(3);
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());
kpluss_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");
}
Thanks in advance!