#include "TPythia8.h" #include "TClonesArray.h" #include "TParticle.h" #include "TCanvas.h" #include "TH2F.h" #include #include #include using namespace Pythia8; using namespace std; //Breit-Wigner function (taken from https://twiki.cern.ch/twiki/pub/CMSPublic/WorkBookHowToFit/BW.C.txt) Double_t mybw(Double_t* x, Double_t* par) { Double_t arg1 = 14.0/22.0; // 2 over pi Double_t arg2 = par[1]*par[1]*par[2]*par[2]; //Gamma=par[1] M=par[2] Double_t arg3 = ((x[0]*x[0]) - (par[2]*par[2]))*((x[0]*x[0]) - (par[2]*par[2])); Double_t arg4 = x[0]*x[0]*x[0]*x[0]*((par[1]*par[1])/(par[2]*par[2])); return par[0]*arg1*arg2/(arg3 + arg4); } void macroName() { gSystem->Load("libEG"); gSystem->Load("libEGPythia8"); // **************************** // *** VARIABLE DECLARATION *** // **************************** TClonesArray* particles = new TClonesArray("TParticle", 1000); TPythia8 * pythia = new TPythia8(); //the mass bounds we use to determine if two particles are from a Z boson or not Double_t mass_bound_high = 99; Double_t mass_bound_low = 80; //TODO: Double-check units! //a histogram of all the masses of the histogram auto mass_hist = new TH1F("Z Boson Mass", "Z Boson Mass;Mass (GeV/c^2);Counts", 50, mass_bound_low, mass_bound_high); //a histogram of the magnitudes of the Z Boson momentum auto p_hist = new TH1F("Z Boson P", "Z Boson Momentum;Momentum (GeV);Counts", 50, 0, 3100); //the 3-momentum of the Z boson TVectorD zboson_p(3); //the transverse momentum of the two decay products Double_t Pt1 = 0; Double_t Pt2 = 0; //the pseudorapidity of the two decay products Double_t n1 = 0; Double_t n2 = 0; //the azimuthal angle of the two decay products Double_t phi1 = 0; Double_t phi2 = 0; //the invariant mass of the Z boson Double_t zboson_m = 0; //magnitude of the momentum of the Z boson Double_t zboson_p_mag = 0; //the maximum momentum (for debugging) Double_t max_p = 0; //the maximum mass (for debugging) Double_t max_m = 0; //resonance width of a Z boson (for the Breit-Wigner fit) Double_t resonance_width = 0.083984; //generally accepted mass for a Z boson (for the Breit-Wigner fit) Double_t Z_mass = 91; //number of events size_t N = 1000; //the number of Z bosons we found size_t zboson_count = 0; // ************* // *** SETUP *** // ************* //Go all the way down to zero transverse momentum pythia->ReadString("PhaseSpace:pTHatMin = 0."); //Drell-Yan process of quark-antiquark annihilation to a virtual photon or a Z boson pythia->ReadString("WeakSingleBoson::ffbar2gmZ = on"); //Compare the Z boson azimuthal angular distribution at low Z boson pT for primordial //transverse momentum on vs off pythia->ReadString("BeamRemnants:primordialKT=on"); //p+p fixed-target collision pythia->Initialize(2212, 2212, 6500, 6500); //print graphing results gStyle->SetOptFit(1111); // ******************** // *** COLLECT DATA *** // ******************** //Event loop for(size_t i = 0; i < N; ++i) { //generate events pythia->GenerateEvent(); //import the particles into our array of particles pythia->ImportParticles(particles,"All"); Int_t np = particles->GetEntriesFast(); //particle loop for (Int_t ip = 0; ip < np; ++ip) { //store the particle we're working with TParticle *part = (TParticle*) particles->At(ip); Int_t ist = part->GetStatusCode(); //Positive codes are final particles if (ist <= 0) continue; Int_t id_i = part->GetPdgCode(); //if this particle is a mu... if (id_i == 13) { //reset mu and anti mu data Pt1 = 0; Pt2 = 0; phi1 = 0; phi2 = 0; n1 = 0; n2 = 0; //reset three-momentum of zboson zboson_p[0] = 0; zboson_p[1] = 0; zboson_p[2] = 0; //get data on the mu Pt1 = part->Pt(); phi1 = part->Phi(); n1 = part->Eta(); //inner loop used to find matching anti mu for (Int_t ip2 = 0; ip2 < np; ++ip2) { TParticle *part2 = (TParticle*) particles->At(ip2); Int_t ist2 = part2->GetStatusCode(); //Positive codes are final particles if (ist2 <= 0) continue; Int_t id_i2 = part2->GetPdgCode(); //if this particle is an anti-mu... if (id_i2 == -13) { //collect data on the anti-mu Pt2 = part2->Pt(); phi2 = part2->Phi(); n2 = part2->Eta(); //calculate the mass of the potential Z boson zboson_m = sqrt(2*Pt1*Pt2*(cosh(n1-n2) - cos(phi1-phi2))); //check if this is a valid Z boson if (zboson_m >= mass_bound_low && zboson_m <= mass_bound_high) { //calculate the Z boson 3-momentum zboson_p[0] += part->Px() + part2->Px(); zboson_p[1] += part->Py() + part2->Py(); zboson_p[2] += part->Pz() + part2->Pz(); //calculate the magnitude of the Z boson momentum zboson_p_mag = sqrt(zboson_p[0] * zboson_p[0] + zboson_p[1] * zboson_p[1] + zboson_p[2] * zboson_p[2]); //fill histograms mass_hist->Fill(zboson_m); p_hist->Fill(zboson_p_mag); //Find some info for debugging: //increase the number of Z bosons we found ++zboson_count; //find the maximum Z boson mass and momentum if (zboson_m > max_m) max_m = zboson_m; if (zboson_p_mag > max_p) max_p = zboson_p_mag; //the appropriate pair has been found, //break the inner loop break; }//if }//if }//for (second particle loop) }//if }//for (first particle loop) }//for (event loop) // ****************** // *** GRAPH DATA *** // ****************** //create canvas auto c = new TCanvas("Z Boson Data", "Z Boson Data", 1200, 600); c->Divide(1, 2); //graph mass histogram and fit Breit-Wigner to it c->cd(1); TF1 *func = new TF1("mybw",mybw,mass_bound_low, mass_bound_high, 3); func->SetParameter(1, resonance_width); //get gamma/resonance width for parameter 1 func->SetParameter(2, mass_hist->GetMean()); //get M/Z mass for parameter 2 auto mass_fit = mass_hist->Fit("mybw", "QR"); //mass_fit->DrawClone("Same"); mass_hist->DrawClone("Same"); //graph momentum histogram c->cd(2); p_hist->DrawClone("Same"); //print the number of found Z bosons and maximum/minimum momentum for debugging purposes cout << "We found " << zboson_count << " Z bosons!\n"; cout << "Maximum mass: " << max_m << " & Maximum momentum: " << max_p << "\n"; }