Edward Williams

January 28, 2016:

Finally successfully created logbook!

Working on Pythia.

March 9, 2016: Unfortunately I have missed some logbook entries; I will be careful not to do this in future. I am now working on a jet clustering exercise with my research group for our project with Dr. Eno. I am making progress on understanding the code; I will add more after our meeting with Dr. Eno today.

Homework 1:

(I produced some graphs, but they were not visible when I put my code onto this page; they are now visible in the file at the bottom)

I have been working on the assignment 1 pythia tutorial. I used git.clone to copy Professor Eno’s code from github, as the tutorial instructed. I had some difficulties using git clone at first, which I solved with the help of classmates. I used

source PYTHIA_setup.csh

in the directory, which allows it to run root. I used make

pythiaTree

and

./pythiaTree >& output.txt

to test the code and see if I had downloaded it correctly; it produced histograms and a tree as intended. I then changed the code to store the particle’s mass in the tree, by adding the following lines (in bold):

…

Float_t Apz[5000];

Float_t Am[5000];

Float_t ApT[5000];

…

… 

 t1->Branch("Apy",Apx,"Apy[nTot]/F");

 t1->Branch("Apz",Apx,"Apz[nTot]/F");

 t1->Branch("Am",Apx,"Am[nTot]/F");

…

…

Apy[nTot-1]=pythia.event[i].py();

Apz[nTot-1]=pythia.event[i].pz();

Am[nTot-1]=pythia.event[i].m();

…

I also changed the number of events run by changing the 100 to 1000 in the for loop: for (int iEvent = 0; iEvent < 100; ++iEvent) {

I changed it back to a smaller number later to speed up the compiling and running during debugging.

After I made these changes, no change was reflected in the histograms or leaves. I struggled with this until I discovered that the lines

make pythiaTree

and

./ pythiaTree >& output.txt

were responsible for compiling and running the file, and had to be used after every change in the code for it to be represented in the graphs. I had thought that they only had to be used once, at the beginning of the assignment. Once I used these after every change, my coding was reflected in the graphs, and the mass and event number change were shown successfully.

The second part of the code gave me difficulties, which I struggled with at length. However, I eventually managed to resolve the difficulties.

I copied and pasted the lines of code from the tutorial in place of the QCD section of the code, as instructed. I looked at the pythia8 documentation at http://home.thep.lu.se/~torbjorn/pythia82html/Welcome.html, but did not find it very helpful; most of the lines of code from the tutorial were not there, although I could find some on other places on the website.

Next, I attempted to plot the invariant mass of the two highest pT electons in the event, hoping to see the bump at 1500 GeV that was expected. To do this, I created new variables at the beginning of the code:

int idOfSecondMax = 0;

int idOfFirstMax = 0;

Float_t Am[5000];

Float_t ApT[5000];

Float_t firstMax = -1000000;

Float_t secondMax = -1000000;

Float_t ipT;

Float_t tmtemp;

Float_t TargetMass[5000]; //invariant masses of the highest-pT electrons in each event

I also created a new histogram to contain the desired invariant masses found for each event:

TH1F *targetmass = new TH1F("targetmass","invariant mass of highest-pT electrons per array",1600, -500, 500);

I found that the histogram was at first too small to show up to 1500 GeV; I had not known that the last two number in the parentheses are the minimum and maximum values for the x-axis of the histogram. I adjusted these accordingly.

I designed a method of finding the two electrons with the highest pT in the event and determining their invariant mass. I put this code after the code to fill the x, y, and z branches, inside the loop that runs once for each particle in each event.

This is my code for this:

            //finding top pT electrons

            if (pythia.event[i].id() ==11 || pythia.event[i].id() == -11) {

              ipT = pythia.event[i].pT();

              if (ipT>secondMax && ipT<firstMax){

                secondMax = ipT;

                idOfSecondMax = i;

              }

              if (ipT>firstMax){

                secondMax = firstMax;

                idOfSecondMax = idOfFirstMax;

                firstMax = ipT;

                idOfFirstMax = i;

              }

            }

        //done finding

This code first uses event[i].id() == 11 or -11 to determine whether or not the particle being analyzed is an electron or positron. I added this code after looking at some other students’ logbooks for this; I had not known how to use event[i].id(), and had previously just been running this code on any charged particle instead of on any electron; this skewed my data significantly until I solved it.

My code then sets the variable ipT (pT of particle i), defined at the top of the code, to the transverse momentum of the particle. Using this variable is neater than retyping “pythia.event[i].pT();” every time, and reduces risk of typos. This variable is reset for every particle.

My code then checks runs two if statements, to check if ipT is

1.)  greater than secondMax but less than firstMax,

