Dark Matter Project

2/17/16

• Today was my first meeting with Dr. Eno and she assigned us the task of making a jet clustering algorithm.

• Before I go on, here are some links

  • https://sites.google.com/site/honr268neno/home

  • file:///Users/jamesnatoli/Downloads/ElementaryExperimentalHiggsPhysics%20(12).pdf

  • https://github.com/saraheno/darkgen

  • http://pdg.lbl.gov/2002/montecarlorpp.pdf

• I'll first go through the basics of Dark Matter and what I learned from Dr. Eno's ppt, the scientific paper, and my talk today with Dr. Eno

• Dark Matter is a proposed alternative form of matter that is non-baryonic, meaning that it is not made out of quarks. It is not just a random theory but it is well known that another form of matter must exist. Here is what is known for 99% sure

  • There is some missing mass in the universe. A lot actually; ~25% of the mass is invisible to us

  • This "other" matter does not interact via E&M or the Strong force or else we would see it

  • It may interact via the weak force and it definitely interacts via gravity

  • It is not quite so obvious as regular matter being attracted to "empty" areas of space...

• Big evidence is in the rotation curves of galaxies; Galaxies, based on the matter that we can see, are less dense on the outer edges and therefore should rotate slower as a function of the radius

• That is not the case however, as experimental data has shown us that the rotational velocity increases rapidly and then flatlines but does not decrease. The best explanation for this is that there is a halo of dark matter at the edges of galaxies that we cannot see

• The video below shows the differences. The left galaxy is what is predicted and the right galaxy is what is observed

  • https://www.youtube.com/watch?v=-qlV5hjIJm0

  • Another piece of evidence is found in gravitational lensing. Galaxy clusters are extremely massive and therefore bend light slightly around them. By studying the distortion of the light, the mass of the cluster can be calculated

  • The problem was that these calculations were much larger than how much matter was actually in that area. 

  • So scientists know that there is some other type of matter that interacts with Standard Model (SM) matter via gravity but not E&M or the strong force

  • Our project is particularly interested in how Dark Matter might interact with itself via a supposed Dark strong force, called Dark Quantum Chromodynamics or Dark QCD

  • Our best guess is that Dark Matter mimics SM matter and is composed of "dark quarks" that interact with a "dark gluon"

  • It is also proposed that a very heavy dark mediator particle could be formed in energetic collisions in the TeV range and decay in to a SM quark and a dark quark.

• • This process is related to and based on the jets that occur in detectors from SM quarks

  • Basically, the quarks created in the collision are flying apart faster than the string force can contain them. The strong force increases as the distance increases so the separation causes more quarks and anti-quarks to form from the energy

  • This process is called Fragmentation

  • It is followed by a process called Hadronization in which the quarks and anti-quarks collect in to color neural mesons, the most basic being a pion which form in to jets

  • This process is what the supposed dark quarks do, but with an added step, outlined below

