Embedded Files
References:
[1] http://pdg.lbl.gov/2007/reviews/montecarlorpp.pdf - particle ID codes for pythia
[2] http://lcgapp.cern.ch/project/docs/pythia8080.pdf - intro to pythia, some data on commands
[3] https://root.cern.ch/root/roottalk/roottalk02/att-1402/01-h200.C - could be helpful in writing code to find highest pT events
[4] http://stackoverflow.com/questions/18723007/mass-invarint-of-bbbar-pair-by-pythia -plot invarient mass
[5] http://home.thep.lu.se/~torbjorn/talks/tutorial81.pdf
[6] https://root.cern.ch/root/HowtoHistogram.html -histograms in root.
[7] https://root.cern.ch/root/html/tutorials/hist/hsum.C.html -how to plot multiple histograms in root
[8] https://root.cern.ch/phpBB3/viewtopic.php?t=7879 -how to plot multiple histograms in root
[9] https://root.cern.ch/doc/v606/classTHistPainter.html -histogram class root
2/11/2016:
Continued from 2/10/2016, we were working past midnight.
The cross section for Z' 7.33e-11 mb, meaning that in run one of the LHC we produce 1466.
This means our scaling factor (the number we will multiply the DY events to normalize it to) is the ratio: (# of DY events)/(#of Z' events) which is: 35470.6684857
Our ratio of signal to Background should be 1/ (# of DY events)/(#of Z' events) = 0.00004829.
The histogram for Z' is:
The histogram for DY (un-scaled) is:
Working with Tal: To make these 2 histograms in the same file we will make a pythia (Pythia pythia;) we will also make a pythia1 (Pythia pythia1;)
the pythia runs for DY and pythia1 runs for Z'. Basically copy all the code but whenever it says pythia type pythia1.
To plot the 2 histograms on the same TCanvas, you just type t1->Draw(); t2->Draw("same");
This Will plot the 2 histograms on the same canvas.
In order to scale the DY histogram up, we just type mass->Fill(pT2 +pT1,scale factor);
Where our scale factor is 35470. This means that insted of increasing the bin count by 1 it increases it by the scale factor.
If we cut on 250 GeV, the background clearly goes away and we are really only left with signal, this is because our background happens with a mass of 92GeV, while our signal has a mass of 1500 GeV so the mass peaks are so separated.
In order for me to continue to use Tal's code after he updated it to git hub and remove my changes type:
git-checkout pythiaTree.cc
git pull origin
To make the histogram draw in the root file just make a TCanvas right under where root makes the root file. Then after we draw the 2 histograms just do c1->Write();
Before any cuts the histogram of background vs. signal is:
(the signal is in grey and the background is in yellow).
After Cutting on 250GeV:
As you can see we have 100% background rejection, so our cut at 250 is infinity better than no cuts, this is becuse the masses are so different they don't even overlap.
2/10/2016:
To get the transverse momentum for only electrons in the: if (pythia.eveent[i].isFinal()) loop ive inserted:
if(pythia.event[i].id()==11){
Pt[nTot-1]=pythia.event[i].pT();
}
11 is the identifier code for electrons, so this loop selects for only stable electrons. For these stable electrons the transverse momentum is loaded into pT.
Tal and Avi have helped me to understand how to plot the mass of the highest energy pT events.
You also need to look at particles with id -11, which are positrons.
I have taken Tal's code from his github https://github.com/robotal/PythiaWork
Basically in his code, in the i for loop, Tal says if the ith particle in the ievent is and electron, if so if it is the highest pT event, make the old highest event the 2nd highest event, and set the new one as the highest. If the event is lower than the highest but higher than the 2nd highest transfer the current one into the 2nd highest one.
Then outside the i loop but still inside the ievent loop Tal fills the mass histogram with the sum of the 2 highest energy events. This is looped over for every ievent, so we will have 1000 events in out histogram.
I find a sigma of 2.605e-6 Barns, this sigma is at the botom of the output.txt file. There is a value of 330 mili barns called sigma interaction, but this number seemed too high for DY, so we decided to use the 2.605e-6 number.
N_data = 2.605 *10^6 fb * 20fb^-1 = 5.2*10^7 = which means that in run 1 of the LHC we produced 50 million DY events that decayed to electrons.
2/4/2016:
For the Code, the masses are wrong, pythia uses units of GeV, so we need to change 150. to 1500. , 90. to 900. and 300. to 3000.
The new code should be:
Let's make a Z' which decays to electron pairs and has a mass of 1500 GeV. In order to do this, remove the lines that tell it to make QCD events and add instead
pythia.readString("NewGaugeBoson:ffbar2gmZZprime = on");
pythia.readString("32:m0 = 1500.");
pythia.readString("32:mMin = 900.");
pythia.readString("32:mMax = 3000.");
pythia.readString("PhaseSpace:mHatMin = 20. ");
pythia.readString("Zprime:gmZmode = 3");
pythia.readString("Beams:eCM = 14000.");
pythia.readString("32:onMode = off");
pythia.readString("32:onIfAny = 11");
The new pT plot, before any cuts looks like:
2/2/2016:
Code from 1/31/2016 now works. I removed a line at the end that said t1->Draw();
DO NOT ADD T1->DRAW():
Here is my mass histogram in root:
I added the code in the tutorial to change the beam to produce Z' which decay to e- pairs. The new mass graph looks like:
Using the (http://home.thep.lu.se/~torbjorn/pythia82html/Welcome.html) as a reference I found out what all the commands meant:
pythia.readString("NewGaugeBoson:ffbar2gmZZprime = on");This line turns on this flag, which is default off. Turns on f fbar scattering to z'0 channel.
pythia.readString("32:m0 = 150."); Note, 32 means Z'0 boson. Sets the center of the mass dist. at 150
pythia.readString("32:mMin = 90."); Sets the minimum of the mass dist. at 90.
pythia.readString("32:mMax = 300."); Sets the max of the mass dist. to 300.
pythia.readString("PhaseSpace:mHatMin = 20. "); Sets the minimum invariant mass to 20.
pythia.readString("Zprime:gmZmode = 3");Sets the ffbar to Z'0 to only include pure Z'0 process, no gammas included.
pythia.readString("Beams:eCM = 14000."); Sets 14TeV center of mass energy for the beams [1]
pythia.readString("32:onMode = off"); Sets the z'0 channel to off [1]
pythia.readString("32:onIfAny = 11"); Switches back on the Z'0 channel for electrons [1] [2]
In order to plot the invarient mass of the highest 2 pT events we create a class similar to old homework #9:
root hist.root
t1->MakeClass("t1") //creates MyClass.C and MyClass.h files that we can add code in to do our analysis.
1/31/2016:
Pythia is a monte carlo software specifically designed for beam physics. You specify the type of beam particle and its center of mass energy and pythia will generate a list of the particles produced and their properties. This can be useful when optomizing cuts. Useful info and commands in pythia can be found here: http://home.thep.lu.se/~torbjorn/pythia82html/Welcome.html
In order to run Sarah's pythia code we must first clone it from her git repository. To clone (copy git code into a newly created folder) her code we type: git clone https://github.com/saraheno/honrpythia.git
This creates a directory called honrpythia, we then type: cd honrpythia
source PYTHIA_setup.csh
make pythiaTree
./pythiaTree >& output.txt
This makes the pythiaTree, runs the code, and puts it in the output.txt file. Shabnams notes: This program uses pythia to make dijet events (quark-quark, quark-gluon, or gluon-gluon final states), requring the transverse momenta of the jets be greater than 20 GeV, in proton-proton collisions with a center of mass energy of 14 TeV. Look at pythiaTree.cc to see where these parameters are set in the code. The program also makes two histograms and a tree.
From the pythia website linked above (in study output, particle properties) I have found this: double Particle::mSel()
the mass of the particle, picked according to a Breit-Wigner distribution for particles with width. It is different each time called, and is therefore only used once per particle to set its mass m().
This suggests that to add the mass to the tree I have to type code like the momentum vectors (px, py, pz) but replace .px() with .m().
I have inserted the following lines of new code (in green) among the old code (in red):
Float_t Apz[5000];
Float_t m[5000];
and,
t1->Branch("Apz",Apx,"Apz[nTot]/F");
t1->Branch("Mass",Apx,"mass[ntot]/F");
and,
Apz[nTot-1]=pythia.event[i].pz();
m[nTot-1]=pythia.event[i].m();
In order to check my work I loaded root with the tree: root hist.root
and once root loaded typed: TBrowser New, and navagated throught to hist.root, t1 (where I stored the mass) and it wasn't there so clearly this code does not work at this point.
To make the code run 10,000 insted of 100 particles I modified the line to say:
for (int iEvent = 0; iEvent <10000; ++iEvent) {
________________________________________________________________________________
11/15/2015:
The toyMC.C monte carlo code works by creating a root tree where the branches are diffrent combonations of momenta of the daughter particles. It randomly generates, a boson, simulates its decay, and smears the decay numbers so they would look like they came from a real detector.
You can write a new class in root with the MakeClass function. Ie. MakeClass("MyClass") this class would be defined by the contents of MyClass.C and MyClass.h, you can add regular root C code to MyClass.C so that the class will run the code.
To actually run the class type .L MyClass.C+
MyClass m
m.Loop() <- this loops over all the root files attached to this class, this is where root comes in so handy because you dont have to worry about file organization.
to do the higgs analysis:
TChain* chain = new TChain("HZZ4LeptonsAnalysisReduced"); // this needs to match the TTree name
this command make a chain of files that are going to be filled with all the relevant collision for Z analysis.
chain->Add("/data/users/jabeen/HONR268N_HiggsAnalysis/Excercise-Electron/output_DoubleElectron_MC.root");
This command actually adds the files to the chain. Note: chained files must have the same structure.
Inside the ana.C class code : for (Long64_t jentry=0; jentry < nentries; jentry++) {
loop is where we can plot histograms and make cuts and such.
In ana.h are the names of all sorts of variables that we can use for our analysis. ie. f_lept1_eta is the pesudorapidity of the most energetic lepton.
The TLorentzVector is a way in root of analyzing any 4 vector system, it takes doubles as its argument. Ie. TLorentzVector el1, el2;
TLorentzVector el1, el2;
el1.SetPtEtaPhiM(f_lept1_pt, f_lept1_eta, f_lept1_phi, 0.0); el2.SetPtEtaPhiM(f_lept2_pt, f_lept2_eta, f_lept2_phi, 0.0);
will create 2 lorenz vectors.
if we add these 2 vectors and plot a histogram using our monte carlo for Z->2e- we get:
As we can see this peaks at the mass of the Z boson.
if in the ana.C in the if (fchain==0) statement we add:
TFile* output = TFile::Open("Dielectron_MC.root", "RECREATE"); // "RECREATE" would produce
this creates a new file, the subsiquent code:
el1_eta->Write(); el2_eta->Write(); el1_pt->Write(); el2_pt->Write(); Z_ee->Write(); output->Close();
Will write these diffrent histograms to the output file, making it easier for us to look at all the graphs.
Once that code is written to see the graphs open the root file: root -l DIelectron_MC.root
el1_eta->Draw()
this will draw the histogram, but it will not unless you have opened the Dielectron_MC.root file first, it no longer works just to have run the ana.C class as stated above, you need to have opened the root file.
To make cuts is just a matter of adding simple if statements. To cut on Pt for example in the for loop with jentries where you fill the histograms put the filling of the histograms inside the following if statement:
if (f_lept1_pt > 60) {
this will cut on tranverse momentum is bigger than 60.
Keep in mind when modifying ana,C the tutorial left out the line: Long64_t ientry = LoadTree(jentry); //the code needs this line to run
11/11/2015:
A tree is a way of organizing files in root. It allows for data to be accessed easily while still maintaining all of the information for every collision.
Withing root if you type TBrowser New it opens up a browser window based GUI that lets you manually click on diffrent leafs in the tree and look at plots for the dirrfernt things in the tree, such as Px, Py, or Pz.
To graph a function in root use the TF1 class (1 variable function graph), it takes the parameters, name, function, lower bound, upper bound. For example to plot x^2 from 0 to 10
{
TF1 *fa1 = new TF1("fa1","pow(x,2)",0,10);
fa1->Draw();
}
To draw error bars on a dot graph in root create a TCanvas, then create a 1D array with X values, a 1D array with the error values, and do the same for y, such that you have 4 arrays. Now create a new TGraphErrors which takes the inputs (number of points, x array, y array, errors x array, errors y array) for example:
const Int_t n=10;
Float_t x[n] = {1, 2, 3, 4, 5, 6};
Float_t y[n] = {3.9, 6.5, 8.1, 13.5, 15.5, 16};
Float_t ex[n] = {0.01, 0.01, 0.01, 0.01, 0.01, 0.01};
Float_t ey[n] = {0.5, 0.5, 0.5, 0.5, 0.5, 0.5};
TGraphErrors *gr = new TGraphErrors(n,x,y,ex,ey);
11/2/2015:
Now using moba xterm insted of putty, better for root
In root c code, TH1F is a custom class (type of variable) that is a histogram: https://root.cern.ch/doc/master/classTH1F.html
TH1F is a special case of the TArrayF class, it is everything that TArrayF is, but with special extra features
Creators are necessary for classes, they detrimine what kinds of data the class can use as arguments and where it is stored.
The syntax for the TH1F class is: TH1F("Name", "Title",# of bins, Lower bin boundry, upper bin boundry)
for example: TH1F *histo1 = new TH1F("histo1","S",50,0.,50.); //50 bins from 0 to 50
Member functions are things you can do with a class, ie, draw your histogram
The 3 lines necessary to create a histogram (taken from HW7 example code:
TH1F *histo1 = new TH1F("histo1","s",50,0.,50.); //50 bins from 0 to 50
histo1->SetMarkerStyle(21); // sets the "look" of the histogram
histo1->Draw(""); //displays the histogram
The code:
TRandom r; //makes random number r
gRandom->SetSeed(); //sets the random seed <--- the number that is the basis for the random # generator
float p1=r.Gaus(10,2); // gaussian distrubution mean of 10, standard devation (sigma) of 2
if in the example code for HW7 i replace float p1=r.Gaus(5.,1.); with float p1=r.Poisson(5); the fit is always very poor, this is because they are not the same distribution.
10/22/2015:
ROOT is an envoirment that has a bunch of useful built in c and python libraries.
Before running root go to ssh option then x11 then enable x11 forwarding before logging in using putty.
To run root type:
cd to the CMSSW_5_3_30 directory
cmsenv
cd testRoot
root -nw
in root every command is in the .command structure
ie. .x testfile.c runs the testfile.c code, in root .quit quits root.
In C++ an array is a variable with may values, to make an array of size 6 for example type int j[6] = {0,1,2,3,4,5}
this will store 0-5 in each of j's six places
The indexes for all arrays begin at 0.
You can make a 2D array like this: int array [1][1] = {0,1,2,3}
now 0-3 are stored based on 2 pieces of info, not just 1 for a regular array, for ex. array[0][1] =1 array [1][0]= 2
To read data from a file do #include<fstream>
ofstream myfile this initilized the command my file, but I could call it whatever I liked
myfile.open("filename") open the file filename
myfile <<"stuff"" inserts the word stuff into the file filename
myfile.close(); this closes the file
If you want to use multiple .cpp files for one piece of code you can write subroutines in different files and include them in each other to run it as one big file. To include another file in your's
#include "file.cpp"
now i can call all subroutines in file.cpp just as if I wrote them in my code.
To generate random numbers do:
#include <stdio.h>
#include <stdlib.h> // this includes commands like rand and srand
#include <time.h>
srand(time(NULL)): // seeds random number with time which is always unique
rand()%100 //this produces a random number 0-99
You can also make classes in C++ which is why it is called object oriented. A class is a variable type much like an int or a double but it can be much more complex. For example when programming a video game you can have a boat class which has subclasses motor and sailboats. Those can both hold info Ie. # of people.
10/8/2015:
Data (contents of variables) in cpp is stored in a memory adress. The location of the adress for int i is &i, the and character refrences a memory location. A pointer is a variable that is the location of another varible. A pointer is set by int*
To refrence the data stored at the location of a pointer, for ex. int* p, type: *p
You can also create a value without storing it in a memory location by new int(5), this can be used to store data in the memory adresss of a specific pointer without using a variable to store that value ie. int* p = new int(5).
Now *p =5 stored in memory adress p.
The code I made for HW5, the last part of the tutorial is:
#include <iostream>
using namespace std;
int main () {
double* p1 = new double(3.14);
double* p2 = p1;
cout << "the adress of *p1 is " <<p1 << "the value of *p1 is "<< *p1<< endl;
cout << "the adress of *p2 is "<< p2 << "the value if *p2 is " <<*p2<<endl;
*p1 = *p1 *2;
cout <<"*p1 is " <<*p1<<endl;
cout <<"*p2 is "<<*p2<<endl;
return 0;
}
10/1/2015:
I used this proof (p.7-8) http://www.hep.shef.ac.uk/edaw/PHY206/Site/2012_course_files/phy206rlec7.pdf
to prove that pseudo rapidity is equal to rapidity for massless particles. This proof relies on the (1-x)^-1/2 taylor series expansion, but I found a way to do it without the expansion.
Notes for HW4 tutorial:
The lines #include<iostream> and using namespace std; include libraries and commands common to all c++ users
The result of a program multiplying an integer with a double, and storing the result in an integer is the result will be truncated. ie. 2*3.3 =6.6 but will display as 6.
++ adds one to the variable and -- subtracts it. a++; essentially take the place of the line a = a +1;
While loops will only do the loop while the argument (boolean variable) is true (1).
For loops will loop until the variable declared in the first part meets the criteria of the second part, and the loop counts up or down depending on the third part. ie. for (int a=0; a<6 ++) { will go from a=0 to a=6, hitting every integer inbetween);
I rewrote the code using for loops:
#include <iostream>
using namespace std;
int main() {
int idebug =0;
for (int i=0; i<10; i++) {
cout<< "i is " <<i<<": "<<endl;
int m=0;
int idebug =1;
for (m; m<=i; m++) {
cout <<m;
}
cout << endl;
}
return 0;
}
9/29/2015:
to make c++ code run on the shell write the code in emacs, ex. emacs -nw ctest.cpp
the .cpp means it is a C++ code file. To fun C++ code type g++ filename.cpp
./a.out
For the var.cpp script (under variables in HW4) an integer defined as an integer times a double will truncate (not round) the result
++ add 1 and -- subtracts 1.
for loops take syntax for (a variable; boundry constraints; ++ or --) {
}
<< means send data to the left >> means send it to the right
the cout command prints stuff to the terminal
ex. cout <<"hello world<<endl; // send the endl command (end the line in a break) to the text hello world which then gets sent to cout which prints it
to add files to git hub type git add filename.extension
git commit -a -m "comments"
git push -u origin master
gdb is a free debugger for c++ files, to use with a file type g++ -g filename.cpp
gdb (this runs gdb)
file a.out (must be a.out and not filename.cpp)
break 12 (inserts a code break at line 12)
run (runs code until predifined code break)
info locals (prints variable values at the break)
OTHER GDB commands:
q (quit)
step (runs just the next line)
list (list the first 10 lines)
list 12 (list lines centered arround line 12)
integers in c++ have a max value ~2e9
9/22/2015:
Relativity notes:
9/21/2015:
The proper structure to add new files to a git repository is:
1. git init
2. git add filename.extentsion
2. git commit -a -m "your comments here"
3. git remote add origin git@github.com/yourusername/repositoryname.git (< only need to do this the first time you push)
4. git push -u origin master
lots of other helpful tutorials on the github website
9/20/2015:
Fixed ssh-add error, needed to use ` instead of '. Finished github tutorial. successfully pushed README.md from directory gittest to github website using the command git push -u origin master.
9/19/2015:
In order to give yourself permission to execute a file, not just read or write, you need to do chmod +x (gives everyone permission to run)
the > symbol sends thing to different files.
the symbol | is used for piping and sends data between 2 different programs with the syntax: program1 | program2. this sends data from program 1 to program 2.
Did condor tutorial linked in tutorial 2.
Could not add ssh key, thus could not finish github tutorial.
09/17/2015:
Notes on tutorial 2:
Wrote code using emacs (testscript.tcsh):
#!/bin/tcsh
#here we can put some comments
set ARG=$1
#here as a bug statement we echo out the value to the screen
echo $ARG
#here we look in some files
find /usr -name "*$ARG*" -print
this code finds files with the argument in the name. Call script by typing tcsh ./testscript.tcsh argument
typed chmod +x testscript.tcsh and now the file runs by typing: ./testscript.tcsh
the chmod command controlls who has access to run the file. chmod +x is the same command as chmod a+x which lets anybody run the file. to let just the user run the file type chmod u+x
created github account jacobsiegel22
Page updated
Report abuse