2.)  greater than firstMax, or

3.)  less than both

By the end of the loop, firstMax should be the maximum pT value for all particles in the loop, and secondMax should be the second greatest pT value. They are both initialized at -1,000,000, to ensure that they are not greater than the pTs of any particle. The first particle that is checked will have a pT greater than -1,000,000, so its firstMax will become its pT and indexOfFirstMax will become i, its index. The second particle will either become secondMax, or replace the first particle as firstMax if it has a greater pT. If a particle has a pT greater than firstMax, it replaces firstMax and firstMax replaces secondMax. In this way, it is sure that the greatest and second-greatest pT electrons for each event will be found.

After this, I originally had code to set the variable tmtemp (meaning “temporary target mass”) to be the invariant mass of the two particles, event[indexOfFirstMax] and event[indexOfSecondMax]. Then the code

 Targetmass -> Fill( tmtemp );

fills the histogram with the data.

I reviewed this code, thought it would work, and tried to compile it; it produced multiple errors. I hunted through my code and discovered that I had initialized a few variables incorrectly, and made a few typos (including an asymmetrical parenthesis and a , instead of a ;); I fixed these issues. After this, it compiled and ran, and produced a histogram, but it did not have a bump at 1500 GeV that I saw. I struggled with this, googling the issue and trying different approaches, but none worked. I eventually looked online at other students’ logbooks and how they had solved this problem; I discovered that:

1.)  Three of the numbers in the lines of code pasted from the tutorial were off by a factor of 10, due to unit issues.

2.)  I should graph the sum of the transverse momenta of the two highest-pT electrons of each event, instead of their invariant mass. This is contrary to what the tutorial said, but apparently produced the right results.

3.)  How to search specifically for electrons, instead of all charged particles, as described above.

When I made these changes, the code successfully produced a graph with a bump at 1500 GeV:

(1)

For part 2 of the homework, involving the Drell-Yan background, I used the lines of code from the tutorial, and consulted other students’ logbooks to see their approach to the problem. So as not to risk losing my code for part 1, I saved this code as pythiaTreeOriginal.cc and made the edits for part 2 in pythiaTree.cc. The code eventually ran, and produced this graph for the background:

(2)

Project Research:

I have found and read these sources, which relate to my group’s research project:

http://arxiv.org/pdf/1306.4676v2.pdf

http://arxiv.org/abs/1508.05327

http://www.pnas.org/content/112/40/12278.full

http://www.phys.virginia.edu/Files/fetch.asp?EXT=Seminars:2693:SlideShow

http://journals.aps.org/prd/abstract/10.1103/PhysRevD.93.043529

http://research.dsu.edu/cetup/documents/2013-talks/dark-matter/DM-25-YangBai.pdf

 (some graphs below do not show up properly; they are in files at the bottom of the page as well.)

Week 3

We discovered a glitch that meant that particles were not removed from the code after their beam and particle distances were calculated, meaning that the code ran indefinitely, analyzing the same particle, until it was fixed.