• • The Dark Mediator must be heavy (because we haven't discovered it yet) and therefore probably has a longer decay time, but it is unknown

  • The jets would "emerge" at some unknown point in the detector, but farther out than a regular jet because of the added step of the dark mediator

  • Our task for this week is to write jet clustering algorithm that will sort through the dark matter event and organize the ~500 particles at the end in to 5-10 jets based on their momentum.

Jet Clustering (2/17/15)

• I imported Dr. Eno's code with the command

  • git clone https://github.com/saraheno/darkgen.git

• It created a new directory called darkgen with all of the necessary files in it

• Make sure to run "source PYTHIA_setup.csh" every time

• When I tried to run ./pythiaTree.cc >& log.txt it said permission denied so thats not good

More Jet Clustering (2/29/15)

• As is turns out, a lot of my issues were very silly ones

• I was going about it the whole wrong way... here is the correct way to do it

• The jet clustering file is a root file so to access the data, you should go through root. The easiest way is to follow the directions from HW9, shown below

Enter root. Then do the following commands.

            TFile *f = new TFile(“output.root”);

            f->ls();

            t1->MakeClass(“MyClass”);

now exit root. You should see two new files in your area, MyClass.h and MyClass.C.

• Instead of "output.root" put the name of the root file and instead of "MyClass" put whatever you want the new files to be called

• This puts everything nice and neat for you. Most of my work will be done in JetCLuster.C

• I am currently adding a loop to the end of that file. It will loop through every particle pair and compare it's distance with the paired particle and the beam by using this equation

• Dij is the distance from on particle to the other and DiB is the distance from the beam to the particle. R is some given value, in our case, 0.4

• The basic set tup is a nested for loop with the algorithm in the loop. We are supposed to beware of the eta-wrap issue where the eta value is 350 instead of -10

• Here are some plots we are supposed to be doing

  •     # of jets vs # of events with that number of jets

  •      pT spectrum

  •      eta spetrum

  •      phi spectrum

  •      # of particles vs # of jets with that number of particles

  •      calculate size of cone with 90% of the energy

  •    These are the names and types of the variables

    •    vector<float>   *eta;

    •    vector<float>   *phi;

    •    vector<float>   *mass;

  •    vector<float>   *pt;

  • They can be accessed with the command "eta[i]"

  • I've written most of the loop but I am struggling with the part where you have to add the particles 4 momenta and then add the new particle to the list. I'm not sure which list and how to access the 4 momenta

Even More Jet Clustering (3/8/16)

• Working with Edward and Mayowa today, we were able to solve some problems

• First, We figured out how to get the 4 momenta. The value we have is pT or the momentum in the transverse plane so we need to convert that in to the total momentum, and then convert that in to x, y, and z components by multiplying by the appropriate spherical values ie. cos(theta), sin(phi)

• The 4 momenta is (i think) [px, py, pz, mass]

• the equations for eta and phi are required for this, and that pT = psin(theta). Use this equation to find the total p, and then use that to find the components

• Ok... I have finally finished to the point where I can try to compile and see what kind of errors I'm looking at

• Update: I have yet to solve the phi wrap problem

• I had to be pretty creative to find the new values of phi, eta, and pT. To find eta i used the equation eta = -ln(tan(theta/2)) and pT = Psin(theta). Also, pT = sqrt(px^2 + py^2) and P = sqrt(px^2 + py^2 + pz^2)

• To find phi, i used px and rearranged the equation to read phi = acos(px/psin(theta))

• I also added a while loop to keep the for loops that evaluate every pair of particles looping as long as there are still entries. I added a variable pleasestop to which the for loops evaluate to correct for the loss of a particle every time the loops evaluate, so basically it decrements ie. pleasestop--;

• After compiling, the error I got was

  • • Symbol fChain is not defined in current scope  JetCluster.C:63:

    • Error: Failed to evaluate fChain->GetEntriesFast()

• So I guess I try to fix that now and work my way through all of the errors that I'm bound to get

• Before i forget, I also used the structs in a way that might be incorrect. Instead of making an array of structs, I made four float arrays in the struct and two objects of type P (the struct name). I chose this way instead of making the variables in the struct just floating point variables and making an array of structs. I'm not sure which way will work, hopefully my way because that is a lot of things to change. 

• The one thing I was worried about with my way was that because I have two struct objects, particles and jets, will they be allowed to have differing values? for instance, can jets.pt[12] be different than particles.pt[12]? hopefully it can and I read some documentation to help

• To help me learn about structs, I literally just googled structs c++ and read all of the links, the cplusplus.com website was most helpful, as it usually is

Debugging the Code 3/9/16

• So yes I'm skipping around a bit, but this section will primarily focus on the long and tedious debugging process. I'll go through every error I get and explain my solution

  1. "Symbol fChain is not defined in current scope  JetCluster.C:63"

     1. Solved by moving the bracket enclosing the Loop() method to the end of my code, it was previously excluding everything i had typed, oops

        The good news is .L JetCluster.C worked, and so did JetCluster m but I got plenty of errors when running m.Loop()

1. "Warning in <TFile::Append>: Replacing existing TH1: histo1 (Potential memory leak). 

2. Error: Illegal array dimension (Ignore subsequent errors) JetCluster.C:65:'

   1. The warning I'll ignore for now because it can still evaluate but the array dimension thing was something I was very worried about... maybe the array size can't be defined as a Long64_t type so i'll try typecasting first.

   2. Ok so something weird happened where everything crashed and it was scary so I'm going to try again...

   3. Ok no crash this time but I got the same error. I'll try typecasting to (int) 

   4. No such luck... same error and the curious thing is that the error is in line 65 which is where the struct is declared...

   5. The weird thing is that i added a line earlier but the error is still in the same spot, now in a commented line for some reason. 

Changes to Make 3/9/16

• so after our meeting with Dr. Eno, I realized I was doing a lot of things wrong. Actually most things. RIP me

• anyways, here are the changes to make

  • Use struct better, as in have the variables in the struct just floats and the struct object be the list

  • Use lists instead of arrays b/c they have a variable size and makes adding/removing so so so much easier

  • A list will also solve the problem of the loop having to go to a smaller number every time

  • Make functions for some of the complicated mathematical computations

  • Make the phi wrap work (it's not as hard as I though it was actually, just a simple if statement. The values are in radians, which is obvious and I should've known, silly me

• I'll go through making these changes, and then get back to debugging because I'm sure I will have a bucket of issues to deal with once I try to compile again.

Making Changes 3/23/16

• Dr. Eno could not meet this week but I've been making progress :) 

• Here are some things that I've learned...

• So the list is great but it's main drawback is that it has no [] operator like an array, so in order to access a specific element, one needs to "iterate" through the list starting at the beginning, or end and --ing or ++ing

• So basically, everytime I need to change the values of phi, pt, eta, or mass, I need to use a for loop and iterate through... I think

• The list does make things easier, essentially I set the struct objects values, and then add the struct to the list. The list contains struct objects with have 4 variables

• I've also made functions for the complicated math problems to make the code a little cleaner and easier

• For the phi wrap problem, I think I am going to use a function to change all the values before I add them to the list

• Note on iterators. 

  • In order to access the member variables stored in the struct, you have to put parenthesis around the iterator pointer, like so (*iterator).member_variable for it to work

  • Example...

• Just another note on this, "it->height" would also do the same thing

• This prints the height values of every struct object stored in the "Person" list

• Something to remind myself of later... I'm wondering if there is a better way to find the smallest Dij and Dib. I'd like to somehow make use of lists because they have a nice sort function so I'm thinking I might make two more lists of structs, with the member variables being the *it or *jt and the float Dib or Dij value but im not sure how that would work. Something to go back later for

UPDATE 2:23am

• I am now at a point where I can compile again, I have a function to fix the phi wrap problem, the struct is all worked out, I think everything is good

• I forgot how to compile so I'll throw it in here so I don't forget again

  • enter root

  • .L <filename.C>

  • <filename> m

  • m.Loop()

• My first compile attempt said too many brackets so ill have to scrounge through and see what happened

• Aaaaaaaaand I tried again and got about a hundred errors so lets see if I can clean them up by fixing the bracket issue

• A lot of them are things not being declared in the scope which is easy to fix... hopefully

Continuing to Debug-3/28/16

• After making the above changes i am now slogging through the errors

• A lot of them dealt with pointers bc as it turns out, the vectors from the ".h" file are actually pointers, they are named *pt, *eta...

• So in order to access them you need to "dereference" the value by calling (*pt)[3] in order to access an element or you can use the arrow (->) in order to access things

• This is similar to the iterator issue

• I also moved all of my functions and variable declarations outside of the loop at the beginning of the file to avoid scope issues

• As Dr. Eno mentioned, the loop

• • is what actually does stuff. Every jentry is an event so everything needs to be in side of the outer for loop

  • Iterator note: iterator can be incremented by ++ing but not by +1 or +8, but instead by the next method which is implemented like so

    • iterarator1 = next( iterator2, 1 )

    • in this case, iterator1 and 2 are iterators and the second argument that next takes is the difference between the iterators so the line above assigns iterator1 to iterator2 + 1 but if the number was a 5, the iterator1 would be assigned to iterator2 + 5

  • Ok so forget everything I said about the next() method because apparently the version of C++ that we have doesn't have the method so I will have to make a new iterator, assign it to the first iterator, increment it, and then assign it to the second.

  • Actually, I think I will make my own next() function and define it before the code. Take that C++

OMG It Works :) - 3/30/16

• Ok so here is how it works

• The errors I was getting were because for some reason, Root doesn't like lists so we changed every list to a vector, just like how the original values were stored in a vector, except this vector has struct objects as elements. This works! Everything had to be changed to be compatible with vectors. They are accessed by the [] operator ie. "vector_name[index].member_variable" where member variable is one of the variables in the struct

• Because vectors can be accessed with [] there is no need for iterators so my methods and all that jazz with iterators can be deleted, but at least I know how to us them!

• One trick with vector.erase() it doesn't just take a number so you need to do vector_name.erase(vector_name.begin() + n) where n is the element you want to access.

• Also, compile using .L<filename>+ for some reason it helps

• I cleaned up my logic a little bit too...

  • Check the beam distance AFTER the jt (inner) loop because you only need one value

  • My Phiwrap Solution was also incorrect, you can't fix every value, there will be a discontinuity somewhere on the unit circle so I changed it to apply to the DELTAPHI value

  • Also, you need to reset the Dibmin, Dijmin, and minindex_i/j values after each run through of the it and jt loops ( but still inside the while loop) to make sure you're getting a value from the new list of particles. 

• I encountered a problem when removing the minindex_i and _j values. If the one you remove first is smaller than the second one then it will move the index for the other one back and throw everything off so we threw in an if statement to see which one is larger and then we remove that one

• Also, I'm not sure if I mentioned this before, but when you calculate the Dib, do the pt^-2 not pt^2 it is essential

• ***VERY IMPORTANT*** If for some reason the .h file starts being annoying and decides not to work (which happens a lot) it might be that the "std" thing is an issue. Things that aren't basic C++ code need special commands so the compiler recognizes them. If you type "using namespace std;" at the top it corrects every instance of this BUT you might also name something that is in this standard library on accident and so it is best to either write "std::" in front of things like vectors or cout's OR you can type "using std::vector" for example at the top to correct every instance of that specific keyword

Making Histograms- 3/30/16

• Ok so now that the code works, we need to learn a little about jets my plotting various histograms

• So far, I've only made one, the number of jets in each event, displayed below

• • Here's the pT  Spectrum of every jet. Most are very low energy as is evident by the plot below

• To easily view the saved files, use eog <filename

More Plots- 4/4/16

• so it took me a while but I realized that the plots were not accurate. The number of Jets should max out at 50... so I wasn't clearing the vector after every event which means the size just kept growing. I added a Jets.clear() and a HighEnergyJets.clear() and a SortingHEJets.clear() in there to fix this and the plots are much better

• I'm still trying to figure out where to make the cut for the High Energy Jets

• • These (^^) are the updated plots

Even More Plots - 4/10/16

• • Here is the plot that Dr. Eno asked us to make with the number of particles in every jet with a pT > 50 GeV

***For this week, compare down quark jets in PythiaTree.cc to dark quark jets and also compare to the jets we have in our code (but be wary of pT...)

Working on Dark Code - 4/19/16

• For this week, our task is to make plots of the dark jets from the pythiatree code and compare those to the plots of the regular jets that we already made. 

• We will also need to make plots of the jets from the Dark Event

• Things to note:

  • We are working in Pythia now so things can be accessed using the array [] operator and we will use Pythia's clustering algorithm

  • Particles have ID numbers as follows 

    • 4900 denotes a dark particle and v is added to name (Dv, gv, 

    • gluons = 21

    • d, u, s, c, b, t = {1, 2, 3, 4, 5, 6}

    • anti (bar) makes it negative

    • pi neutral = 111

    • rho neutral = 113

    • the charged versions of the pion and rho meson are 211, 213 and are positive or negative based on their charge

    • Creates multiple events - g g collison, q qbar both going to Dv and Dvbar

    • 22 (gamma) is a photon

• To run the code, do 

  • source PYTHIA_setup.csh

  • ./pythiaTree.exe >& output.txt

  • it takes a few minutes to run

• Make the changes in pythiaTree.cc and then root creates the .exe file for you and that is what gets run

• I have copied the original pythiaTree.cc to another file just in case I make any fatal mistakes and need to restart (lets hope that doesn't happen)

• aSlowJet and trigSlowJet are the jet clustering algorithms in Pythia

• Here is the link to the SlowJet class reference 

  • http://home.thep.lu.se/~torbjorn/doxygen/classPythia8_1_1SlowJet.html#ab99ae8739bea53f50e5c46e1e7fc720a

Generating Plots -4/25/16

• I have generated some plots... I think they're correct

• I'm using the variables dq1sj and d1sj which are index holders for the dark quarks and down quark jets in each event. I fill the histograms with the pT (for now) of those jets.

• I'm working on overlaying the two histograms using the ->Draw("same") command with limited luck

• To change the line color so "SetLineColor(n)" duh

• Working on the best format to view because it's kind of confusing how it is now with red and blue lines on white background

• Also, there is a dq2sj and a d2sj which i'm not sure what they are but it could give me more entries to help smooth out the plots, i might try adding them in and then just commenting out the line

Comparing Jets - 5/4/16

• As we expected the, pT and number of particles in the jets were similar (*there is one discrepancy due to, we believe, a bug in the code*) 

• What we think will be different is the impact parameter and the closest approach because the dark particles should be forming farther out in the detector

• • This is a diagram of the closest approach an 

Final Plots -5/4/16

• The semester is winding down and I'm currently writing the final paper for this whole research project

• Our final results confirmed our initial hypothesis that the dark quark jets would appear much later in the detector than the down quark jets

• To clarify, the down quarks basically appeared directly on the beam line whereas the dark quarks were created much farther out from the beam line, up to 4m before there are an insignificant amount left

• Here are the final plots

• • As you can clearly see, the closest approach for the down quarks is basically zero

• • The x axis for the plot below is in mm