There was another error in the code

              px = ( ( (Particles[minindex_i].pt * cos((Particles[minindex_i].phi) ) ) + (Particles[minindex_j].pt * cos(Particles[minindex_j].phi) ) ) );

              py = ( ( (Particles[minindex_i].pt * sin((Particles[minindex_i].phi) ) ) + (Particles[minindex_j].pt * sin( Particles[minindex_j].phi) ) ) );

              pz = ( ( (Particles[minindex_i].pt * sinh(Particles[minindex_i].eta) ) ) + (Particles[minindex_j].pt * sinh(Particles[minindex_j].eta) ) );

where some i’s and j’s were replaced with each other, skewing the results.

Week 4/5

We have been working on finding and eliminating errors in our codes, and expanding what it can do. My code was calculating the energies of particles incorrectly.  I replaced the code with

      for (int wh = 0; wh < numentry - 1; wh++)

        {

           px = (*pt)[wh]*cos((*phi)[wh]);

           py = (*pt)[wh]*sin((*phi)[wh]);

           pz = (*pt)[wh]*sinh((*eta)[wh]);

          item.pt = (*pt)[wh];

          item.phi = (*phi)[wh];

          item.eta = (*eta)[wh];

          item.mass = (*mass)[wh];

          item.energy = pow((pow(px,2)+pow(py,2)+pow(pz,2)+pow(item.mass,2)),0.5);

          Particles.push_back( item );

        }

This correctly extracts the mass and momentum from the particle data struct, and calculates the energies correctly, using the E2 = (px2 +py2 + pz2 + m2)0.5 equation.

Week 6

The code produces histograms successfully, but the histograms for phi and eta look asymmetrical and not how they should, meaning that there must be an unfound bug in that part of the code. The error was caused by an eta written in place of a phi in the code,  which was found and fixed. The method of calculating these angles was also changed; the old method allowed multiple values for the angles where only one should have been possible, and may have been contributing to the problem. We are now supposed to make a histogram showing the number of particles per jet for the jets of relatively high pT’s.

Week 7/8

This code creates a new histogram to display the number of particles in each event:

TH1F *JetNumParticles = new TH1F("JetNumParticles", "Number of Particles in Jets with  pT greater than 50 GeV", 100, 0, 100);

  JetNumParticles -> SetMarkerStyle(4);

A particlecount variable was also added to the initial struct.

This is set to 1 for each particle in the event, and at the end of each jet loop the number of particles in it is calculated and saved to a histogram in the line

              JetNumParticles -> Fill( Jets.back().particlecount );

The line

  JetNumParticles -> Draw("");

was also added to the end of the code.

The code now successfully produces the histogram, showing the number of particles per jet for high pT’s.

(1)

This code only contains 100 events; it was set that way to accelerate debugging. The histogram from the code using the complete number of events is:

(2)

Week 8

For next meeting, Profesor Eno asked us to look at the Dark QCD code that she put on github and determine what each part of it does; also to create histograms to compare dark jets and normal jets. This code will use new events generated using pythia.

Week 9:

            We have finished our code for analyzing and comparing dark and normal jets. The code creates random events (using Pythia) producing both dark and normal quarks and shows their properties in histograms. We found that in some respects, such as the transverse momentum, the dark and normal matter jets are very similar (see image at the bottom to see a histogram comparing them; the spike at the very beginning is caused by a slight glitch in the code that our code relies on.) Professor Eno also explained closest approach (shortest distance between a charged particle moving in a circular path in a magnetic field and the point of initial collision) to us, and gave us a calculation for us to integrate into our code for next time. The closest approach should be much farther for dark particles than for regular matter particles.

Week 10:

            We have successfully added closest approach calculation into our code and produced histograms based on it. As we expected, the closest approach for normal matter daughter particles is approximately equivalent to zero (they are created at the initial collision point), while the closest approach for dark daughter particles varies much more and is, on average, much farther (the dark mediator particle is created at the collision point, then moves through the detector until it decays into dark daughter particles). See histograms at the bottom of the page. The code now produces all the histograms that we need for the poster.