diff --git a/.gitignore b/.gitignore
index 0c26e82..fea4887 100644
--- a/.gitignore
+++ b/.gitignore
@@ -66,3 +66,10 @@
 *.out
 *.app
 
+
+# my stuff
+**/job*
+**/*.gz 
+**/*.tgz
+**/*~
+**/*.root
diff --git a/carbon/MPVcorrection_carbon0mm05.txt b/carbon/MPVcorrection_carbon0mm05.txt
new file mode 100644
index 0000000..dc5574d
--- /dev/null
+++ b/carbon/MPVcorrection_carbon0mm05.txt
@@ -0,0 +1,26 @@
+0 0.987199
+1 0.981107
+2 0.976363
+3 0.971173
+4 0.96674
+5 0.962492
+6 0.960485
+7 0.956097
+8 0.9526
+9 0.95099
+10 0.947653
+11 0.945952
+12 0.943821
+13 0.941682
+14 0.940207
+15 0.939598
+16 0.937827
+17 0.935669
+18 0.934759
+19 0.933702
+20 0.932179
+21 0.930267
+22 0.931688
+23 0.930076
+24 0.928356
+25 0.927973
diff --git a/carbon/Makefile b/carbon/Makefile
new file mode 100644
index 0000000..c16461b
--- /dev/null
+++ b/carbon/Makefile
@@ -0,0 +1,38 @@
+
+ROOTCFLAGS   := $(shell root-config --cflags)
+ROOTLIBS     := $(shell root-config --libs)
+ROOTGLIBS    := $(shell root-config --glibs)
+
+GSLCFLAGS   := $(shell gsl-config --cflags)
+GSLLIBS     := $(shell gsl-config --libs)
+GSLGLIBS    := $(shell gsl-config --glibs)
+
+
+LDFLAGS       = -O
+LIBS         += $(ROOTLIBS) $(GSLLIBS) -L/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/lib
+CFLAGS       += $(ROOTCFLAGS) $(GSLCFLAGS)  -I/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/include
+
+#  Not sure why Minuit isn't being included -- put in by hand
+#
+LIBS         += -lMinuit -ltbb
+
+
+all: eventdatroot
+
+
+
+eventdatroot: eventdatroot.o
+	g++ -o eventdatroot  eventdatroot.o $(LDFLAGS) $(LIBS)
+
+%.o: %.c
+	g++ ${CFLAGS} -c  -g $<
+
+
+test:
+	@echo $(ROOTCFLAGS) 
+	@echo $(LDFLAGS) 
+	@echo $(LIBS)
+
+clean:	
+	-rm -f  *.o
+
diff --git a/carbon/carbon0mm5_SPratio.pdf b/carbon/carbon0mm5_SPratio.pdf
new file mode 100644
index 0000000..03a5f9c
Binary files /dev/null and b/carbon/carbon0mm5_SPratio.pdf differ
diff --git a/carbon/eventdatroot b/carbon/eventdatroot
new file mode 100755
index 0000000..4676f31
Binary files /dev/null and b/carbon/eventdatroot differ
diff --git a/carbon/eventdatroot.c b/carbon/eventdatroot.c
new file mode 100644
index 0000000..4bad079
--- /dev/null
+++ b/carbon/eventdatroot.c
@@ -0,0 +1,137 @@
+#define eventdatroot_cxx
+#include "eventdatroot.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+//#include <TROOT.h>
+#include <TF1.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+
+
+using namespace std;
+
+void eventdatroot::Loop(Float_t scale)
+{
+  //   In a ROOT session, you can do:
+  //      root> .L eventdatroot.C
+  //      root> eventdatroot t
+  //      root> t.GetEntry(12); // Fill t data members with entry number 12
+  //      root> t.Show();       // Show values of entry 12
+  //      root> t.Show(16);     // Read and show values of entry 16
+  //      root> t.Loop();       // Loop on all entries
+  //
+
+  //     This is the loop skeleton where:
+  //    jentry is the global entry number in the chain
+  //    ientry is the entry number in the current Tree
+  //  Note that the argument to GetEntry must be:
+  //    jentry for TChain::GetEntry
+  //    ientry for TTree::GetEntry and TBranch::GetEntry
+  //
+  //       To read only selected branches, Insert statements like:
+  // METHOD1:
+  //    fChain->SetBranchStatus("*",0);  // disable all branches
+  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
+  // METHOD2: replace line
+  //    fChain->GetEntry(jentry);       //read all branches
+  //by  b_branchname->GetEntry(ientry); //read only this branch
+
+  Float_t ratio_electron, ratio_primary, etotal, spratio;
+  // TF1 * f_sp_ps = new TF1("f_sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in polystyrene [MeV cm2 /g]
+  // f_sp_ps->SetParameters(361.936, -0.892255, 1.19133);  //protons between 50 and 220 MeV
+  TF1 * f_c12sp_ps = new TF1("f_c12sp_ps","[0]*pow(x,[1])+[2]", 80, 480); //stopping power of carbon in polystyrene [MeV cm2 /g]
+  f_c12sp_ps->SetParameters(14538.9, -0.922423,49.2859);  //carbon between 80 and 480 MeV
+
+  
+  if (fChain == 0) return;
+
+  Long64_t nentries = fChain->GetEntriesFast();
+
+  Long64_t nbytes = 0, nb = 0;
+  for (Long64_t jentry=0; jentry<nentries;jentry++) {
+    Long64_t ientry = LoadTree(jentry);
+    if (ientry < 0) break;
+    nb = fChain->GetEntry(jentry);   nbytes += nb;
+    // if (Cut(ientry) < 0) continue;
+    //    cout << jentry << endl;
+
+    h_endist_0->Fill(ENDIST[0]); // edep by ionisation
+    h_endist_1->Fill(ENDIST[1]); // edep by pi0, e-, e+ and #gamma
+    h_endist_2->Fill(ENDIST[2]); // edep by nuclear recoils and heavy fragments
+    h_endist_3->Fill(ENDIST[3]); // edep by part. #lt threshold
+    h_endist_4->Fill(ENDIST[4]); // E leaving the syste
+    h_endist_5->Fill(ENDIST[5]); // E carried by discarded particles
+    h_endist_6->Fill(ENDIST[6]); // resid, excit. E after evap.
+    h_endist_7->Fill(ENDIST[7]); // edep by low-energy neutrons
+    h_endist_8->Fill(ENDIST[8]); // E of part. #gt time limit
+    h_endist_9->Fill(ENDIST[9]); // E lost in endothermic nuclear reactions
+    h_endist_10->Fill(ENDIST[10]); // E lost in endothermic low-E n-reaction
+    h_endist_11->Fill(ENDIST[11]); // missing E
+
+    etotal = 0.;
+    for (unsigned int i = 0; i<3; i++ ){
+      etotal += ENDIST[i];
+    }
+
+    ratio_electron = ENDIST[1] / etotal;
+    ratio_primary = ENDIST[0] / etotal;
+
+    h_ratio_e->Fill(ratio_electron);
+    h_ratio_p->Fill(ratio_primary);
+    h_sum_ep->Fill(ENDIST[1]+ENDIST[0]);
+
+    h_let_ep->Fill( (ENDIST[1]+ENDIST[0])/scale/36. ); // [GeV/cm/Z^2]
+    h_let_e->Fill( (ENDIST[1])/scale/36. );
+    h_let_p->Fill( (ENDIST[0])/scale/36. );
+
+    spratio = ( ( (ENDIST[0]+ENDIST[1])/scale*1000./1.06) /  f_c12sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000./12. ) ); //ratio of deposited energy to PSTAR stopping power (MeV/cm/density) / MeVcm2/g
+    h_spratio->Fill(spratio);
+
+    //   cout <<  ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06)  << " " << f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) << " " << spratio << endl;
+
+  }
+}
+
+
+
+int main(int argc, const char* argv[])
+{
+  Int_t thickkey = atoi(argv[2]);
+  Float_t THICKNESS[6] = {0.025, 0.050, 0.100, 0.150, 0.500, 1.000}; //#cm
+  // Float_t THICKNESS[5] = {0.050, 0.100, 0.150, 1.000, 10.00}; //#cm
+
+  ///open the input and output root files
+  char  finname[50];
+  sprintf(finname, "%s.root",argv[1]);
+  TFile * f = new TFile(finname,"OPEN");
+
+  char  foutname[50];
+  sprintf(foutname, "%s_out.root",argv[1]);
+  cout << foutname << endl;
+  TFile * fout = new TFile(foutname, "RECREATE");
+
+  //run the analysis loop above
+  eventdatroot t;
+  TTree * tree;
+  f->GetObject("EVENTDAT",tree);  
+  t.Init(tree);
+  t.Loop( THICKNESS[thickkey] );
+  t.Closefile(fout);
+
+  return 1;
+
+}
diff --git a/carbon/eventdatroot.h b/carbon/eventdatroot.h
new file mode 100644
index 0000000..f4cb35a
--- /dev/null
+++ b/carbon/eventdatroot.h
@@ -0,0 +1,199 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Mar  8 13:05:30 2019 by ROOT version 6.12/06
+// from TTree EVENTDAT/FLUKA Course Exercise	DATE:  3/ 4/19,  TIME: 11:54:46
+// found on file: jobs/runjob1001_eventdata.root
+//////////////////////////////////////////////////////////
+
+#ifndef eventdatroot_h
+#define eventdatroot_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1.h>
+
+// Header file for the classes stored in the TTree if any.
+
+class eventdatroot {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Float_t         DATA_ALL_PART[6];
+   Float_t         DATA_BEAMPART[6];
+   Float_t         DATA_ENERGY[6];
+   Float_t         DATA_id211[6];
+   Int_t           SEEDS[4];
+   Float_t         ENDIST[12];
+
+   // List of branches
+   TBranch        *b_DATA;   //!
+   TBranch        *b_seed;   //!
+   TBranch        *b_endist;   //!
+
+      
+   TH1F * h_endist_0;
+   TH1F * h_endist_1;
+   TH1F * h_endist_2;
+   TH1F * h_endist_3;
+   TH1F * h_endist_4;
+   TH1F * h_endist_5;
+   TH1F * h_endist_6;
+   TH1F * h_endist_7;
+   TH1F * h_endist_8;
+   TH1F * h_endist_9;
+   TH1F * h_endist_10;
+   TH1F * h_endist_11;
+
+   TH1F * h_ratio_e;
+   TH1F * h_ratio_p;
+   TH1F * h_sum_ep;
+   
+   TH1F * h_let_ep;
+   TH1F * h_let_e;
+   TH1F * h_let_p;
+
+   TH1F * h_spratio;
+
+   eventdatroot(TTree *tree=0);
+   virtual ~eventdatroot();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop(Float_t scale);
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+   virtual Int_t    Closefile(TFile * file);
+};
+
+#endif
+
+#ifdef eventdatroot_cxx
+eventdatroot::eventdatroot(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+  //  if (tree == 0) {
+     //   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("jobs/runjob1001_eventdata.root");
+      //     if (!f || !f->IsOpen()) {
+		//      f = new TFile("jobs/runjob1001_eventdata.root");
+	 //     }
+      //    f->GetObject("EVENTDAT",tree);
+
+      //   }
+   //   Init(tree);
+}
+
+eventdatroot::~eventdatroot()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t eventdatroot::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t eventdatroot::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void eventdatroot::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("DATA", DATA_ALL_PART, &b_DATA);
+   fChain->SetBranchAddress("SEEDS", SEEDS, &b_seed);
+   fChain->SetBranchAddress("ENDIST", ENDIST, &b_endist);
+
+   h_endist_0 = new TH1F("h_endist_0","edep by ionisation",400,0.0,0.05);
+   h_endist_1 = new TH1F("h_endist_1","edep by pi0, e-, e+ and #gamma",400,0.0,0.03);
+   h_endist_2 = new TH1F("h_endist_2","edep by nuclear recoils and heavy fragments",400,0.0,0.001);
+   h_endist_3 = new TH1F("h_endist_3","edep by part. #lt threshold",200,0.0,0.001);
+   h_endist_4 = new TH1F("h_endist_4","E leaving the system",200,0.0,1.1);
+   h_endist_5 = new TH1F("h_endist_5","E carried by discarded particles",200,0.0,0.1);
+   h_endist_6 = new TH1F("h_endist_6","resid, excit. E after evap.",200,0.0,0.1);
+   h_endist_7 = new TH1F("h_endist_7","edep by low-energy neutrons",200,0.0,0.1);
+   h_endist_8 = new TH1F("h_endist_8","E of part. #gt time limit",200,0.0,0.1);
+   h_endist_9 = new TH1F("h_endist_9","E lost in endothermic nuclear reactions",200,0.0,0.1);
+   h_endist_10 = new TH1F("h_endist_10","E lost in endothermic low-E n-reactions",200,0.0,0.1);
+   h_endist_11 = new TH1F("h_endist_11","missing E",400,0.0,0.01);
+
+   h_ratio_e = new TH1F("h_ratio_e","fraction of edep by pi0, e-, e+ and #gamma  ",400,0.0,1.0);
+   h_ratio_p = new TH1F("h_ratio_p","fraction of edep by ionisation",400,0.0,1.0);
+   h_sum_ep = new TH1F("h_sum_ep","sum of edep by ion, pi0, e-, e+ and #gamma,",400,0.0,0.1);
+
+   h_let_ep = new TH1F("h_let_ep","LET/Z^{2} (GeV/cm) by ion, pi0, e-, e+ and #gamma,",200,0.0,0.05);
+   h_let_e = new TH1F("h_let_e","LET/Z^{2} (GeV/cm) by pi0, e-, e+ and #gamma,",200,0.0,0.02);
+   h_let_p = new TH1F("h_let_p","LET/Z^{2} (GeV/cm) by ion",200,0.0,0.02);
+
+   h_spratio = new TH1F("h_spratio","total LET/SP (MeV/cm/#rho)/SP(MeVcm2/g)",200,0.5,3.0);
+
+   
+   Notify();
+}
+
+Bool_t eventdatroot::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void eventdatroot::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t eventdatroot::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+
+Int_t eventdatroot::Closefile(TFile * file)
+{
+  file->cd();
+  file->Write();
+  file->Close();
+
+  return 1;
+}
+
+#endif // #ifdef eventdatroot_cxx
diff --git a/carbon/langaus.C b/carbon/langaus.C
new file mode 100644
index 0000000..55c9fd2
--- /dev/null
+++ b/carbon/langaus.C
@@ -0,0 +1,348 @@
+ /// \ingroup tutorial_fit
+ /// \notebook
+ /// Convoluted Landau and Gaussian Fitting Function
+ ///         (using ROOT's Landau and Gauss functions)
+ ///
+ ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
+ ///
+ ///  to execute this example, do:
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C
+ /// ~~~
+ ///
+ /// or
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C++
+ /// ~~~
+ ///
+ /// \macro_image
+ /// \macro_output
+ /// \macro_code
+ ///
+ /// \authors H.Pernegger, Markus Friedl
+ 
+ #include "TH1.h"
+ #include "TF1.h"
+ #include "TROOT.h"
+ #include "TStyle.h"
+ #include "TMath.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+ 
+ Double_t langaufun(Double_t *x, Double_t *par) {
+ 
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+ 
+       // Numeric constants
+       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+       Double_t mpshift  = -0.22278298;       // Landau maximum location
+ 
+       // Control constants
+       Double_t np = 200.0;      // number of convolution steps
+       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
+ 
+       // Variables
+       Double_t xx;
+       Double_t mpc;
+       Double_t fland;
+       Double_t sum = 0.0;
+       Double_t xlow,xupp;
+       Double_t step;
+       Double_t i;
+ 
+ 
+       // MP shift correction
+       mpc = par[1] - mpshift * par[0];
+ 
+       // Range of convolution integral
+       xlow = x[0] - sc * par[3];
+       xupp = x[0] + 2*sc * par[3];
+ 
+       step = (xupp-xlow) / np;
+ 
+       // Convolution integral of Landau and Gaussian by sum
+       for(i=1.0; i<=np/2; i++) {
+          xx = xlow + (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+ 
+          xx = xupp - (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+       }
+ 
+       return (par[2] * step * sum * invsq2pi / par[3]);
+ }
+ 
+ 
+ 
+ TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
+ {
+    // Once again, here are the Landau * Gaussian parameters:
+    //   par[0]=Width (scale) parameter of Landau density
+    //   par[1]=Most Probable (MP, location) parameter of Landau density
+    //   par[2]=Total area (integral -inf to inf, normalization constant)
+    //   par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    // Variables for langaufit call:
+    //   his             histogram to fit
+    //   fitrange[2]     lo and hi boundaries of fit range
+    //   startvalues[4]  reasonable start values for the fit
+    //   parlimitslo[4]  lower parameter limits
+    //   parlimitshi[4]  upper parameter limits
+    //   fitparams[4]    returns the final fit parameters
+    //   fiterrors[4]    returns the final fit errors
+    //   ChiSqr          returns the chi square
+    //   NDF             returns ndf
+ 
+    Int_t i;
+    Char_t FunName[100];
+ 
+    sprintf(FunName,"Fitfcn_%s",his->GetName());
+ 
+    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
+    if (ffitold) delete ffitold;
+ 
+    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
+    ffit->SetParameters(startvalues);
+    ffit->SetParNames("Width","MP","Area","GSigma");
+ 
+    for (i=0; i<4; i++) {
+       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
+    }
+ 
+    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
+ 
+    ffit->GetParameters(fitparams);    // obtain fit parameters
+    for (i=0; i<4; i++) {
+       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
+    }
+    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
+    NDF[0] = ffit->GetNDF();           // obtain ndf
+ 
+    return (ffit);              // return fit function
+ 
+ }
+ 
+ 
+ Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
+ 
+    // Seaches for the location (x value) at the maximum of the
+    // Landau-Gaussian convolute and its full width at half-maximum.
+    //
+    // The search is probably not very efficient, but it's a first try.
+ 
+    Double_t p,x,fy,fxr,fxl;
+    Double_t step;
+    Double_t l,lold;
+    Int_t i = 0;
+    Int_t MAXCALLS = 10000;
+ 
+ 
+    // Search for maximum
+ 
+    p = params[1] - 0.1 * params[0];
+    step = 0.05 * params[0];
+    lold = -2.0;
+    l    = -1.0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = langaufun(&x,params);
+ 
+       if (l < lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-1);
+ 
+    maxx = x;
+ 
+    fy = l/2;
+ 
+ 
+    // Search for right x location of fy
+ 
+    p = maxx + params[0];
+    step = params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-2);
+ 
+    fxr = x;
+ 
+ 
+    // Search for left x location of fy
+ 
+    p = maxx - 0.5 * params[0];
+    step = -params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-3);
+ 
+ 
+    fxl = x;
+ 
+    FWHM = fxr - fxl;
+    return (0);
+ }
+ 
+void langaus() {
+    // Fill Histogram
+   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
+                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
+                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
+                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
+                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
+    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
+   */
+   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
+   TH1D * hSNR;
+
+    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
+    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+   
+    
+    // hSNR = h_beamSignal_b0[13];
+    // Fitting SNR histo
+    printf("Fitting...\n");
+    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
+    
+    // Setting fit range and start values
+    Double_t fr[2];
+    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
+    Double_t chisqr;
+    Int_t    ndf;
+    TF1 *fitsnr;
+    Double_t SNRPeak, SNRFWHM;
+    char rootfilename[50] = "";
+    Double_t norm;
+    ofstream myfile;
+    myfile.open ("MPVcorrection_carbon0mm05.txt");
+
+    
+    int j = 1;
+    for (int i = 0; i<26;i++){
+      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile *rootFile = new TFile(rootfilename,"OPEN");
+
+      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
+      norm = hSNR->GetEntries();
+      hSNR->Scale(1/norm);
+      
+      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
+      fr[1]=hSNR->GetMean() + 3*hSNR->GetRMS();
+ 
+      pllo[0]=0.0006; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01;
+      plhi[0]=0.01; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.05;
+      sv[0]=0.002;
+      //sv[1]=2.1;
+      sv[2]=0.01;
+      //   sv[3]=0.03;
+      sv[3] = hSNR->GetRMS()/2.;
+      sv[1] = hSNR->GetMean();
+    
+      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
+
+      
+      langaupro(fp,SNRPeak,SNRFWHM);
+      graph_x[i] = beta_carbon[i];
+      graph_y[i] = fp[1];
+      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
+      
+      printf("Fitting done\nPlotting results...\n");
+            myfile << i << " " << graph_y[i] << endl;
+
+      // Global style settings
+      gStyle->SetOptStat(1111);
+      gStyle->SetOptFit(111);
+      //gStyle->SetLabelSize(0.03,"x");
+      //gStyle->SetLabelSize(0.03,"y");
+      
+      // hSNR->GetXaxis()->SetRange(0,70);
+      hSNR->Draw();
+      fitsnr->Draw("lsame");
+      c1->Update();
+      //  c1->WaitPrimitive();
+      hSNR->Delete();
+      rootFile->Close();
+    }
+     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+   
+    graph_1->Draw("A*");
+    c2->SaveAs("carbon0mm5_SPratio.pdf");
+    myfile.close();
+
+    
+ }
diff --git a/carbon/plots.C b/carbon/plots.C
new file mode 100644
index 0000000..65b2e2c
--- /dev/null
+++ b/carbon/plots.C
@@ -0,0 +1,32 @@
+{
+
+  TGraph * gr_letratio = new TGraph();
+
+    char  finname[50];
+    int j = 0;
+    
+    Double_t x, q;
+    q = 0.5; // 0.5 for "median"
+    
+    
+    //  for (int j = 0; j<5;j++){
+    for (int i = 0; i<26;i++){
+
+      sprintf(finname, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile * f = new TFile(finname,"OPEN");
+      TH1F * h = (TH1F*)f->Get("h_spratio");
+      cout << i << " " << h->GetMean() << " +/- " << h->GetRMS() << endl;
+      h->ComputeIntegral(); // just a precaution
+      h->GetQuantiles(1, &x, &q);
+      std::cout << "median = " << x << std::endl;
+
+      
+      h->Delete();
+      f->Close();
+      
+    }
+    //  }
+  
+
+}
+
diff --git a/carbon/runJob.sh b/carbon/runJob.sh
new file mode 100755
index 0000000..d8d1dd2
--- /dev/null
+++ b/carbon/runJob.sh
@@ -0,0 +1,64 @@
+#!/bin/bash
+
+TYPE="HEAVYION"
+#energy=(0.04812 0.05982 0.07003 0.07917 0.08753 0.09530 0.10261 0.10956 0.11620 0.12257 0.12872 0.13465 0.14041 0.14601 0.15150 0.15689 0.16221 0.16746 0.17262 0.17770 0.18481 0.19187 0.19876 0.20548 0.21206 0.21851); #GeV/u, protons
+energy=(0.08883 0.11058 0.12979 0.14713 0.16309 0.17801 0.19213 0.20560 0.21852 0.23098 0.24303 0.25471 0.26608 0.27719 0.28810 0.29887 0.30952 0.32007 0.33048 0.34077 0.35522 0.36964 0.38378 0.39766 0.41132 0.42477); #GeV/u, carbon
+#energy=(0.05057 0.06186 0.07173 0.08064 0.08885 0.09652 0.10376 0.11064 0.11723 0.12355 0.12964 0.13555 0.14127 0.14684 0.15226 0.15756 0.16273 0.16780 0.17277 0.17764 0.18456 0.19154 0.19836 0.20503 0.21157 0.21798); #GeV/u, helium
+
+#THICKNESS=(0.025 0.050 0.075 0.100 1.25 0.150 0.2500 0.500) #cm
+THICKNESS=(0.025 0.050 0.100 0.150 0.500 1.000 0.125) #cm
+
+HOME=/work/leverington/fluka_gfortran7/project/carbon
+ZED=6.
+ATN=12.
+
+for j in 0 1 2 3 4 6
+do 
+    mkdir -p $HOME/jobs$j #make the directory if it doesn't exist
+    JOB_HOME=$HOME/jobs$j 
+    rm -r  $JOB_HOME/* #clean up the directory
+    cd $JOB_HOME
+    
+    
+    for i in {0..25} #26
+    do
+	# make a copy of the template
+	cp $HOME/runjobtemplate.inp  $JOB_HOME/runjob$i.inp
+	touch $JOB_HOME/runjob$i.inp
+	
+	#replace template placeholders with values
+	sed -i 's/ENERGYIN/'"-${energy[i]}"'/g' $JOB_HOME/runjob$i.inp #the negative sign indicates energy, not momentum
+	sed -i 's/TYPE/'"$TYPE"'/g' $JOB_HOME/runjob$i.inp  #beam particle type
+	sed -i 's/THICK/'"${THICKNESS[j]}"'/g' $JOB_HOME/runjob$i.inp  #thickness of polystyrene layer
+	sed -i 's/ZED/'"$ZED"'/g' $JOB_HOME/runjob$i.inp  #atomic charge
+	sed -i 's/ATN/'"$ATN"'/g' $JOB_HOME/runjob$i.inp  #atomic number
+
+	echo "Job: $j-$i   Ion type: $TYPE    Energy: ${energy[i]} GeV/u    Plane thickness: ${THICKNESS[j]} cm"
+	echo ""
+	
+	touch $JOB_HOME/fluka$j-$i.sh  #script to submit to the queue
+	echo "#!/bin/bash" >> $JOB_HOME/fluka$j-$i.sh
+	echo ". /local/env.sh" >> $JOB_HOME/fluka$j-$i.sh
+#	echo ". /cvmfs/lhcb.cern.ch/lib/lhcb/LBSCRIPTS/LBSCRIPTS_v9r2p2/InstallArea/scripts/LbLogin.sh -c x86_64-slc6-gcc7-opt" >> $JOB_HOME/fluka$j-$i.sh
+	echo "source /home/lhcb/leverington/setup.sh" >> $JOB_HOME/fluka$j-$i.sh
+        
+        #execute this command(s)
+	echo "cd "$JOB_HOME >> $JOB_HOME/fluka$j-$i.sh
+	echo "$FLUPRO/flutil/rfluka -e $FLUPRO/flutil/flukadpm3 -N0 -M1 runjob"$i >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "/work/leverington/readfluka/tools/eventdat2root runjob"$i"001_eventdata"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "$HOME/eventdatroot runjob"$i"001_eventdata "$j  >> $JOB_HOME/fluka$j-$i.sh
+	
+#	qsub -l os=slc6  -l ujl=20 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+	qsub -l os=slc6 -l hio=10 -l ujl=40 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+
+#	$FLUPRO/flutil/rfluka -e /home/leverington/software/fluka/flutil/flukadpm3 -N0 -M1 $JOB_HOME/runjob$i  #run 
+#	sleep 2
+#	eventdat2root runjob$i\001_eventdata  #convert to root file
+#	sleep 2
+#	$HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+ 	sleep 1
+    done
+    sleep 30
+done
diff --git a/carbon/runanalysis.sh b/carbon/runanalysis.sh
new file mode 100755
index 0000000..e48cf7e
--- /dev/null
+++ b/carbon/runanalysis.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+HOME=/work/leverington/fluka_gfortran7/project/carbon
+
+
+for j in {1..5}
+do
+    JOB_HOME=$HOME/jobs$j
+    cd $JOB_HOME
+
+	     for i in {0..25} #26
+	     do
+		 
+		 $HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+		 sleep 0.1 
+	     done
+done
diff --git a/carbon/runjobtemplate.inp b/carbon/runjobtemplate.inp
new file mode 100644
index 0000000..806aad9
--- /dev/null
+++ b/carbon/runjobtemplate.inp
@@ -0,0 +1,95 @@
+* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
+TITLE
+FLUKA Course Exercise
+*
+* use names everywhere and free format for geometry
+DEFAULTS                                                              HADROTHE
+*
+* beam definitions
+BEAM      ENERGYIN  0.0       -1.7      -0.8      -0.8                TYPE
+HI-PROPE     ZED       ATN
+BEAMPOS                             -0.1
+DELTARAY       1E-06       0.0       0.0  HYDROGEN  @LASTMAT          NOPRINT
+PART-THR  -0.00001  PROTON    PROTON                 0.0   
+*
+* Geometry
+* --------
+GEOBEGIN                                                              COMBNAME
+    0    0             Cylindrical Target
+*
+* Bodies
+* ------
+*
+* Blackhole to include geometry
+SPH BLK        0.0 0.0 0.0 10000.
+* Void sphere
+SPH VOID       0.0 0.0 0.0 1000.
+* Infinite cylinder
+ZCC TARG       0.0 0.0 5.
+*  planes cutting the cylinder
+XYP ZTlow      0.0
+XYP T1seg      THICK
+XYP T2seg      2.
+XYP ZThigh     10.
+*  additional plane for scoring
+XYP ZTscor     9.9999
+END
+*
+* Regions
+* -------
+*
+* Blackhole
+BLKHOLE      5 +BLK -VOID
+*
+* Target segment 1
+TARGS1       5 +TARG -ZTlow +T1seg
+* Target segment 2
+TARGS2       5 +TARG -T1seg +T2seg
+* Target segment 3
+TARGS3       5 +TARG -T2seg +ZTscor
+* Scoring region
+TARGS4       5 +TARG -ZTscor +ZThigh
+*
+* Air around target
+INAIR        5 | +VOID -TARG
+               | +VOID +ZTlow
+               | +VOID -ZThigh
+END
+*
+GEOEND
+*
+* Materials definition
+* --------------------
+MATERIAL                        0.001965                              CO2
+COMPOUND          1.    CARBON        2.    OXYGEN                    CO2
+*
+* Assign materials
+* ----------------
+ASSIGNMA    BLCKHOLE   BLKHOLE
+ASSIGNMA    POLYSTYR    TARGS1
+ASSIGNMA      VACUUM    TARGS2
+ASSIGNMA      VACUUM    TARGS3
+ASSIGNMA      VACUUM    TARGS4
+ASSIGNMA      VACUUM     INAIR
+*
+* Exercise: heavy ions
+* --------------------
+*
+*  charge spectrum of ions
+USRYIELD       1422.  HEAVYION      -68.    TARGS1    TARGS2        1.FragZ1
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS2    TARGS3        1.FragZ2
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS3    TARGS4        1.FragZ3
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+*
+*  LET in water of ions and charged particles
+USRYIELD       2223.  HEAVYION      -69.    TARGS3    TARGS4        1.LETHI
+USRYIELD         20.       0.0      200.       9.5       2.5     2703. &
+USRYIELD       2223.  ALL-CHAR      -69.    TARGS3    TARGS4        1.LETCh
+USRYIELD         20.       0.0      200.       9.5      -2.5     2703. &
+EVENTDAT        -21.                                                  eventdata
+*
+RANDOMIZ          1.
+START        100000.                                     0.0
+STOP
diff --git a/helium/MPVcorrection_helium0mm05.pdf b/helium/MPVcorrection_helium0mm05.pdf
new file mode 100644
index 0000000..142d046
Binary files /dev/null and b/helium/MPVcorrection_helium0mm05.pdf differ
diff --git a/helium/MPVcorrection_helium0mm05.txt b/helium/MPVcorrection_helium0mm05.txt
new file mode 100644
index 0000000..35ab77e
--- /dev/null
+++ b/helium/MPVcorrection_helium0mm05.txt
@@ -0,0 +1,26 @@
+0 1.03895
+1 1.02252
+2 1.01559
+3 1.01178
+4 1.00999
+5 1.00774
+6 1.00678
+7 1.00476
+8 1.00382
+9 1.0028
+10 1.00133
+11 1.00115
+12 0.999925
+13 0.999349
+14 0.997534
+15 0.996118
+16 0.99561
+17 0.995045
+18 0.992608
+19 0.991975
+20 0.991158
+21 0.990081
+22 0.988478
+23 0.986514
+24 0.985795
+25 0.984543
diff --git a/helium/Makefile b/helium/Makefile
new file mode 100644
index 0000000..c16461b
--- /dev/null
+++ b/helium/Makefile
@@ -0,0 +1,38 @@
+
+ROOTCFLAGS   := $(shell root-config --cflags)
+ROOTLIBS     := $(shell root-config --libs)
+ROOTGLIBS    := $(shell root-config --glibs)
+
+GSLCFLAGS   := $(shell gsl-config --cflags)
+GSLLIBS     := $(shell gsl-config --libs)
+GSLGLIBS    := $(shell gsl-config --glibs)
+
+
+LDFLAGS       = -O
+LIBS         += $(ROOTLIBS) $(GSLLIBS) -L/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/lib
+CFLAGS       += $(ROOTCFLAGS) $(GSLCFLAGS)  -I/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/include
+
+#  Not sure why Minuit isn't being included -- put in by hand
+#
+LIBS         += -lMinuit -ltbb
+
+
+all: eventdatroot
+
+
+
+eventdatroot: eventdatroot.o
+	g++ -o eventdatroot  eventdatroot.o $(LDFLAGS) $(LIBS)
+
+%.o: %.c
+	g++ ${CFLAGS} -c  -g $<
+
+
+test:
+	@echo $(ROOTCFLAGS) 
+	@echo $(LDFLAGS) 
+	@echo $(LIBS)
+
+clean:	
+	-rm -f  *.o
+
diff --git a/helium/eventdatroot b/helium/eventdatroot
new file mode 100755
index 0000000..3fc4a09
Binary files /dev/null and b/helium/eventdatroot differ
diff --git a/helium/eventdatroot.c b/helium/eventdatroot.c
new file mode 100644
index 0000000..ceeb87b
--- /dev/null
+++ b/helium/eventdatroot.c
@@ -0,0 +1,143 @@
+#define eventdatroot_cxx
+#include "eventdatroot.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+//#include <TROOT.h>
+#include <TF1.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+
+
+using namespace std;
+
+void eventdatroot::Loop(Float_t scale)
+{
+  //   In a ROOT session, you can do:
+  //      root> .L eventdatroot.C
+  //      root> eventdatroot t
+  //      root> t.GetEntry(12); // Fill t data members with entry number 12
+  //      root> t.Show();       // Show values of entry 12
+  //      root> t.Show(16);     // Read and show values of entry 16
+  //      root> t.Loop();       // Loop on all entries
+  //
+
+  //     This is the loop skeleton where:
+  //    jentry is the global entry number in the chain
+  //    ientry is the entry number in the current Tree
+  //  Note that the argument to GetEntry must be:
+  //    jentry for TChain::GetEntry
+  //    ientry for TTree::GetEntry and TBranch::GetEntry
+  //
+  //       To read only selected branches, Insert statements like:
+  // METHOD1:
+  //    fChain->SetBranchStatus("*",0);  // disable all branches
+  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
+  // METHOD2: replace line
+  //    fChain->GetEntry(jentry);       //read all branches
+  //by  b_branchname->GetEntry(ientry); //read only this branch
+
+  Float_t ratio_electron, ratio_primary, etotal, spratio;
+  // TF1 * f_sp_ps = new TF1("f_sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in polystyrene [MeV cm2 /g]
+  // f_sp_ps->SetParameters(361.936, -0.892255, 1.19133);  //protons between 50 and 220 MeV
+  //  TF1 * f_c12sp_ps = new TF1("f_c12sp_ps","[0]*pow(x,[1])+[2]", 80, 480); //stopping power of carbon in polystyrene [MeV cm2 /g]
+  //  f_c12sp_ps->SetParameters(14538.9, -0.922423,49.2859);  //carbon between 80 and 480 MeV
+
+  TF1 * f_h2sp_ps = new TF1("f_h2sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of helium in polystyrene [MeV cm2 /g]
+  f_h2sp_ps->SetParameters(1387.08, -0.882686,4.60833);  //heelium between 50 and 300 MeV
+
+
+  
+  if (fChain == 0) return;
+
+  Long64_t nentries = fChain->GetEntriesFast();
+
+  Long64_t nbytes = 0, nb = 0;
+  for (Long64_t jentry=0; jentry<nentries;jentry++) {
+    Long64_t ientry = LoadTree(jentry);
+    if (ientry < 0) break;
+    nb = fChain->GetEntry(jentry);   nbytes += nb;
+    // if (Cut(ientry) < 0) continue;
+    //    cout << jentry << endl;
+
+    h_endist_0->Fill(ENDIST[0]); // edep by ionisation
+    h_endist_1->Fill(ENDIST[1]); // edep by pi0, e-, e+ and #gamma
+    h_endist_2->Fill(ENDIST[2]); // edep by nuclear recoils and heavy fragments
+    h_endist_3->Fill(ENDIST[3]); // edep by part. #lt threshold
+    h_endist_4->Fill(ENDIST[4]); // E leaving the syste
+    h_endist_5->Fill(ENDIST[5]); // E carried by discarded particles
+    h_endist_6->Fill(ENDIST[6]); // resid, excit. E after evap.
+    h_endist_7->Fill(ENDIST[7]); // edep by low-energy neutrons
+    h_endist_8->Fill(ENDIST[8]); // E of part. #gt time limit
+    h_endist_9->Fill(ENDIST[9]); // E lost in endothermic nuclear reactions
+    h_endist_10->Fill(ENDIST[10]); // E lost in endothermic low-E n-reaction
+    h_endist_11->Fill(ENDIST[11]); // missing E
+
+    etotal = 0.;
+    for (unsigned int i = 0; i<3; i++ ){
+      etotal += ENDIST[i];
+    }
+
+    ratio_electron = ENDIST[1] / etotal;
+    ratio_primary = ENDIST[0] / etotal;
+
+    h_ratio_e->Fill(ratio_electron);
+    h_ratio_p->Fill(ratio_primary);
+    h_sum_ep->Fill(ENDIST[1]+ENDIST[0]);
+
+    h_let_ep->Fill( (ENDIST[1]+ENDIST[0])/scale/4. ); //GeV/cm/Z^2
+    h_let_e->Fill( (ENDIST[1])/scale/4. );
+    h_let_p->Fill( (ENDIST[0])/scale/4. );
+
+    spratio = ( ( (ENDIST[0]+ENDIST[1])/scale*1000./1.06) /  f_h2sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000./4. ) ); //ratio of deposited energy to MSTAR stopping power (MeV/cm/density) / MeVcm2/g
+    h_spratio->Fill(spratio);
+
+    //   cout <<  ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06)  << " " << f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) << " " << spratio << endl;
+
+  }
+}
+
+
+
+int main(int argc, const char* argv[])
+{
+  Int_t thickkey = atoi(argv[2]);
+  // Float_t THICKNESS[6] = {0.025, 0.050, 0.075, 0.100, 0.0125, 0.150, 0.2500, 0.500}; //#cm
+  Float_t THICKNESS[6] = {0.025, 0.050, 0.100, 0.150,  0.500, 0.125 }; //#cm
+
+  // Float_t THICKNESS[5] = {0.050, 0.100, 0.150, 1.000, 10.00}; //#cm
+
+  ///open the input and output root files
+  char  finname[50];
+  sprintf(finname, "%s.root",argv[1]);
+  TFile * f = new TFile(finname,"OPEN");
+
+  char  foutname[50];
+  sprintf(foutname, "%s_out.root",argv[1]);
+  cout << foutname << endl;
+  TFile * fout = new TFile(foutname, "RECREATE");
+
+  //run the analysis loop above
+  eventdatroot t;
+  TTree * tree;
+  f->GetObject("EVENTDAT",tree);  
+  t.Init(tree);
+  t.Loop( THICKNESS[thickkey] );
+  t.Closefile(fout);
+
+  return 1;
+
+}
diff --git a/helium/eventdatroot.h b/helium/eventdatroot.h
new file mode 100644
index 0000000..f4cb35a
--- /dev/null
+++ b/helium/eventdatroot.h
@@ -0,0 +1,199 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Mar  8 13:05:30 2019 by ROOT version 6.12/06
+// from TTree EVENTDAT/FLUKA Course Exercise	DATE:  3/ 4/19,  TIME: 11:54:46
+// found on file: jobs/runjob1001_eventdata.root
+//////////////////////////////////////////////////////////
+
+#ifndef eventdatroot_h
+#define eventdatroot_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1.h>
+
+// Header file for the classes stored in the TTree if any.
+
+class eventdatroot {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Float_t         DATA_ALL_PART[6];
+   Float_t         DATA_BEAMPART[6];
+   Float_t         DATA_ENERGY[6];
+   Float_t         DATA_id211[6];
+   Int_t           SEEDS[4];
+   Float_t         ENDIST[12];
+
+   // List of branches
+   TBranch        *b_DATA;   //!
+   TBranch        *b_seed;   //!
+   TBranch        *b_endist;   //!
+
+      
+   TH1F * h_endist_0;
+   TH1F * h_endist_1;
+   TH1F * h_endist_2;
+   TH1F * h_endist_3;
+   TH1F * h_endist_4;
+   TH1F * h_endist_5;
+   TH1F * h_endist_6;
+   TH1F * h_endist_7;
+   TH1F * h_endist_8;
+   TH1F * h_endist_9;
+   TH1F * h_endist_10;
+   TH1F * h_endist_11;
+
+   TH1F * h_ratio_e;
+   TH1F * h_ratio_p;
+   TH1F * h_sum_ep;
+   
+   TH1F * h_let_ep;
+   TH1F * h_let_e;
+   TH1F * h_let_p;
+
+   TH1F * h_spratio;
+
+   eventdatroot(TTree *tree=0);
+   virtual ~eventdatroot();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop(Float_t scale);
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+   virtual Int_t    Closefile(TFile * file);
+};
+
+#endif
+
+#ifdef eventdatroot_cxx
+eventdatroot::eventdatroot(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+  //  if (tree == 0) {
+     //   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("jobs/runjob1001_eventdata.root");
+      //     if (!f || !f->IsOpen()) {
+		//      f = new TFile("jobs/runjob1001_eventdata.root");
+	 //     }
+      //    f->GetObject("EVENTDAT",tree);
+
+      //   }
+   //   Init(tree);
+}
+
+eventdatroot::~eventdatroot()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t eventdatroot::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t eventdatroot::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void eventdatroot::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("DATA", DATA_ALL_PART, &b_DATA);
+   fChain->SetBranchAddress("SEEDS", SEEDS, &b_seed);
+   fChain->SetBranchAddress("ENDIST", ENDIST, &b_endist);
+
+   h_endist_0 = new TH1F("h_endist_0","edep by ionisation",400,0.0,0.05);
+   h_endist_1 = new TH1F("h_endist_1","edep by pi0, e-, e+ and #gamma",400,0.0,0.03);
+   h_endist_2 = new TH1F("h_endist_2","edep by nuclear recoils and heavy fragments",400,0.0,0.001);
+   h_endist_3 = new TH1F("h_endist_3","edep by part. #lt threshold",200,0.0,0.001);
+   h_endist_4 = new TH1F("h_endist_4","E leaving the system",200,0.0,1.1);
+   h_endist_5 = new TH1F("h_endist_5","E carried by discarded particles",200,0.0,0.1);
+   h_endist_6 = new TH1F("h_endist_6","resid, excit. E after evap.",200,0.0,0.1);
+   h_endist_7 = new TH1F("h_endist_7","edep by low-energy neutrons",200,0.0,0.1);
+   h_endist_8 = new TH1F("h_endist_8","E of part. #gt time limit",200,0.0,0.1);
+   h_endist_9 = new TH1F("h_endist_9","E lost in endothermic nuclear reactions",200,0.0,0.1);
+   h_endist_10 = new TH1F("h_endist_10","E lost in endothermic low-E n-reactions",200,0.0,0.1);
+   h_endist_11 = new TH1F("h_endist_11","missing E",400,0.0,0.01);
+
+   h_ratio_e = new TH1F("h_ratio_e","fraction of edep by pi0, e-, e+ and #gamma  ",400,0.0,1.0);
+   h_ratio_p = new TH1F("h_ratio_p","fraction of edep by ionisation",400,0.0,1.0);
+   h_sum_ep = new TH1F("h_sum_ep","sum of edep by ion, pi0, e-, e+ and #gamma,",400,0.0,0.1);
+
+   h_let_ep = new TH1F("h_let_ep","LET/Z^{2} (GeV/cm) by ion, pi0, e-, e+ and #gamma,",200,0.0,0.05);
+   h_let_e = new TH1F("h_let_e","LET/Z^{2} (GeV/cm) by pi0, e-, e+ and #gamma,",200,0.0,0.02);
+   h_let_p = new TH1F("h_let_p","LET/Z^{2} (GeV/cm) by ion",200,0.0,0.02);
+
+   h_spratio = new TH1F("h_spratio","total LET/SP (MeV/cm/#rho)/SP(MeVcm2/g)",200,0.5,3.0);
+
+   
+   Notify();
+}
+
+Bool_t eventdatroot::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void eventdatroot::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t eventdatroot::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+
+Int_t eventdatroot::Closefile(TFile * file)
+{
+  file->cd();
+  file->Write();
+  file->Close();
+
+  return 1;
+}
+
+#endif // #ifdef eventdatroot_cxx
diff --git a/helium/langaus.C b/helium/langaus.C
new file mode 100644
index 0000000..2bcbd45
--- /dev/null
+++ b/helium/langaus.C
@@ -0,0 +1,355 @@
+ /// \ingroup tutorial_fit
+ /// \notebook
+ /// Convoluted Landau and Gaussian Fitting Function
+ ///         (using ROOT's Landau and Gauss functions)
+ ///
+ ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
+ ///
+ ///  to execute this example, do:
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C
+ /// ~~~
+ ///
+ /// or
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C++
+ /// ~~~
+ ///
+ /// \macro_image
+ /// \macro_output
+ /// \macro_code
+ ///
+ /// \authors H.Pernegger, Markus Friedl
+ 
+ #include "TH1.h"
+ #include "TF1.h"
+ #include "TROOT.h"
+ #include "TStyle.h"
+ #include "TMath.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+ 
+ Double_t langaufun(Double_t *x, Double_t *par) {
+ 
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+ 
+       // Numeric constants
+       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+       Double_t mpshift  = -0.22278298;       // Landau maximum location
+ 
+       // Control constants
+       Double_t np = 200.0;      // number of convolution steps
+       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
+ 
+       // Variables
+       Double_t xx;
+       Double_t mpc;
+       Double_t fland;
+       Double_t sum = 0.0;
+       Double_t xlow,xupp;
+       Double_t step;
+       Double_t i;
+ 
+ 
+       // MP shift correction
+       mpc = par[1] - mpshift * par[0];
+ 
+       // Range of convolution integral
+       xlow = x[0] - sc * par[3];
+       xupp = x[0] + 2*sc * par[3];
+ 
+       step = (xupp-xlow) / np;
+ 
+       // Convolution integral of Landau and Gaussian by sum
+       for(i=1.0; i<=np/2; i++) {
+          xx = xlow + (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+ 
+          xx = xupp - (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+       }
+ 
+       return (par[2] * step * sum * invsq2pi / par[3]);
+ }
+ 
+ 
+ 
+ TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
+ {
+    // Once again, here are the Landau * Gaussian parameters:
+    //   par[0]=Width (scale) parameter of Landau density
+    //   par[1]=Most Probable (MP, location) parameter of Landau density
+    //   par[2]=Total area (integral -inf to inf, normalization constant)
+    //   par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    // Variables for langaufit call:
+    //   his             histogram to fit
+    //   fitrange[2]     lo and hi boundaries of fit range
+    //   startvalues[4]  reasonable start values for the fit
+    //   parlimitslo[4]  lower parameter limits
+    //   parlimitshi[4]  upper parameter limits
+    //   fitparams[4]    returns the final fit parameters
+    //   fiterrors[4]    returns the final fit errors
+    //   ChiSqr          returns the chi square
+    //   NDF             returns ndf
+ 
+    Int_t i;
+    Char_t FunName[100];
+ 
+    sprintf(FunName,"Fitfcn_%s",his->GetName());
+ 
+    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
+    if (ffitold) delete ffitold;
+ 
+    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
+    ffit->SetParameters(startvalues);
+    ffit->SetParNames("Width","MP","Area","GSigma");
+ 
+    for (i=0; i<4; i++) {
+       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
+    }
+ 
+    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
+ 
+    ffit->GetParameters(fitparams);    // obtain fit parameters
+    for (i=0; i<4; i++) {
+       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
+    }
+    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
+    NDF[0] = ffit->GetNDF();           // obtain ndf
+ 
+    return (ffit);              // return fit function
+ 
+ }
+ 
+ 
+ Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
+ 
+    // Seaches for the location (x value) at the maximum of the
+    // Landau-Gaussian convolute and its full width at half-maximum.
+    //
+    // The search is probably not very efficient, but it's a first try.
+ 
+    Double_t p,x,fy,fxr,fxl;
+    Double_t step;
+    Double_t l,lold;
+    Int_t i = 0;
+    Int_t MAXCALLS = 10000;
+ 
+ 
+    // Search for maximum
+ 
+    p = params[1] - 0.1 * params[0];
+    step = 0.05 * params[0];
+    lold = -2.0;
+    l    = -1.0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = langaufun(&x,params);
+ 
+       if (l < lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-1);
+ 
+    maxx = x;
+ 
+    fy = l/2;
+ 
+ 
+    // Search for right x location of fy
+ 
+    p = maxx + params[0];
+    step = params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-2);
+ 
+    fxr = x;
+ 
+ 
+    // Search for left x location of fy
+ 
+    p = maxx - 0.5 * params[0];
+    step = -params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-3);
+ 
+ 
+    fxl = x;
+ 
+    FWHM = fxr - fxl;
+    return (0);
+ }
+ 
+void langaus() {
+    // Fill Histogram
+   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
+                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
+                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
+                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
+                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
+    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
+   */
+   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
+   TH1D * hSNR;
+
+    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
+    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+   
+    
+    // hSNR = h_beamSignal_b0[13];
+    // Fitting SNR histo
+    printf("Fitting...\n");
+    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
+    
+    // Setting fit range and start values
+    Double_t fr[2];
+    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
+    Double_t chisqr;
+    Int_t    ndf;
+    TF1 *fitsnr;
+    Double_t SNRPeak, SNRFWHM;
+    char rootfilename[50] = "";
+    char saveplotname[50] = "";
+    char fout_mpv_name[50] = "";
+    Double_t norm;
+    
+    ofstream myfile, fout_mpv;
+    int j = 1;
+
+    myfile.open ("MPVcorrection_helium0mm05.txt");
+    sprintf(fout_mpv_name, "jobs%i/plots/mpv.txt",j);
+     fout_mpv.open (fout_mpv_name);
+    
+    for (int i = 0; i<26;i++){
+      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile *rootFile = new TFile(rootfilename,"OPEN");
+
+      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
+      norm = hSNR->GetEntries();
+      hSNR->Scale(1/norm);
+      
+      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
+      fr[1]=hSNR->GetMean() + 3*hSNR->GetRMS();
+ 
+      pllo[0]=0.0006; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01;
+      plhi[0]=0.01; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.05;
+      sv[0]= hSNR->GetRMS()/10.;;
+      //sv[1]=2.1;
+      sv[2]=0.01;
+      //   sv[3]=0.03;
+      sv[3] = hSNR->GetRMS()/2.;
+      sv[1] = hSNR->GetMean();
+    
+      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
+
+      
+      langaupro(fp,SNRPeak,SNRFWHM);
+      graph_x[i] = beta_helium[i];
+      graph_y[i] = fp[1];
+      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
+      fout_mpv << graph_x[i] << " " << graph_y[i] << " " << graph_yerr[i] << endl;
+      printf("Fitting done\nPlotting results...\n");
+            myfile << i << " " << graph_y[i] << endl;
+
+      // Global style settings
+      gStyle->SetOptStat(1111);
+      gStyle->SetOptFit(111);
+      //gStyle->SetLabelSize(0.03,"x");
+      //gStyle->SetLabelSize(0.03,"y");
+      
+      // hSNR->GetXaxis()->SetRange(0,70);
+      hSNR->Draw();
+      fitsnr->Draw("lsame");
+      c1->Update();
+      sprintf(saveplotname, "jobs%i/plots/runjob%i001_h_spratio.pdf",j,i);
+      c1->SaveAs(saveplotname);
+      //  c1->WaitPrimitive();
+      hSNR->Delete();
+      rootFile->Close();
+    }
+     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+    graph_1->Draw("A*");
+    c2->SaveAs("MPVcorrection_helium0mm05.pdf");
+    myfile.close();
+    fout_mpv.close();
+
+    
+ }
diff --git a/helium/plots.C b/helium/plots.C
new file mode 100644
index 0000000..8b0643d
--- /dev/null
+++ b/helium/plots.C
@@ -0,0 +1,63 @@
+{
+
+
+   Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+    
+  TGraph * gr_letratio = new TGraph();
+
+    
+    char  finname[50] = "";
+    char fout_mean_name[50] = "";
+    char fout_median_name[50] = "";
+    char fout_mpv_name[50] = "";
+
+    ofstream fout_mean, fout_median;
+    int j = 1;
+        sprintf(fout_mean_name, "jobs%i/plots/mean.txt",j);
+      sprintf(fout_median_name, "jobs%i/plots/median.txt",j);
+      sprintf(fout_mpv_name, "jobs%i/plots/mpv.txt",j);
+
+    
+    fout_median.open (fout_median_name);
+    fout_mean.open (fout_mean_name);
+  
+    Double_t x, q;
+    q = 0.5; // 0.5 for "median"
+    
+    
+    //  for (int j = 0; j<5;j++){
+    for (int i = 0; i<26;i++){
+
+      sprintf(finname, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile * f = new TFile(finname,"OPEN");
+      TH1F * h = (TH1F*)f->Get("h_spratio");
+      cout << i << " " << h->GetMean() << " +/- " << h->GetRMS() << endl;
+      fout_mean  << beta_helium[i] << " " << h->GetMean() << " " << 0.012*h->GetMean() << endl;
+      h->ComputeIntegral(); // just a precaution
+      h->GetQuantiles(1, &x, &q);
+      std::cout << "median = " << x << std::endl;
+      fout_median  << beta_helium[i] << " " << x << " " << 0.012*x << endl;
+
+      
+      h->Delete();
+      f->Close();
+      
+    }
+    //  }
+      TGraphErrors * graph_1 = new TGraphErrors(fout_mean_name,"%lg %lg %lg" );
+    TGraphErrors * graph_2 = new TGraphErrors(fout_median_name,"%lg %lg %lg" );
+    TGraphErrors * graph_3 = new TGraphErrors(fout_mpv_name,"%lg %lg %lg" );
+
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+    graph_3->Draw("A*");
+    graph_1->Draw("*");
+    graph_2->Draw("*");
+    fout_mean.close();
+    fout_median.close();
+
+}
+
diff --git a/helium/runJob.sh b/helium/runJob.sh
new file mode 100755
index 0000000..e4e0c29
--- /dev/null
+++ b/helium/runJob.sh
@@ -0,0 +1,63 @@
+#!/bin/bash
+
+TYPE="4-HELIUM"
+#energy=(0.04812 0.05982 0.07003 0.07917 0.08753 0.09530 0.10261 0.10956 0.11620 0.12257 0.12872 0.13465 0.14041 0.14601 0.15150 0.15689 0.16221 0.16746 0.17262 0.17770 0.18481 0.19187 0.19876 0.20548 0.21206 0.21851); #GeV/u, protons
+#energy=(0.08883 0.11058 0.12979 0.14713 0.16309 0.17801 0.19213 0.20560 0.21852 0.23098 0.24303 0.25471 0.26608 0.27719 0.28810 0.29887 0.30952 0.32007 0.33048 0.34077 0.35522 0.36964 0.38378 0.39766 0.41132 0.42477); #GeV/u, carbon
+energy=(0.05057 0.06186 0.07173 0.08064 0.08885 0.09652 0.10376 0.11064 0.11723 0.12355 0.12964 0.13555 0.14127 0.14684 0.15226 0.15756 0.16273 0.16780 0.17277 0.17764 0.18456 0.19154 0.19836 0.20503 0.21157 0.21798); #GeV/u, helium
+
+#THICKNESS=(0.025 0.050 0.075 0.100 1.25 0.150 0.2500 0.500) #cm
+THICKNESS=(0.025 0.050 0.100 0.150 0.500 1.000 0.125) #cm
+
+HOME=/work/leverington/fluka_gfortran7/project/helium
+ZED=2.
+ATN=4.
+
+for j in {0..6}
+do 
+    mkdir -p $HOME/jobs$j #make the directory if it doesn't exist
+    JOB_HOME=$HOME/jobs$j 
+    rm -r  $JOB_HOME/* #clean up the directory
+    cd $JOB_HOME
+    
+    
+    for i in {0..25} #26
+    do
+	# make a copy of the template
+	cp $HOME/runjobtemplate.inp  $JOB_HOME/runjob$i.inp
+	touch $JOB_HOME/runjob$i.inp
+	
+	#replace template placeholders with values
+	sed -i 's/ENERGYIN/'"-${energy[i]}"'/g' $JOB_HOME/runjob$i.inp #the negative sign indicates energy, not momentum
+	sed -i 's/TYPE/'"$TYPE"'/g' $JOB_HOME/runjob$i.inp  #beam particle type
+	sed -i 's/THICK/'"${THICKNESS[j]}"'/g' $JOB_HOME/runjob$i.inp  #thickness of polystyrene layer
+	sed -i 's/ZED/'"$ZED"'/g' $JOB_HOME/runjob$i.inp  #atomic charge
+	sed -i 's/ATN/'"$ATN"'/g' $JOB_HOME/runjob$i.inp  #atomic number
+
+	echo "Job: $j-$i   Ion type: $TYPE    Energy: ${energy[i]} GeV/u    Plane thickness: ${THICKNESS[j]} cm"
+	echo ""
+	
+	touch $JOB_HOME/fluka$j-$i.sh  #script to submit to the queue
+	echo "#!/bin/bash" >> $JOB_HOME/fluka$j-$i.sh
+	echo ". /local/env.sh" >> $JOB_HOME/fluka$j-$i.sh
+#	echo ". /cvmfs/lhcb.cern.ch/lib/lhcb/LBSCRIPTS/LBSCRIPTS_v9r2p2/InstallArea/scripts/LbLogin.sh -c x86_64-slc6-gcc7-opt" >> $JOB_HOME/fluka$j-$i.sh
+	echo "source /home/lhcb/leverington/setup.sh" >> $JOB_HOME/fluka$j-$i.sh
+        
+        #execute this command(s)
+	echo "cd "$JOB_HOME >> $JOB_HOME/fluka$j-$i.sh
+	echo "$FLUPRO/flutil/rfluka -e $FLUPRO/flutil/flukadpm3 -N0 -M1 runjob"$i >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "/work/leverington/readfluka/tools/eventdat2root runjob"$i"001_eventdata"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "$HOME/eventdatroot runjob"$i"001_eventdata "$j  >> $JOB_HOME/fluka$j-$i.sh
+	
+#	qsub -l os=slc6  -l ujl=20 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+	qsub -l os=slc6 -l hio=10 -l ujl=40 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+
+#	$FLUPRO/flutil/rfluka -e /home/leverington/software/fluka/flutil/flukadpm3 -N0 -M1 $JOB_HOME/runjob$i  #run 
+#	sleep 2
+#	eventdat2root runjob$i\001_eventdata  #convert to root file
+#	sleep 2
+#	$HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+ 	sleep 1
+    done
+done
diff --git a/helium/runanalysis.sh b/helium/runanalysis.sh
new file mode 100755
index 0000000..81b60ac
--- /dev/null
+++ b/helium/runanalysis.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+HOME=/work/leverington/fluka_gfortran7/project/helium
+
+
+for j in {0..5}
+do
+    JOB_HOME=$HOME/jobs$j
+    cd $JOB_HOME
+
+	     for i in {0..25} #26
+	     do
+		 
+		 $HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+		 sleep 0.1 
+	     done
+done
diff --git a/helium/runjobtemplate.inp b/helium/runjobtemplate.inp
new file mode 100644
index 0000000..868cc24
--- /dev/null
+++ b/helium/runjobtemplate.inp
@@ -0,0 +1,95 @@
+* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
+TITLE
+FLUKA Course Exercise
+*
+* use names everywhere and free format for geometry
+DEFAULTS                                                              HADROTHE
+*
+* beam definitions
+BEAM      ENERGYIN  0.0       -1.7      -0.8      -0.8                TYPE
+*HI-PROPE     ZED       ATN
+BEAMPOS                             -0.1
+DELTARAY       1E-06       0.0       0.0  HYDROGEN  @LASTMAT          NOPRINT
+*PART-THR  -0.00001  PROTON    PROTON                 0.0   
+*
+* Geometry
+* --------
+GEOBEGIN                                                              COMBNAME
+    0    0             Cylindrical Target
+*
+* Bodies
+* ------
+*
+* Blackhole to include geometry
+SPH BLK        0.0 0.0 0.0 10000.
+* Void sphere
+SPH VOID       0.0 0.0 0.0 1000.
+* Infinite cylinder
+ZCC TARG       0.0 0.0 5.
+*  planes cutting the cylinder
+XYP ZTlow      0.0
+XYP T1seg      THICK
+XYP T2seg      2.
+XYP ZThigh     10.
+*  additional plane for scoring
+XYP ZTscor     9.9999
+END
+*
+* Regions
+* -------
+*
+* Blackhole
+BLKHOLE      5 +BLK -VOID
+*
+* Target segment 1
+TARGS1       5 +TARG -ZTlow +T1seg
+* Target segment 2
+TARGS2       5 +TARG -T1seg +T2seg
+* Target segment 3
+TARGS3       5 +TARG -T2seg +ZTscor
+* Scoring region
+TARGS4       5 +TARG -ZTscor +ZThigh
+*
+* Air around target
+INAIR        5 | +VOID -TARG
+               | +VOID +ZTlow
+               | +VOID -ZThigh
+END
+*
+GEOEND
+*
+* Materials definition
+* --------------------
+MATERIAL                        0.001965                              CO2
+COMPOUND          1.    CARBON        2.    OXYGEN                    CO2
+*
+* Assign materials
+* ----------------
+ASSIGNMA    BLCKHOLE   BLKHOLE
+ASSIGNMA    POLYSTYR    TARGS1
+ASSIGNMA      VACUUM    TARGS2
+ASSIGNMA      VACUUM    TARGS3
+ASSIGNMA      VACUUM    TARGS4
+ASSIGNMA      VACUUM     INAIR
+*
+* Exercise: heavy ions
+* --------------------
+*
+*  charge spectrum of ions
+USRYIELD       1422.  HEAVYION      -68.    TARGS1    TARGS2        1.FragZ1
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS2    TARGS3        1.FragZ2
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS3    TARGS4        1.FragZ3
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+*
+*  LET in water of ions and charged particles
+USRYIELD       2223.  HEAVYION      -69.    TARGS3    TARGS4        1.LETHI
+USRYIELD         20.       0.0      200.       9.5       2.5     2703. &
+USRYIELD       2223.  ALL-CHAR      -69.    TARGS3    TARGS4        1.LETCh
+USRYIELD         20.       0.0      200.       9.5      -2.5     2703. &
+EVENTDAT        -21.                                                  eventdata
+*
+RANDOMIZ          1.
+START        100000.                                     0.0
+STOP
diff --git a/oxygen/MPVcorrection_oxygen0mm05.txt b/oxygen/MPVcorrection_oxygen0mm05.txt
new file mode 100644
index 0000000..a189427
--- /dev/null
+++ b/oxygen/MPVcorrection_oxygen0mm05.txt
@@ -0,0 +1,25 @@
+0 0.986872
+1 0.97926
+2 0.975462
+3 0.969715
+4 0.964747
+5 0.961425
+6 0.957963
+7 0.955564
+8 0.954657
+9 0.951364
+11 0.948485
+12 0.946011
+13 0.944629
+14 0.944191
+15 0.941793
+16 0.941145
+17 0.940926
+18 0.940841
+19 0.938605
+20 0.94052
+21 0.938271
+22 0.936244
+23 0.939211
+24 0.936627
+25 0.935963
diff --git a/oxygen/Makefile b/oxygen/Makefile
new file mode 100644
index 0000000..c16461b
--- /dev/null
+++ b/oxygen/Makefile
@@ -0,0 +1,38 @@
+
+ROOTCFLAGS   := $(shell root-config --cflags)
+ROOTLIBS     := $(shell root-config --libs)
+ROOTGLIBS    := $(shell root-config --glibs)
+
+GSLCFLAGS   := $(shell gsl-config --cflags)
+GSLLIBS     := $(shell gsl-config --libs)
+GSLGLIBS    := $(shell gsl-config --glibs)
+
+
+LDFLAGS       = -O
+LIBS         += $(ROOTLIBS) $(GSLLIBS) -L/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/lib
+CFLAGS       += $(ROOTCFLAGS) $(GSLCFLAGS)  -I/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/include
+
+#  Not sure why Minuit isn't being included -- put in by hand
+#
+LIBS         += -lMinuit -ltbb
+
+
+all: eventdatroot
+
+
+
+eventdatroot: eventdatroot.o
+	g++ -o eventdatroot  eventdatroot.o $(LDFLAGS) $(LIBS)
+
+%.o: %.c
+	g++ ${CFLAGS} -c  -g $<
+
+
+test:
+	@echo $(ROOTCFLAGS) 
+	@echo $(LDFLAGS) 
+	@echo $(LIBS)
+
+clean:	
+	-rm -f  *.o
+
diff --git a/oxygen/eventdatroot b/oxygen/eventdatroot
new file mode 100755
index 0000000..a62921f
Binary files /dev/null and b/oxygen/eventdatroot differ
diff --git a/oxygen/eventdatroot.c b/oxygen/eventdatroot.c
new file mode 100644
index 0000000..9470138
--- /dev/null
+++ b/oxygen/eventdatroot.c
@@ -0,0 +1,139 @@
+#define eventdatroot_cxx
+#include "eventdatroot.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+//#include <TROOT.h>
+#include <TF1.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+
+
+using namespace std;
+
+void eventdatroot::Loop(Float_t scale)
+{
+  //   In a ROOT session, you can do:
+  //      root> .L eventdatroot.C
+  //      root> eventdatroot t
+  //      root> t.GetEntry(12); // Fill t data members with entry number 12
+  //      root> t.Show();       // Show values of entry 12
+  //      root> t.Show(16);     // Read and show values of entry 16
+  //      root> t.Loop();       // Loop on all entries
+  //
+
+  //     This is the loop skeleton where:
+  //    jentry is the global entry number in the chain
+  //    ientry is the entry number in the current Tree
+  //  Note that the argument to GetEntry must be:
+  //    jentry for TChain::GetEntry
+  //    ientry for TTree::GetEntry and TBranch::GetEntry
+  //
+  //       To read only selected branches, Insert statements like:
+  // METHOD1:
+  //    fChain->SetBranchStatus("*",0);  // disable all branches
+  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
+  // METHOD2: replace line
+  //    fChain->GetEntry(jentry);       //read all branches
+  //by  b_branchname->GetEntry(ientry); //read only this branch
+
+  Float_t ratio_electron, ratio_primary, etotal, spratio;
+  // TF1 * f_sp_ps = new TF1("f_sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in polystyrene [MeV cm2 /g]
+  // f_sp_ps->SetParameters(361.936, -0.892255, 1.19133);  //protons between 50 and 220 MeV
+  // TF1 * f_c12sp_ps = new TF1("f_c12sp_ps","[0]*pow(x,[1])+[2]", 80, 480); //stopping power of carbon in polystyrene [MeV cm2 /g]
+  // f_c12sp_ps->SetParameters(14538.9, -0.922423,49.2859);  //carbon between 80 and 480 MeV
+TF1 * f_o16sp_ps = new TF1("f_o16sp_ps","[0]*pow(x,[1])+[2]", 80,480); //stopping power of oxygen in polystyrene [MeV cm2 /g]
+f_o16sp_ps->SetParameters(26465.6, -0.928309,88.6728);  //oxygen between 80 and 480 MeV
+
+  
+  if (fChain == 0) return;
+
+  Long64_t nentries = fChain->GetEntriesFast();
+
+  Long64_t nbytes = 0, nb = 0;
+  for (Long64_t jentry=0; jentry<nentries;jentry++) {
+    Long64_t ientry = LoadTree(jentry);
+    if (ientry < 0) break;
+    nb = fChain->GetEntry(jentry);   nbytes += nb;
+    // if (Cut(ientry) < 0) continue;
+    //    cout << jentry << endl;
+
+    h_endist_0->Fill(ENDIST[0]); // edep by ionisation
+    h_endist_1->Fill(ENDIST[1]); // edep by pi0, e-, e+ and #gamma
+    h_endist_2->Fill(ENDIST[2]); // edep by nuclear recoils and heavy fragments
+    h_endist_3->Fill(ENDIST[3]); // edep by part. #lt threshold
+    h_endist_4->Fill(ENDIST[4]); // E leaving the syste
+    h_endist_5->Fill(ENDIST[5]); // E carried by discarded particles
+    h_endist_6->Fill(ENDIST[6]); // resid, excit. E after evap.
+    h_endist_7->Fill(ENDIST[7]); // edep by low-energy neutrons
+    h_endist_8->Fill(ENDIST[8]); // E of part. #gt time limit
+    h_endist_9->Fill(ENDIST[9]); // E lost in endothermic nuclear reactions
+    h_endist_10->Fill(ENDIST[10]); // E lost in endothermic low-E n-reaction
+    h_endist_11->Fill(ENDIST[11]); // missing E
+
+    etotal = 0.;
+    for (unsigned int i = 0; i<3; i++ ){
+      etotal += ENDIST[i];
+    }
+
+    ratio_electron = ENDIST[1] / etotal;
+    ratio_primary = ENDIST[0] / etotal;
+
+    h_ratio_e->Fill(ratio_electron);
+    h_ratio_p->Fill(ratio_primary);
+    h_sum_ep->Fill(ENDIST[1]+ENDIST[0]);
+
+    h_let_ep->Fill( (ENDIST[1]+ENDIST[0])/scale/64. ); // [GeV/cm/Z^2]
+    h_let_e->Fill( (ENDIST[1])/scale/64. );
+    h_let_p->Fill( (ENDIST[0])/scale/64. );
+
+    spratio = ( ( (ENDIST[0]+ENDIST[1])/scale*1000./1.06) /  f_o16sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000./16. ) ); //ratio of deposited energy to PSTAR stopping power (MeV/cm/density) / MeVcm2/g
+    h_spratio->Fill(spratio);
+
+    //   cout <<  ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06)  << " " << f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) << " " << spratio << endl;
+
+  }
+}
+
+
+
+int main(int argc, const char* argv[])
+{
+  Int_t thickkey = atoi(argv[2]);
+  Float_t THICKNESS[6] = {0.025, 0.050, 0.100, 0.150, 0.500, 1.000}; //#cm
+  // Float_t THICKNESS[5] = {0.050, 0.100, 0.150, 1.000, 10.00}; //#cm
+
+  ///open the input and output root files
+  char  finname[50];
+  sprintf(finname, "%s.root",argv[1]);
+  TFile * f = new TFile(finname,"OPEN");
+
+  char  foutname[50];
+  sprintf(foutname, "%s_out.root",argv[1]);
+  cout << foutname << endl;
+  TFile * fout = new TFile(foutname, "RECREATE");
+
+  //run the analysis loop above
+  eventdatroot t;
+  TTree * tree;
+  f->GetObject("EVENTDAT",tree);  
+  t.Init(tree);
+  t.Loop( THICKNESS[thickkey] );
+  t.Closefile(fout);
+
+  return 1;
+
+}
diff --git a/oxygen/eventdatroot.h b/oxygen/eventdatroot.h
new file mode 100644
index 0000000..f4cb35a
--- /dev/null
+++ b/oxygen/eventdatroot.h
@@ -0,0 +1,199 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Mar  8 13:05:30 2019 by ROOT version 6.12/06
+// from TTree EVENTDAT/FLUKA Course Exercise	DATE:  3/ 4/19,  TIME: 11:54:46
+// found on file: jobs/runjob1001_eventdata.root
+//////////////////////////////////////////////////////////
+
+#ifndef eventdatroot_h
+#define eventdatroot_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1.h>
+
+// Header file for the classes stored in the TTree if any.
+
+class eventdatroot {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Float_t         DATA_ALL_PART[6];
+   Float_t         DATA_BEAMPART[6];
+   Float_t         DATA_ENERGY[6];
+   Float_t         DATA_id211[6];
+   Int_t           SEEDS[4];
+   Float_t         ENDIST[12];
+
+   // List of branches
+   TBranch        *b_DATA;   //!
+   TBranch        *b_seed;   //!
+   TBranch        *b_endist;   //!
+
+      
+   TH1F * h_endist_0;
+   TH1F * h_endist_1;
+   TH1F * h_endist_2;
+   TH1F * h_endist_3;
+   TH1F * h_endist_4;
+   TH1F * h_endist_5;
+   TH1F * h_endist_6;
+   TH1F * h_endist_7;
+   TH1F * h_endist_8;
+   TH1F * h_endist_9;
+   TH1F * h_endist_10;
+   TH1F * h_endist_11;
+
+   TH1F * h_ratio_e;
+   TH1F * h_ratio_p;
+   TH1F * h_sum_ep;
+   
+   TH1F * h_let_ep;
+   TH1F * h_let_e;
+   TH1F * h_let_p;
+
+   TH1F * h_spratio;
+
+   eventdatroot(TTree *tree=0);
+   virtual ~eventdatroot();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop(Float_t scale);
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+   virtual Int_t    Closefile(TFile * file);
+};
+
+#endif
+
+#ifdef eventdatroot_cxx
+eventdatroot::eventdatroot(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+  //  if (tree == 0) {
+     //   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("jobs/runjob1001_eventdata.root");
+      //     if (!f || !f->IsOpen()) {
+		//      f = new TFile("jobs/runjob1001_eventdata.root");
+	 //     }
+      //    f->GetObject("EVENTDAT",tree);
+
+      //   }
+   //   Init(tree);
+}
+
+eventdatroot::~eventdatroot()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t eventdatroot::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t eventdatroot::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void eventdatroot::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("DATA", DATA_ALL_PART, &b_DATA);
+   fChain->SetBranchAddress("SEEDS", SEEDS, &b_seed);
+   fChain->SetBranchAddress("ENDIST", ENDIST, &b_endist);
+
+   h_endist_0 = new TH1F("h_endist_0","edep by ionisation",400,0.0,0.05);
+   h_endist_1 = new TH1F("h_endist_1","edep by pi0, e-, e+ and #gamma",400,0.0,0.03);
+   h_endist_2 = new TH1F("h_endist_2","edep by nuclear recoils and heavy fragments",400,0.0,0.001);
+   h_endist_3 = new TH1F("h_endist_3","edep by part. #lt threshold",200,0.0,0.001);
+   h_endist_4 = new TH1F("h_endist_4","E leaving the system",200,0.0,1.1);
+   h_endist_5 = new TH1F("h_endist_5","E carried by discarded particles",200,0.0,0.1);
+   h_endist_6 = new TH1F("h_endist_6","resid, excit. E after evap.",200,0.0,0.1);
+   h_endist_7 = new TH1F("h_endist_7","edep by low-energy neutrons",200,0.0,0.1);
+   h_endist_8 = new TH1F("h_endist_8","E of part. #gt time limit",200,0.0,0.1);
+   h_endist_9 = new TH1F("h_endist_9","E lost in endothermic nuclear reactions",200,0.0,0.1);
+   h_endist_10 = new TH1F("h_endist_10","E lost in endothermic low-E n-reactions",200,0.0,0.1);
+   h_endist_11 = new TH1F("h_endist_11","missing E",400,0.0,0.01);
+
+   h_ratio_e = new TH1F("h_ratio_e","fraction of edep by pi0, e-, e+ and #gamma  ",400,0.0,1.0);
+   h_ratio_p = new TH1F("h_ratio_p","fraction of edep by ionisation",400,0.0,1.0);
+   h_sum_ep = new TH1F("h_sum_ep","sum of edep by ion, pi0, e-, e+ and #gamma,",400,0.0,0.1);
+
+   h_let_ep = new TH1F("h_let_ep","LET/Z^{2} (GeV/cm) by ion, pi0, e-, e+ and #gamma,",200,0.0,0.05);
+   h_let_e = new TH1F("h_let_e","LET/Z^{2} (GeV/cm) by pi0, e-, e+ and #gamma,",200,0.0,0.02);
+   h_let_p = new TH1F("h_let_p","LET/Z^{2} (GeV/cm) by ion",200,0.0,0.02);
+
+   h_spratio = new TH1F("h_spratio","total LET/SP (MeV/cm/#rho)/SP(MeVcm2/g)",200,0.5,3.0);
+
+   
+   Notify();
+}
+
+Bool_t eventdatroot::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void eventdatroot::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t eventdatroot::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+
+Int_t eventdatroot::Closefile(TFile * file)
+{
+  file->cd();
+  file->Write();
+  file->Close();
+
+  return 1;
+}
+
+#endif // #ifdef eventdatroot_cxx
diff --git a/oxygen/langaus.C b/oxygen/langaus.C
new file mode 100644
index 0000000..7fd72ae
--- /dev/null
+++ b/oxygen/langaus.C
@@ -0,0 +1,348 @@
+ /// \ingroup tutorial_fit
+ /// \notebook
+ /// Convoluted Landau and Gaussian Fitting Function
+ ///         (using ROOT's Landau and Gauss functions)
+ ///
+ ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
+ ///
+ ///  to execute this example, do:
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C
+ /// ~~~
+ ///
+ /// or
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C++
+ /// ~~~
+ ///
+ /// \macro_image
+ /// \macro_output
+ /// \macro_code
+ ///
+ /// \authors H.Pernegger, Markus Friedl
+ 
+ #include "TH1.h"
+ #include "TF1.h"
+ #include "TROOT.h"
+ #include "TStyle.h"
+ #include "TMath.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+ 
+ Double_t langaufun(Double_t *x, Double_t *par) {
+ 
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+ 
+       // Numeric constants
+       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+       Double_t mpshift  = -0.22278298;       // Landau maximum location
+ 
+       // Control constants
+       Double_t np = 200.0;      // number of convolution steps
+       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
+ 
+       // Variables
+       Double_t xx;
+       Double_t mpc;
+       Double_t fland;
+       Double_t sum = 0.0;
+       Double_t xlow,xupp;
+       Double_t step;
+       Double_t i;
+ 
+ 
+       // MP shift correction
+       mpc = par[1] - mpshift * par[0];
+ 
+       // Range of convolution integral
+       xlow = x[0] - sc * par[3];
+       xupp = x[0] + 2*sc * par[3];
+ 
+       step = (xupp-xlow) / np;
+ 
+       // Convolution integral of Landau and Gaussian by sum
+       for(i=1.0; i<=np/2; i++) {
+          xx = xlow + (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+ 
+          xx = xupp - (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+       }
+ 
+       return (par[2] * step * sum * invsq2pi / par[3]);
+ }
+ 
+ 
+ 
+ TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
+ {
+    // Once again, here are the Landau * Gaussian parameters:
+    //   par[0]=Width (scale) parameter of Landau density
+    //   par[1]=Most Probable (MP, location) parameter of Landau density
+    //   par[2]=Total area (integral -inf to inf, normalization constant)
+    //   par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    // Variables for langaufit call:
+    //   his             histogram to fit
+    //   fitrange[2]     lo and hi boundaries of fit range
+    //   startvalues[4]  reasonable start values for the fit
+    //   parlimitslo[4]  lower parameter limits
+    //   parlimitshi[4]  upper parameter limits
+    //   fitparams[4]    returns the final fit parameters
+    //   fiterrors[4]    returns the final fit errors
+    //   ChiSqr          returns the chi square
+    //   NDF             returns ndf
+ 
+    Int_t i;
+    Char_t FunName[100];
+ 
+    sprintf(FunName,"Fitfcn_%s",his->GetName());
+ 
+    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
+    if (ffitold) delete ffitold;
+ 
+    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
+    ffit->SetParameters(startvalues);
+    ffit->SetParNames("Width","MP","Area","GSigma");
+ 
+    for (i=0; i<4; i++) {
+       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
+    }
+ 
+    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
+ 
+    ffit->GetParameters(fitparams);    // obtain fit parameters
+    for (i=0; i<4; i++) {
+       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
+    }
+    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
+    NDF[0] = ffit->GetNDF();           // obtain ndf
+ 
+    return (ffit);              // return fit function
+ 
+ }
+ 
+ 
+ Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
+ 
+    // Seaches for the location (x value) at the maximum of the
+    // Landau-Gaussian convolute and its full width at half-maximum.
+    //
+    // The search is probably not very efficient, but it's a first try.
+ 
+    Double_t p,x,fy,fxr,fxl;
+    Double_t step;
+    Double_t l,lold;
+    Int_t i = 0;
+    Int_t MAXCALLS = 10000;
+ 
+ 
+    // Search for maximum
+ 
+    p = params[1] - 0.1 * params[0];
+    step = 0.05 * params[0];
+    lold = -2.0;
+    l    = -1.0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = langaufun(&x,params);
+ 
+       if (l < lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-1);
+ 
+    maxx = x;
+ 
+    fy = l/2;
+ 
+ 
+    // Search for right x location of fy
+ 
+    p = maxx + params[0];
+    step = params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-2);
+ 
+    fxr = x;
+ 
+ 
+    // Search for left x location of fy
+ 
+    p = maxx - 0.5 * params[0];
+    step = -params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-3);
+ 
+ 
+    fxl = x;
+ 
+    FWHM = fxr - fxl;
+    return (0);
+ }
+ 
+void langaus() {
+    // Fill Histogram
+   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
+                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
+                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
+                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
+                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
+    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
+   */
+   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
+   TH1D * hSNR;
+
+    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
+    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+   
+    
+    // hSNR = h_beamSignal_b0[13];
+    // Fitting SNR histo
+    printf("Fitting...\n");
+    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
+    
+    // Setting fit range and start values
+    Double_t fr[2];
+    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
+    Double_t chisqr;
+    Int_t    ndf;
+    TF1 *fitsnr;
+    Double_t SNRPeak, SNRFWHM;
+    char rootfilename[50] = "";
+    Double_t norm;
+    ofstream myfile;
+    myfile.open ("MPVcorrection_oxygen0mm05.txt");
+
+    
+    int j = 1;
+    for (int i = 0; i<26;i++){
+      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile *rootFile = new TFile(rootfilename,"OPEN");
+      if(!rootFile->GetListOfKeys()->Contains("h_spratio")) continue;
+      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
+      norm = hSNR->GetEntries();
+      hSNR->Scale(1/norm);
+      
+      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
+      fr[1]=hSNR->GetMean() + 3*hSNR->GetRMS();
+ 
+      pllo[0]=0.0006; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01;
+      plhi[0]=0.01; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.05;
+      sv[0]=0.002;
+      //sv[1]=2.1;
+      sv[2]=0.01;
+      //   sv[3]=0.03;
+      sv[3] = hSNR->GetRMS()/2.;
+      sv[1] = hSNR->GetMean();
+    
+      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
+
+      
+      langaupro(fp,SNRPeak,SNRFWHM);
+      graph_x[i] = beta_oxygen[i];
+      graph_y[i] = fp[1];
+      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
+      
+      printf("Fitting done\nPlotting results...\n");
+            myfile << i << " " << graph_y[i] << endl;
+
+      // Global style settings
+      gStyle->SetOptStat(1111);
+      gStyle->SetOptFit(111);
+      //gStyle->SetLabelSize(0.03,"x");
+      //gStyle->SetLabelSize(0.03,"y");
+      
+      // hSNR->GetXaxis()->SetRange(0,70);
+      hSNR->Draw();
+      fitsnr->Draw("lsame");
+      c1->Update();
+      c1->WaitPrimitive();
+      hSNR->Delete();
+      rootFile->Close();
+    }
+     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+   
+    graph_1->Draw("A*");
+    c2->SaveAs("oxygen0mm5_SPratio.pdf");
+    myfile.close();
+
+    
+ }
diff --git a/oxygen/oxygen0mm5_SPratio.pdf b/oxygen/oxygen0mm5_SPratio.pdf
new file mode 100644
index 0000000..98edec5
Binary files /dev/null and b/oxygen/oxygen0mm5_SPratio.pdf differ
diff --git a/oxygen/plots.C b/oxygen/plots.C
new file mode 100644
index 0000000..65b2e2c
--- /dev/null
+++ b/oxygen/plots.C
@@ -0,0 +1,32 @@
+{
+
+  TGraph * gr_letratio = new TGraph();
+
+    char  finname[50];
+    int j = 0;
+    
+    Double_t x, q;
+    q = 0.5; // 0.5 for "median"
+    
+    
+    //  for (int j = 0; j<5;j++){
+    for (int i = 0; i<26;i++){
+
+      sprintf(finname, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile * f = new TFile(finname,"OPEN");
+      TH1F * h = (TH1F*)f->Get("h_spratio");
+      cout << i << " " << h->GetMean() << " +/- " << h->GetRMS() << endl;
+      h->ComputeIntegral(); // just a precaution
+      h->GetQuantiles(1, &x, &q);
+      std::cout << "median = " << x << std::endl;
+
+      
+      h->Delete();
+      f->Close();
+      
+    }
+    //  }
+  
+
+}
+
diff --git a/oxygen/runJob.sh b/oxygen/runJob.sh
new file mode 100755
index 0000000..4940256
--- /dev/null
+++ b/oxygen/runJob.sh
@@ -0,0 +1,65 @@
+#!/bin/bash
+
+TYPE="HEAVYION"
+#energy=(0.04812 0.05982 0.07003 0.07917 0.08753 0.09530 0.10261 0.10956 0.11620 0.12257 0.12872 0.13465 0.14041 0.14601 0.15150 0.15689 0.16221 0.16746 0.17262 0.17770 0.18481 0.19187 0.19876 0.20548 0.21206 0.21851); #GeV/u, protons
+#energy=(0.08883 0.11058 0.12979 0.14713 0.16309 0.17801 0.19213 0.20560 0.21852 0.23098 0.24303 0.25471 0.26608 0.27719 0.28810 0.29887 0.30952 0.32007 0.33048 0.34077 0.35522 0.36964 0.38378 0.39766 0.41132 0.42477); #GeV/u, carbon
+#energy=(0.05057 0.06186 0.07173 0.08064 0.08885 0.09652 0.10376 0.11064 0.11723 0.12355 0.12964 0.13555 0.14127 0.14684 0.15226 0.15756 0.16273 0.16780 0.17277 0.17764 0.18456 0.19154 0.19836 0.20503 0.21157 0.21798); #GeV/u, helium
+energy=(0.10377 0.12965 0.15242 0.17298 0.19196 0.20963 0.22653 0.24258 0.25811 0.27304 0.2875 0.30155 0.31532 0.32861 0.34184 0.35488 0.36794 0.38048 0.39309 0.40565 0.42323 0.44082 0.45811 0.47528 0.49193 0.50838 ); #GeV/u, oxygen
+
+#THICKNESS=(0.025 0.050 0.075 0.100 1.25 0.150 0.2500 0.500) #cm
+THICKNESS=(0.025 0.050 0.100 0.150 0.500 1.000 0.125) #cm
+
+HOME=/work/leverington/fluka_gfortran7/project/oxygen
+ZED=8.
+ATN=16.
+
+for j in 0 1 2 3 4 6
+do 
+    mkdir -p $HOME/jobs$j #make the directory if it doesn't exist
+    JOB_HOME=$HOME/jobs$j 
+    rm -r  $JOB_HOME/* #clean up the directory
+    cd $JOB_HOME
+    
+    
+    for i in {0..25} #26
+    do
+	# make a copy of the template
+	cp $HOME/runjobtemplate.inp  $JOB_HOME/runjob$i.inp
+	touch $JOB_HOME/runjob$i.inp
+	
+	#replace template placeholders with values
+	sed -i 's/ENERGYIN/'"-${energy[i]}"'/g' $JOB_HOME/runjob$i.inp #the negative sign indicates energy, not momentum
+	sed -i 's/TYPE/'"$TYPE"'/g' $JOB_HOME/runjob$i.inp  #beam particle type
+	sed -i 's/THICK/'"${THICKNESS[j]}"'/g' $JOB_HOME/runjob$i.inp  #thickness of polystyrene layer
+	sed -i 's/ZED/'"$ZED"'/g' $JOB_HOME/runjob$i.inp  #atomic charge
+	sed -i 's/ATN/'"$ATN"'/g' $JOB_HOME/runjob$i.inp  #atomic number
+
+	echo "Job: $j-$i   Ion type: $TYPE    Energy: ${energy[i]} GeV/u    Plane thickness: ${THICKNESS[j]} cm"
+	echo ""
+	
+	touch $JOB_HOME/fluka$j-$i.sh  #script to submit to the queue
+	echo "#!/bin/bash" >> $JOB_HOME/fluka$j-$i.sh
+	echo ". /local/env.sh" >> $JOB_HOME/fluka$j-$i.sh
+#	echo ". /cvmfs/lhcb.cern.ch/lib/lhcb/LBSCRIPTS/LBSCRIPTS_v9r2p2/InstallArea/scripts/LbLogin.sh -c x86_64-slc6-gcc7-opt" >> $JOB_HOME/fluka$j-$i.sh
+	echo "source /home/lhcb/leverington/setup.sh" >> $JOB_HOME/fluka$j-$i.sh
+        
+        #execute this command(s)
+	echo "cd "$JOB_HOME >> $JOB_HOME/fluka$j-$i.sh
+	echo "$FLUPRO/flutil/rfluka -e $FLUPRO/flutil/flukadpm3 -N0 -M1 runjob"$i >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "/work/leverington/readfluka/tools/eventdat2root runjob"$i"001_eventdata"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "$HOME/eventdatroot runjob"$i"001_eventdata "$j  >> $JOB_HOME/fluka$j-$i.sh
+	
+#	qsub -l os=slc6  -l ujl=20 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+	qsub -l os=slc6 -l hio=10 -l ujl=40 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+
+#	$FLUPRO/flutil/rfluka -e /home/leverington/software/fluka/flutil/flukadpm3 -N0 -M1 $JOB_HOME/runjob$i  #run 
+#	sleep 2
+#	eventdat2root runjob$i\001_eventdata  #convert to root file
+#	sleep 2
+#	$HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+ 	sleep 1
+    done
+    sleep 30
+done
diff --git a/oxygen/runanalysis.sh b/oxygen/runanalysis.sh
new file mode 100755
index 0000000..5deb6de
--- /dev/null
+++ b/oxygen/runanalysis.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+HOME=/work/leverington/fluka_gfortran7/project/oxygen
+
+
+for j in {0..5}
+do
+    JOB_HOME=$HOME/jobs$j
+    cd $JOB_HOME
+
+	     for i in {0..25} #26
+	     do
+		 
+		 $HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+		 sleep 0.1 
+	     done
+done
diff --git a/oxygen/runjobtemplate.inp b/oxygen/runjobtemplate.inp
new file mode 100644
index 0000000..806aad9
--- /dev/null
+++ b/oxygen/runjobtemplate.inp
@@ -0,0 +1,95 @@
+* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
+TITLE
+FLUKA Course Exercise
+*
+* use names everywhere and free format for geometry
+DEFAULTS                                                              HADROTHE
+*
+* beam definitions
+BEAM      ENERGYIN  0.0       -1.7      -0.8      -0.8                TYPE
+HI-PROPE     ZED       ATN
+BEAMPOS                             -0.1
+DELTARAY       1E-06       0.0       0.0  HYDROGEN  @LASTMAT          NOPRINT
+PART-THR  -0.00001  PROTON    PROTON                 0.0   
+*
+* Geometry
+* --------
+GEOBEGIN                                                              COMBNAME
+    0    0             Cylindrical Target
+*
+* Bodies
+* ------
+*
+* Blackhole to include geometry
+SPH BLK        0.0 0.0 0.0 10000.
+* Void sphere
+SPH VOID       0.0 0.0 0.0 1000.
+* Infinite cylinder
+ZCC TARG       0.0 0.0 5.
+*  planes cutting the cylinder
+XYP ZTlow      0.0
+XYP T1seg      THICK
+XYP T2seg      2.
+XYP ZThigh     10.
+*  additional plane for scoring
+XYP ZTscor     9.9999
+END
+*
+* Regions
+* -------
+*
+* Blackhole
+BLKHOLE      5 +BLK -VOID
+*
+* Target segment 1
+TARGS1       5 +TARG -ZTlow +T1seg
+* Target segment 2
+TARGS2       5 +TARG -T1seg +T2seg
+* Target segment 3
+TARGS3       5 +TARG -T2seg +ZTscor
+* Scoring region
+TARGS4       5 +TARG -ZTscor +ZThigh
+*
+* Air around target
+INAIR        5 | +VOID -TARG
+               | +VOID +ZTlow
+               | +VOID -ZThigh
+END
+*
+GEOEND
+*
+* Materials definition
+* --------------------
+MATERIAL                        0.001965                              CO2
+COMPOUND          1.    CARBON        2.    OXYGEN                    CO2
+*
+* Assign materials
+* ----------------
+ASSIGNMA    BLCKHOLE   BLKHOLE
+ASSIGNMA    POLYSTYR    TARGS1
+ASSIGNMA      VACUUM    TARGS2
+ASSIGNMA      VACUUM    TARGS3
+ASSIGNMA      VACUUM    TARGS4
+ASSIGNMA      VACUUM     INAIR
+*
+* Exercise: heavy ions
+* --------------------
+*
+*  charge spectrum of ions
+USRYIELD       1422.  HEAVYION      -68.    TARGS1    TARGS2        1.FragZ1
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS2    TARGS3        1.FragZ2
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS3    TARGS4        1.FragZ3
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+*
+*  LET in water of ions and charged particles
+USRYIELD       2223.  HEAVYION      -69.    TARGS3    TARGS4        1.LETHI
+USRYIELD         20.       0.0      200.       9.5       2.5     2703. &
+USRYIELD       2223.  ALL-CHAR      -69.    TARGS3    TARGS4        1.LETCh
+USRYIELD         20.       0.0      200.       9.5      -2.5     2703. &
+EVENTDAT        -21.                                                  eventdata
+*
+RANDOMIZ          1.
+START        100000.                                     0.0
+STOP
diff --git a/pion/MPVcorrection_proton0mm05 b/pion/MPVcorrection_proton0mm05
new file mode 100644
index 0000000..6e11235
--- /dev/null
+++ b/pion/MPVcorrection_proton0mm05
@@ -0,0 +1,116 @@
+%!PS-Adobe-2.0
+%%Title: MPVcorrection_proton0mm05: c2
+%%Creator: ROOT Version 6.12/06
+%%CreationDate: Mon Apr 15 17:01:35 2019
+%%Orientation: Landscape
+%%DocumentNeededResources: ProcSet (FontSetInit)
+%%EndComments
+%%BeginProlog
+/s {stroke} def /l {lineto} def /m {moveto} def /t {translate} def
+/r {rotate} def /rl {roll}  def /R {repeat} def
+/d {rlineto} def /rm {rmoveto} def /gr {grestore} def /f {eofill} def
+/c {setrgbcolor} def /black {0 setgray} def /sd {setdash} def
+/cl {closepath} def /sf {scalefont setfont} def /lw {setlinewidth} def
+/box {m dup 0 exch d exch 0 d 0 exch neg d cl} def
+/NC{systemdict begin initclip end}def/C{NC box clip newpath}def
+/bl {box s} def /bf {gsave box gsave f grestore 1 lw [] 0 sd s grestore} def /Y { 0 exch d} def /X { 0 d} def 
+/K {{pop pop 0 moveto} exch kshow} bind def
+/ita {/ang 15 def gsave [1 0 ang dup sin exch cos div 1 0 0] concat} def 
+/mp {newpath /y exch def /x exch def} def
+/side {[w .77 mul w .23 mul] .385 w mul sd w 0 l currentpoint t -144 r} def
+/mr {mp x y w2 0 360 arc} def /m24 {mr s} def /m20 {mr f} def
+/mb {mp x y w2 add m w2 neg 0 d 0 w neg d w 0 d 0 w d cl} def
+/mt {mp x y w2 add m w2 neg w neg d w 0 d cl} def
+/w4 {w 4 div} def
+/w6 {w 6 div} def
+/w8 {w 8 div} def
+/m21 {mb f} def /m25 {mb s} def /m22 {mt f} def /m26{mt s} def
+/m23 {mp x y w2 sub m w2 w d w neg 0 d cl f} def
+/m27 {mp x y w2 add m w3 neg w2 neg d w3 w2 neg d w3 w2 d cl s} def
+/m28 {mp x w2 sub y w2 sub w3 add m w3 0 d  0 w3 neg d w3 0 d 0 w3 d w3 0 d  0 w3 d w3 neg 0 d 0 w3 d w3 neg 0 d 0 w3 neg d w3 neg 0 d cl s } def
+/m29 {mp gsave x w2 sub y w2 add w3 sub m currentpoint t 4 {side} repeat cl fill gr} def
+/m30 {mp gsave x w2 sub y w2 add w3 sub m currentpoint t 4 {side} repeat cl s gr} def
+/m31 {mp x y w2 sub m 0 w d x w2 sub y m w 0 d x w2 sub y w2 add m w w neg d x w2 sub y w2 sub m w w d s} def
+/m32 {mp x y w2 sub m w2 w d w neg 0 d cl s} def
+/m33 {mp x y w2 add m w3 neg w2 neg d w3 w2 neg d w3 w2 d cl f} def
+/m34 {mp x w2 sub y w2 sub w3 add m w3 0 d  0 w3 neg d w3 0 d 0 w3 d w3 0 d  0 w3 d w3 neg 0 d 0 w3 d w3 neg 0 d 0 w3 neg d w3 neg 0 d cl f } def
+/m35 {mp x y w2 add m w2 neg w2 neg d w2 w2 neg d w2 w2 d w2 neg w2 d x y w2 sub m 0 w d x w2 sub y m w 0 d s} def
+/m36 {mb x w2 sub y w2 add m w w neg d x w2 sub y w2 sub m w w d s} def
+/m37 {mp x y m w4 neg w2 d w4 neg w2 neg d w2 0 d  w4 neg w2 neg d w2 0 d w4 neg w2 d w2 0 d w4 neg w2 d w4 neg w2 neg d cl s} def
+/m38 {mp x w4 sub y w2 add m w4 neg w4 neg d 0 w2 neg d w4 w4 neg d w2 0 d w4 w4 d 0 w2 d w4 neg w4 d w2 neg 0 d x y w2 sub m 0 w d x w2 sub y m w 0 d cl s} def
+/m39 {mp x y m w4 neg w2 d w4 neg w2 neg d w2 0 d  w4 neg w2 neg d w2 0 d w4 neg w2 d w2 0 d w4 neg w2 d w4 neg w2 neg d cl f} def
+/m40 {mp x y m w4 w2 d w4 w4 neg d w2 neg w4 neg d w2 w4 neg d w4 neg w4 neg d w4 neg w2 d w4 neg w2 neg d w4 neg w4 d w2 w4 d w2 neg w4 d w4 w4 d w4 w2 neg d cl s} def
+/m41 {mp x y m w4 w2 d w4 w4 neg d w2 neg w4 neg d w2 w4 neg d w4 neg w4 neg d w4 neg w2 d w4 neg w2 neg d w4 neg w4 d w2 w4 d w2 neg w4 d w4 w4 d w4 w2 neg d cl f} def
+/m42 {mp x y w2 add m w8 neg w2 -3 4 div mul d w2 -3 4 div mul w8 neg d w2 3 4 div mul w8 neg d w8 w2 -3 4 div mul d w8 w2 3 4 div mul d w2 3 4 div mul w8 d w2 -3 4 div mul w8 d w8 neg w2 3 4 div mul d cl s} def
+/m43 {mp x y w2 add m w8 neg w2 -3 4 div mul d w2 -3 4 div mul w8 neg d w2 3 4 div mul w8 neg d w8 w2 -3 4 div mul d w8 w2 3 4 div mul d w2 3 4 div mul w8 d w2 -3 4 div mul w8 d w8 neg w2 3 4 div mul d cl f} def
+/m44 {mp x y m w6 neg w2 d w2 2 3 div mul 0 d w6 neg w2 neg d w2 w6 d 0 w2 -2 3 div mul d w2 neg w6 d w6 w2 neg d w2 -2 3 div mul 0 d w6 w2 d w2 neg w6 neg d 0 w2 2 3 div mul d w2 w6 neg d cl s} def
+/m45 {mp x y m w6 neg w2 d w2 2 3 div mul 0 d w6 neg w2 neg d w2 w6 d 0 w2 -2 3 div mul d w2 neg w6 d w6 w2 neg d w2 -2 3 div mul 0 d w6 w2 d w2 neg w6 neg d 0 w2 2 3 div mul d w2 w6 neg d cl f} def
+/m46 {mp x y w4 add m w4 neg w4 d w4 neg w4 neg d  w4 w4 neg d w4 neg w4 neg d w4 w4 neg d w4 w4 d w4 w4 neg d w4 w4 d w4 neg w4 d w4 w4 d w4 neg w4 d w4 neg w4 neg d cl s} def
+/m47 {mp x y w4 add m w4 neg w4 d w4 neg w4 neg d w4 w4 neg d w4 neg w4 neg d  w4 w4 neg d w4 w4 d w4 w4 neg d w4 w4 d w4 neg w4 d w4 w4 d w4 neg w4 d w4 neg w4 neg d cl f} def
+/m48 {mp x y w4 add m w4 neg w4 d w4 neg w4 neg d w4 w4 neg d  w4 neg w4 neg d w4 w4 neg d w4 w4 d w4 w4 neg d w4 w4 d w4 neg w4 d w4 w4 d w4 neg w4 d w4 neg w4 neg d  w4 w4 neg d w4 neg w4 neg d w4 neg w4 d w4 w4 d cl f} def
+/m49 {mp x w2 sub w3 add y w2 sub w3 add m  0 w3 neg d w3 0 d 0 w3 d w3 0 d 0 w3 d w3 neg 0 d 0 w3 d w3 neg 0 d 0 w3 neg d w3 neg 0 d 0 w3 neg d w3 0 d 0 w3 d w3 0 d 0 w3 neg d w3 neg 0 d cl f } def
+/m2 {mp x y w2 sub m 0 w d x w2 sub y m w 0 d s} def
+/m5 {mp x w2 sub y w2 sub m w w d x w2 sub y w2 add m w w neg d s} def
+%%IncludeResource: ProcSet (FontSetInit)
+%%IncludeResource: font Times-Roman
+%%IncludeResource: font Times-Italic
+%%IncludeResource: font Times-Bold
+%%IncludeResource: font Times-BoldItalic
+%%IncludeResource: font Helvetica
+%%IncludeResource: font Helvetica-Oblique
+%%IncludeResource: font Helvetica-Bold
+%%IncludeResource: font Helvetica-BoldOblique
+%%IncludeResource: font Courier
+%%IncludeResource: font Courier-Oblique
+%%IncludeResource: font Courier-Bold
+%%IncludeResource: font Courier-BoldOblique
+%%IncludeResource: font Symbol
+%%IncludeResource: font ZapfDingbats
+/reEncode {exch findfont dup length dict begin {1 index /FID eq  {pop pop} {def} ifelse } forall /Encoding exch def currentdict end dup /FontName get exch definefont pop } def [/Times-Bold /Times-Italic /Times-BoldItalic /Helvetica /Helvetica-Oblique
+ /Helvetica-Bold /Helvetica-BoldOblique /Courier /Courier-Oblique /Courier-Bold /Courier-BoldOblique /Times-Roman /AvantGarde-Book /AvantGarde-BookOblique /AvantGarde-Demi /AvantGarde-DemiOblique /Bookman-Demi /Bookman-DemiItalic /Bookman-Light
+ /Bookman-LightItalic /Helvetica-Narrow /Helvetica-Narrow-Bold /Helvetica-Narrow-BoldOblique /Helvetica-Narrow-Oblique /NewCenturySchlbk-Roman /NewCenturySchlbk-Bold /NewCenturySchlbk-BoldItalic /NewCenturySchlbk-Italic /Palatino-Bold
+ /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman ] {ISOLatin1Encoding reEncode } forall/Zone {/iy exch def /ix exch def  ix 1 sub  3144 mul  1 iy sub  2224 mul t} def
+%%EndProlog
+%%BeginSetup
+%%EndSetup
+newpath  gsave  90 r 0 -594 t  28 20 t .25 .25 scale  gsave 
+%%Page: 1 1
+ gsave  gsave 
+ 1 1 Zone
+ gsave  0 0 t black[  ] 0 sd 3 lw 1 1 1 c 2948 2117 0 0 bf black 1 1 1 c 2358 1693 295 212 bf black 2358 1693 295 212 bl 295 212 m 2358 X s 426 262 m -50 Y s 503 237 m -25 Y s 579 237 m -25 Y s 656 237 m -25 Y s 733 237 m -25 Y s 809 262 m -50 Y s
+ 886 237 m -25 Y s 963 237 m -25 Y s 1039 237 m -25 Y s 1116 237 m -25 Y s 1193 262 m -50 Y s 1269 237 m -25 Y s 1346 237 m -25 Y s 1423 237 m -25 Y s 1499 237 m -25 Y s 1576 262 m -50 Y s 1653 237 m -25 Y s 1729 237 m -25 Y s 1806 237 m -25 Y s 1883
+ 237 m -25 Y s 1959 262 m -50 Y s 2036 237 m -25 Y s 2113 237 m -25 Y s 2189 237 m -25 Y s 2266 237 m -25 Y s 2343 262 m -50 Y s 426 262 m -50 Y s 349 237 m -25 Y s 2343 262 m -50 Y s 2420 237 m -25 Y s 2496 237 m -25 Y s 2573 237 m -25 Y s 2650 237
+ m -25 Y s
+ gsave  2948 2117 0 0 C 376.816 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.3) show NC gr 
+ gsave  2948 2117 0 0 C 738.855 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.35) show NC gr 
+ gsave  2948 2117 0 0 C 1141.53 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.4) show NC gr 
+ gsave  2948 2117 0 0 C 1503.57 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.45) show NC gr 
+ gsave  2948 2117 0 0 C 1909.94 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.5) show NC gr 
+ gsave  2948 2117 0 0 C 2271.98 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.55) show NC gr  295 212 m 1693 Y s 366 256 m -71 X s 330 303 m -35 X s 330 351 m -35 X s 330 398 m -35 X s 366 445 m -71 X s 330 493 m -35 X s 330 540 m -35 X s 330
+ 588 m -35 X s 366 635 m -71 X s 330 682 m -35 X s 330 730 m -35 X s 330 777 m -35 X s 366 824 m -71 X s 330 872 m -35 X s 330 919 m -35 X s 330 967 m -35 X s 366 1014 m -71 X s 330 1061 m -35 X s 330 1109 m -35 X s 330 1156 m -35 X s 366 1203 m -71
+ X s 330 1251 m -35 X s 330 1298 m -35 X s 330 1346 m -35 X s 366 1393 m -71 X s 330 1440 m -35 X s 330 1488 m -35 X s 330 1535 m -35 X s 366 1582 m -71 X s 330 1630 m -35 X s 330 1677 m -35 X s 330 1725 m -35 X s 366 1772 m -71 X s 366 256 m -71 X s
+ 366 1772 m -71 X s 330 1819 m -35 X s 330 1867 m -35 X s
+ gsave  2948 2117 0 0 C 181.019 232.739 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.8) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 421.147 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.82) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 609.555 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.84) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 801.658 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.86) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 990.066 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.88) show NC gr 
+ gsave  2948 2117 0 0 C 181.019 1178.47 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.9) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1370.58 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.92) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1558.98 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.94) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1747.39 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.96) show NC gr  1 1 1 c black
+ gsave  2948 2117 0 0 C 1333.63 2013.38 t 0 r /Helvetica findfont 103.44 sf 0 0 m (Graph) show NC gr  1 1 1 c black /w 30 def /w2 {w 2 div} def /w3 {w 3 div} def 491 1655 740 1444 933 1271 1091 1163 1225 1077 1342 1016 1445 968 1539 903 1625 840 1703
+ 796 1777 763 1845 736 1908 711 1969 672 2026 651 2080 618 2132 577 2183 591 2231 568 2277 524 2340 514 2400 491 2457 447 23 { m31} R 491 1655 m 109 Y s 484 1764 m 15 X s 491 1655 m -109 Y s 484 1546 m 15 X s 740 1444 m 106 Y s 733 1550 m 15 X s 740
+ 1444 m -106 Y s 733 1338 m 15 X s 933 1271 m 104 Y s 926 1375 m 15 X s 933 1271 m -104 Y s 926 1167 m 15 X s 1091 1163 m 102 Y s 1084 1265 m 14 X s 1091 1163 m -103 Y s 1084 1060 m 14 X s 1225 1077 m 101 Y s 1218 1178 m 14 X s 1225 1077 m -102 Y s
+ 1218 975 m 14 X s 1342 1016 m 101 Y s 1334 1117 m 15 X s 1342 1016 m -101 Y s 1334 915 m 15 X s 1445 968 m 101 Y s 1438 1069 m 15 X s 1445 968 m -100 Y s 1438 868 m 15 X s 1539 903 m 99 Y s 1532 1002 m 15 X s 1539 903 m -100 Y s 1532 803 m 15 X s
+ 1625 840 m 99 Y s 1617 939 m 15 X s 1625 840 m -98 Y s 1617 742 m 15 X s 1703 796 m 98 Y s 1696 894 m 15 X s 1703 796 m -98 Y s 1696 698 m 15 X s 1777 763 m 97 Y s 1769 860 m 15 X s 1777 763 m -98 Y s 1769 665 m 15 X s 1845 736 m 98 Y s 1837 834 m
+ 15 X s 1845 736 m -97 Y s 1837 639 m 15 X s 1908 711 m 97 Y s 1901 808 m 15 X s 1908 711 m -97 Y s 1901 614 m 15 X s 1969 672 m 96 Y s 1961 768 m 15 X s 1969 672 m -97 Y s 1961 575 m 15 X s 2026 651 m 96 Y s 2018 747 m 15 X s 2026 651 m -97 Y s 2018
+ 554 m 15 X s 2080 618 m 96 Y s 2073 714 m 15 X s 2080 618 m -96 Y s 2073 522 m 15 X s 2132 577 m 95 Y s 2125 672 m 15 X s 2132 577 m -96 Y s 2125 481 m 15 X s 2183 591 m 95 Y s 2175 686 m 15 X s 2183 591 m -96 Y s 2175 495 m 15 X s 2231 568 m 95 Y s
+ 2223 663 m 15 X s 2231 568 m -96 Y s 2223 472 m 15 X s 2277 524 m 94 Y s 2269 618 m 15 X s 2277 524 m -95 Y s 2269 429 m 15 X s 2340 514 m 95 Y s 2332 609 m 15 X s 2340 514 m -95 Y s 2332 419 m 15 X s 2400 491 m 94 Y s 2392 585 m 15 X s 2400 491 m
+ -95 Y s 2392 396 m 15 X s 2457 447 m 94 Y s 2449 541 m 15 X s 2457 447 m -94 Y s 2449 353 m 15 X s 1 1 1 c black
+ gsave  2948 2117 0 0 C 1333.63 2013.38 t 0 r /Helvetica findfont 103.44 sf 0 0 m (Graph) show NC gr  gr showpage
+ gr 
+%%Trailer
+%%Pages:  1
+ gr  gr  gr 
+%%EOF
diff --git a/pion/MPVcorrection_proton0mm05.txt b/pion/MPVcorrection_proton0mm05.txt
new file mode 100644
index 0000000..1515740
--- /dev/null
+++ b/pion/MPVcorrection_proton0mm05.txt
@@ -0,0 +1,26 @@
+0 0.94765
+1 0.925397
+2 0.907096
+3 0.895693
+4 0.886649
+5 0.880236
+6 0.875197
+7 0.868248
+8 0.861681
+9 0.857026
+10 0.853484
+11 0.850676
+12 0.848045
+13 0.843891
+14 0.841664
+15 0.838244
+16 0.83387
+17 0.835321
+18 0.832924
+19 0.828248
+20 0.827223
+21 0.824773
+22 0.82014
+23 0.818965
+24 0.816504
+25 0.815992
diff --git a/pion/Makefile b/pion/Makefile
new file mode 100644
index 0000000..c16461b
--- /dev/null
+++ b/pion/Makefile
@@ -0,0 +1,38 @@
+
+ROOTCFLAGS   := $(shell root-config --cflags)
+ROOTLIBS     := $(shell root-config --libs)
+ROOTGLIBS    := $(shell root-config --glibs)
+
+GSLCFLAGS   := $(shell gsl-config --cflags)
+GSLLIBS     := $(shell gsl-config --libs)
+GSLGLIBS    := $(shell gsl-config --glibs)
+
+
+LDFLAGS       = -O
+LIBS         += $(ROOTLIBS) $(GSLLIBS) -L/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/lib
+CFLAGS       += $(ROOTCFLAGS) $(GSLCFLAGS)  -I/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/include
+
+#  Not sure why Minuit isn't being included -- put in by hand
+#
+LIBS         += -lMinuit -ltbb
+
+
+all: eventdatroot
+
+
+
+eventdatroot: eventdatroot.o
+	g++ -o eventdatroot  eventdatroot.o $(LDFLAGS) $(LIBS)
+
+%.o: %.c
+	g++ ${CFLAGS} -c  -g $<
+
+
+test:
+	@echo $(ROOTCFLAGS) 
+	@echo $(LDFLAGS) 
+	@echo $(LIBS)
+
+clean:	
+	-rm -f  *.o
+
diff --git a/pion/eventdatroot b/pion/eventdatroot
new file mode 100755
index 0000000..e65d083
Binary files /dev/null and b/pion/eventdatroot differ
diff --git a/pion/eventdatroot.c b/pion/eventdatroot.c
new file mode 100644
index 0000000..a13394e
--- /dev/null
+++ b/pion/eventdatroot.c
@@ -0,0 +1,135 @@
+#define eventdatroot_cxx
+#include "eventdatroot.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+//#include <TROOT.h>
+#include <TF1.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+
+
+using namespace std;
+
+void eventdatroot::Loop(Float_t scale)
+{
+  //   In a ROOT session, you can do:
+  //      root> .L eventdatroot.C
+  //      root> eventdatroot t
+  //      root> t.GetEntry(12); // Fill t data members with entry number 12
+  //      root> t.Show();       // Show values of entry 12
+  //      root> t.Show(16);     // Read and show values of entry 16
+  //      root> t.Loop();       // Loop on all entries
+  //
+
+  //     This is the loop skeleton where:
+  //    jentry is the global entry number in the chain
+  //    ientry is the entry number in the current Tree
+  //  Note that the argument to GetEntry must be:
+  //    jentry for TChain::GetEntry
+  //    ientry for TTree::GetEntry and TBranch::GetEntry
+  //
+  //       To read only selected branches, Insert statements like:
+  // METHOD1:
+  //    fChain->SetBranchStatus("*",0);  // disable all branches
+  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
+  // METHOD2: replace line
+  //    fChain->GetEntry(jentry);       //read all branches
+  //by  b_branchname->GetEntry(ientry); //read only this branch
+
+  Float_t ratio_electron, ratio_primary, etotal, spratio;
+  TF1 * f_sp_ps = new TF1("f_sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in polystyrene [MeV cm2 /g]
+  f_sp_ps->SetParameters(361.936, -0.892255, 1.19133);  //protons between 50 and 220 MeV
+
+  
+  if (fChain == 0) return;
+
+  Long64_t nentries = fChain->GetEntriesFast();
+
+  Long64_t nbytes = 0, nb = 0;
+  for (Long64_t jentry=0; jentry<nentries;jentry++) {
+    Long64_t ientry = LoadTree(jentry);
+    if (ientry < 0) break;
+    nb = fChain->GetEntry(jentry);   nbytes += nb;
+    // if (Cut(ientry) < 0) continue;
+    //    cout << jentry << endl;
+
+    h_endist_0->Fill(ENDIST[0]); // edep by ionisation
+    h_endist_1->Fill(ENDIST[1]); // edep by pi0, e-, e+ and #gamma
+    h_endist_2->Fill(ENDIST[2]); // edep by nuclear recoils and heavy fragments
+    h_endist_3->Fill(ENDIST[3]); // edep by part. #lt threshold
+    h_endist_4->Fill(ENDIST[4]); // E leaving the syste
+    h_endist_5->Fill(ENDIST[5]); // E carried by discarded particles
+    h_endist_6->Fill(ENDIST[6]); // resid, excit. E after evap.
+    h_endist_7->Fill(ENDIST[7]); // edep by low-energy neutrons
+    h_endist_8->Fill(ENDIST[8]); // E of part. #gt time limit
+    h_endist_9->Fill(ENDIST[9]); // E lost in endothermic nuclear reactions
+    h_endist_10->Fill(ENDIST[10]); // E lost in endothermic low-E n-reaction
+    h_endist_11->Fill(ENDIST[11]); // missing E
+
+    etotal = 0.;
+    for (unsigned int i = 0; i<3; i++ ){
+      etotal += ENDIST[i];
+    }
+
+    ratio_electron = ENDIST[1] / etotal;
+    ratio_primary = ENDIST[0] / etotal;
+
+    h_ratio_e->Fill(ratio_electron);
+    h_ratio_p->Fill(ratio_primary);
+    h_sum_ep->Fill(ENDIST[1]+ENDIST[0]);
+
+    h_let_ep->Fill( (ENDIST[1]+ENDIST[0])/scale );
+    h_let_e->Fill( (ENDIST[1])/scale );
+    h_let_p->Fill( (ENDIST[0])/scale );
+
+    spratio = ( ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06) / f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) ); //ratio of deposited energy to PSTAR stopping power (MeV/cm/density) / MeVcm2/g
+    h_spratio->Fill(spratio);
+
+    //   cout <<  ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06)  << " " << f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) << " " << spratio << endl;
+
+  }
+}
+
+
+
+int main(int argc, const char* argv[])
+{
+  Int_t thickkey = atoi(argv[2]);
+  Float_t THICKNESS[6] = {0.025, 0.050, 0.100, 0.150, 0.500, 1.000}; //#cm
+  // Float_t THICKNESS[5] = {0.050, 0.100, 0.150, 1.000, 10.00}; //#cm
+
+  ///open the input and output root files
+  char  finname[50];
+  sprintf(finname, "%s.root",argv[1]);
+  TFile * f = new TFile(finname,"OPEN");
+
+  char  foutname[50];
+  sprintf(foutname, "%s_out.root",argv[1]);
+  cout << foutname << endl;
+  TFile * fout = new TFile(foutname, "RECREATE");
+
+  //run the analysis loop above
+  eventdatroot t;
+  TTree * tree;
+  f->GetObject("EVENTDAT",tree);  
+  t.Init(tree);
+  t.Loop( THICKNESS[thickkey] );
+  t.Closefile(fout);
+
+  return 1;
+
+}
diff --git a/pion/eventdatroot.h b/pion/eventdatroot.h
new file mode 100644
index 0000000..0252c6c
--- /dev/null
+++ b/pion/eventdatroot.h
@@ -0,0 +1,199 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Mar  8 13:05:30 2019 by ROOT version 6.12/06
+// from TTree EVENTDAT/FLUKA Course Exercise	DATE:  3/ 4/19,  TIME: 11:54:46
+// found on file: jobs/runjob1001_eventdata.root
+//////////////////////////////////////////////////////////
+
+#ifndef eventdatroot_h
+#define eventdatroot_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1.h>
+
+// Header file for the classes stored in the TTree if any.
+
+class eventdatroot {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Float_t         DATA_ALL_PART[6];
+   Float_t         DATA_BEAMPART[6];
+   Float_t         DATA_ENERGY[6];
+   Float_t         DATA_id211[6];
+   Int_t           SEEDS[4];
+   Float_t         ENDIST[12];
+
+   // List of branches
+   TBranch        *b_DATA;   //!
+   TBranch        *b_seed;   //!
+   TBranch        *b_endist;   //!
+
+      
+   TH1F * h_endist_0;
+   TH1F * h_endist_1;
+   TH1F * h_endist_2;
+   TH1F * h_endist_3;
+   TH1F * h_endist_4;
+   TH1F * h_endist_5;
+   TH1F * h_endist_6;
+   TH1F * h_endist_7;
+   TH1F * h_endist_8;
+   TH1F * h_endist_9;
+   TH1F * h_endist_10;
+   TH1F * h_endist_11;
+
+   TH1F * h_ratio_e;
+   TH1F * h_ratio_p;
+   TH1F * h_sum_ep;
+   
+   TH1F * h_let_ep;
+   TH1F * h_let_e;
+   TH1F * h_let_p;
+
+   TH1F * h_spratio;
+
+   eventdatroot(TTree *tree=0);
+   virtual ~eventdatroot();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop(Float_t scale);
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+   virtual Int_t    Closefile(TFile * file);
+};
+
+#endif
+
+#ifdef eventdatroot_cxx
+eventdatroot::eventdatroot(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+  //  if (tree == 0) {
+     //   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("jobs/runjob1001_eventdata.root");
+      //     if (!f || !f->IsOpen()) {
+		//      f = new TFile("jobs/runjob1001_eventdata.root");
+	 //     }
+      //    f->GetObject("EVENTDAT",tree);
+
+      //   }
+   //   Init(tree);
+}
+
+eventdatroot::~eventdatroot()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t eventdatroot::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t eventdatroot::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void eventdatroot::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("DATA", DATA_ALL_PART, &b_DATA);
+   fChain->SetBranchAddress("SEEDS", SEEDS, &b_seed);
+   fChain->SetBranchAddress("ENDIST", ENDIST, &b_endist);
+
+   h_endist_0 = new TH1F("h_endist_0","edep by ionisation",200,0.0,0.001);
+   h_endist_1 = new TH1F("h_endist_1","edep by pi0, e-, e+ and #gamma",200,0.0,0.001);
+   h_endist_2 = new TH1F("h_endist_2","edep by nuclear recoils and heavy fragments",200,0.0,0.001);
+   h_endist_3 = new TH1F("h_endist_3","edep by part. #lt threshold",200,0.0,0.001);
+   h_endist_4 = new TH1F("h_endist_4","E leaving the system",200,0.0,0.1);
+   h_endist_5 = new TH1F("h_endist_5","E carried by discarded particles",200,0.0,0.1);
+   h_endist_6 = new TH1F("h_endist_6","resid, excit. E after evap.",200,0.0,0.1);
+   h_endist_7 = new TH1F("h_endist_7","edep by low-energy neutrons",200,0.0,0.1);
+   h_endist_8 = new TH1F("h_endist_8","E of part. #gt time limit",200,0.0,0.1);
+   h_endist_9 = new TH1F("h_endist_9","E lost in endothermic nuclear reactions",200,0.0,0.1);
+   h_endist_10 = new TH1F("h_endist_10","E lost in endothermic low-E n-reactions",200,0.0,0.1);
+   h_endist_11 = new TH1F("h_endist_11","missing E",200,0.0,0.01);
+
+   h_ratio_e = new TH1F("h_ratio_e","fraction of edep by pi0, e-, e+ and #gamma  ",200,0.0,1.0);
+   h_ratio_p = new TH1F("h_ratio_p","fraction of edep by ionisation",200,0.0,1.0);
+   h_sum_ep = new TH1F("h_sum_ep","sum of edep by ion, pi0, e-, e+ and #gamma,",200,0.0,0.0005);
+
+   h_let_ep = new TH1F("h_let_ep","LET (GeV/cm) by ion, pi0, e-, e+ and #gamma,",200,0.0,0.01);
+   h_let_e = new TH1F("h_let_e","LET (GeV/cm) by pi0, e-, e+ and #gamma,",200,0.0,0.005);
+   h_let_p = new TH1F("h_let_p","LET (GeV/cm) by ion",200,0.0,0.005);
+
+   h_spratio = new TH1F("h_spratio","total LET/SP (MeV/cm/#rho)/SP(MeVcm2/g)",200,0.5,3.0);
+
+   
+   Notify();
+}
+
+Bool_t eventdatroot::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void eventdatroot::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t eventdatroot::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+
+Int_t eventdatroot::Closefile(TFile * file)
+{
+  file->cd();
+  file->Write();
+  file->Close();
+
+  return 1;
+}
+
+#endif // #ifdef eventdatroot_cxx
diff --git a/pion/langaus.C b/pion/langaus.C
new file mode 100644
index 0000000..a15e388
--- /dev/null
+++ b/pion/langaus.C
@@ -0,0 +1,343 @@
+ /// \ingroup tutorial_fit
+ /// \notebook
+ /// Convoluted Landau and Gaussian Fitting Function
+ ///         (using ROOT's Landau and Gauss functions)
+ ///
+ ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
+ ///
+ ///  to execute this example, do:
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C
+ /// ~~~
+ ///
+ /// or
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C++
+ /// ~~~
+ ///
+ /// \macro_image
+ /// \macro_output
+ /// \macro_code
+ ///
+ /// \authors H.Pernegger, Markus Friedl
+ 
+ #include "TH1.h"
+ #include "TF1.h"
+ #include "TROOT.h"
+ #include "TStyle.h"
+ #include "TMath.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+ 
+ Double_t langaufun(Double_t *x, Double_t *par) {
+ 
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+ 
+       // Numeric constants
+       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+       Double_t mpshift  = -0.22278298;       // Landau maximum location
+ 
+       // Control constants
+       Double_t np = 200.0;      // number of convolution steps
+       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
+ 
+       // Variables
+       Double_t xx;
+       Double_t mpc;
+       Double_t fland;
+       Double_t sum = 0.0;
+       Double_t xlow,xupp;
+       Double_t step;
+       Double_t i;
+ 
+ 
+       // MP shift correction
+       mpc = par[1] - mpshift * par[0];
+ 
+       // Range of convolution integral
+       xlow = x[0] - sc * par[3];
+       xupp = x[0] + 2*sc * par[3];
+ 
+       step = (xupp-xlow) / np;
+ 
+       // Convolution integral of Landau and Gaussian by sum
+       for(i=1.0; i<=np/2; i++) {
+          xx = xlow + (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+ 
+          xx = xupp - (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+       }
+ 
+       return (par[2] * step * sum * invsq2pi / par[3]);
+ }
+ 
+ 
+ 
+ TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
+ {
+    // Once again, here are the Landau * Gaussian parameters:
+    //   par[0]=Width (scale) parameter of Landau density
+    //   par[1]=Most Probable (MP, location) parameter of Landau density
+    //   par[2]=Total area (integral -inf to inf, normalization constant)
+    //   par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    // Variables for langaufit call:
+    //   his             histogram to fit
+    //   fitrange[2]     lo and hi boundaries of fit range
+    //   startvalues[4]  reasonable start values for the fit
+    //   parlimitslo[4]  lower parameter limits
+    //   parlimitshi[4]  upper parameter limits
+    //   fitparams[4]    returns the final fit parameters
+    //   fiterrors[4]    returns the final fit errors
+    //   ChiSqr          returns the chi square
+    //   NDF             returns ndf
+ 
+    Int_t i;
+    Char_t FunName[100];
+ 
+    sprintf(FunName,"Fitfcn_%s",his->GetName());
+ 
+    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
+    if (ffitold) delete ffitold;
+ 
+    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
+    ffit->SetParameters(startvalues);
+    ffit->SetParNames("Width","MP","Area","GSigma");
+ 
+    for (i=0; i<4; i++) {
+       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
+    }
+ 
+    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
+ 
+    ffit->GetParameters(fitparams);    // obtain fit parameters
+    for (i=0; i<4; i++) {
+       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
+    }
+    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
+    NDF[0] = ffit->GetNDF();           // obtain ndf
+ 
+    return (ffit);              // return fit function
+ 
+ }
+ 
+ 
+ Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
+ 
+    // Seaches for the location (x value) at the maximum of the
+    // Landau-Gaussian convolute and its full width at half-maximum.
+    //
+    // The search is probably not very efficient, but it's a first try.
+ 
+    Double_t p,x,fy,fxr,fxl;
+    Double_t step;
+    Double_t l,lold;
+    Int_t i = 0;
+    Int_t MAXCALLS = 10000;
+ 
+ 
+    // Search for maximum
+ 
+    p = params[1] - 0.1 * params[0];
+    step = 0.05 * params[0];
+    lold = -2.0;
+    l    = -1.0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = langaufun(&x,params);
+ 
+       if (l < lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-1);
+ 
+    maxx = x;
+ 
+    fy = l/2;
+ 
+ 
+    // Search for right x location of fy
+ 
+    p = maxx + params[0];
+    step = params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-2);
+ 
+    fxr = x;
+ 
+ 
+    // Search for left x location of fy
+ 
+    p = maxx - 0.5 * params[0];
+    step = -params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-3);
+ 
+ 
+    fxl = x;
+ 
+    FWHM = fxr - fxl;
+    return (0);
+ }
+ 
+void langaus() {
+    // Fill Histogram
+   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
+                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
+                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
+                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
+                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
+    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
+   */
+   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
+   TH1D * hSNR;
+
+    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
+    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+   
+    
+    // hSNR = h_beamSignal_b0[13];
+    // Fitting SNR histo
+    printf("Fitting...\n");
+    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
+    
+    // Setting fit range and start values
+    Double_t fr[2];
+    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
+    Double_t chisqr;
+    Int_t    ndf;
+    TF1 *fitsnr;
+    Double_t SNRPeak, SNRFWHM;
+    char rootfilename[50] = "";
+    int j = 1;
+     Double_t norm;
+            ofstream myfile;
+       myfile.open ("MPVcorrection_proton0mm05.txt");
+    for (int i = 0; i<26;i++){
+
+      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile *rootFile = new TFile(rootfilename,"OPEN");
+
+      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
+      norm = hSNR->GetEntries();
+      hSNR->Scale(1/norm);
+      
+      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
+      fr[1]=hSNR->GetMean() + 3*hSNR->GetRMS();
+ 
+      pllo[0]=0.00001; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01; 
+      plhi[0]=0.1; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.1;
+      sv[0]=0.05;
+      //sv[1]=2.1;
+      sv[2]=0.01;
+      sv[3]=0.02;
+      sv[1] = hSNR->GetMean();
+    
+      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
+
+      
+      langaupro(fp,SNRPeak,SNRFWHM);
+      graph_x[i] = beta_proton[i];
+      graph_y[i] = fp[1];
+      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
+      
+      printf("Fitting done\nPlotting results...\n");
+      myfile << i << " " << graph_y[i] << endl;
+      // Global style settings
+      gStyle->SetOptStat(1111);
+      gStyle->SetOptFit(111);
+      //gStyle->SetLabelSize(0.03,"x");
+      //gStyle->SetLabelSize(0.03,"y");
+      
+      // hSNR->GetXaxis()->SetRange(0,70);
+      hSNR->Draw();
+      fitsnr->Draw("lsame");
+      c1->Update();
+      //c1->WaitPrimitive();
+      hSNR->Delete();
+      rootFile->Close();
+    }
+     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+    graph_1->Draw("A*");
+    c2->SaveAs("MPVcorrection_proton0mm05");
+    myfile.close();
+    
+ }
diff --git a/pion/plots.C b/pion/plots.C
new file mode 100644
index 0000000..b49fdcb
--- /dev/null
+++ b/pion/plots.C
@@ -0,0 +1,32 @@
+{
+
+  TGraph * gr_letratio = new TGraph();
+
+    char  finname[50];
+    int j = 1;
+    
+    Double_t x, q;
+    q = 0.5; // 0.5 for "median"
+    
+    
+    //  for (int j = 0; j<5;j++){
+    for (int i = 0; i<26;i++){
+
+      sprintf(finname, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile * f = new TFile(finname,"OPEN");
+      TH1F * h = (TH1F*)f->Get("h_spratio");
+      cout << i << " " << h->GetMean() << " +/- " << h->GetRMS() << endl;
+      h->ComputeIntegral(); // just a precaution
+      h->GetQuantiles(1, &x, &q);
+      std::cout << "median = " << x << std::endl;
+
+      
+      h->Delete();
+      f->Close();
+      
+    }
+    //  }
+  
+
+}
+
diff --git a/pion/runJob.sh b/pion/runJob.sh
new file mode 100755
index 0000000..649000b
--- /dev/null
+++ b/pion/runJob.sh
@@ -0,0 +1,64 @@
+#!/bin/bash
+
+TYPE="PION+"
+#energy=(0.04812 0.05982 0.07003 0.07917 0.08753 0.09530 0.10261 0.10956 0.11620 0.12257 0.12872 0.13465 0.14041 0.14601 0.15150 0.15689 0.16221 0.16746 0.17262 0.17770 0.18481 0.19187 0.19876 0.20548 0.21206 0.21851); #GeV/u, protons
+#energy=(0.08883 0.11058 0.12979 0.14713 0.16309 0.17801 0.19213 0.20560 0.21852 0.23098 0.24303 0.25471 0.26608 0.27719 0.28810 0.29887 0.30952 0.32007 0.33048 0.34077 0.35522 0.36964 0.38378 0.39766 0.41132 0.42477); #GeV/u, carbon
+#energy=(0.05057 0.06186 0.07173 0.08064 0.08885 0.09652 0.10376 0.11064 0.11723 0.12355 0.12964 0.13555 0.14127 0.14684 0.15226 0.15756 0.16273 0.16780 0.17277 0.17764 0.18456 0.19154 0.19836 0.20503 0.21157 0.21798); #GeV/u, helium
+energy=(180.0000 0)
+#THICKNESS=(0.025 0.050 0.075 0.100 1.25 0.150 0.2500 0.500) #cm
+THICKNESS=(0.025 0.050 0.100 0.150 0.500 1.000 0.125) #cm
+
+HOME=/work/leverington/fluka_gfortran7/project/pion
+#ZED=2.
+#ATN=4.
+
+for j in {0..6}
+do 
+    mkdir -p $HOME/jobs$j #make the directory if it doesn't exist
+    JOB_HOME=$HOME/jobs$j 
+    rm -r  $JOB_HOME/* #clean up the directory
+    cd $JOB_HOME
+    
+    
+    for i in 0 #26
+    do
+	# make a copy of the template
+	cp $HOME/runjobtemplate.inp  $JOB_HOME/runjob$i.inp
+	touch $JOB_HOME/runjob$i.inp
+	
+	#replace template placeholders with values
+	sed -i 's/ENERGYIN/'"${energy[i]}"'/g' $JOB_HOME/runjob$i.inp
+	#for -${energy[i]} the negative sign indicates energy, not momentum (use momentum for pions)
+	sed -i 's/TYPE/'"$TYPE"'/g' $JOB_HOME/runjob$i.inp  #beam particle type
+	sed -i 's/THICK/'"${THICKNESS[j]}"'/g' $JOB_HOME/runjob$i.inp  #thickness of polystyrene layer
+#	sed -i 's/ZED/'"$ZED"'/g' $JOB_HOME/runjob$i.inp  #atomic charge
+#	sed -i 's/ATN/'"$ATN"'/g' $JOB_HOME/runjob$i.inp  #atomic number
+
+	echo "Job: $j-$i   Ion type: $TYPE    Energy: ${energy[i]} GeV/u    Plane thickness: ${THICKNESS[j]} cm"
+	echo ""
+	
+	touch $JOB_HOME/fluka$j-$i.sh  #script to submit to the queue
+	echo "#!/bin/bash" >> $JOB_HOME/fluka$j-$i.sh
+	echo ". /local/env.sh" >> $JOB_HOME/fluka$j-$i.sh
+#	echo ". /cvmfs/lhcb.cern.ch/lib/lhcb/LBSCRIPTS/LBSCRIPTS_v9r2p2/InstallArea/scripts/LbLogin.sh -c x86_64-slc6-gcc7-opt" >> $JOB_HOME/fluka$j-$i.sh
+	echo "source /home/lhcb/leverington/setup.sh" >> $JOB_HOME/fluka$j-$i.sh
+        
+        #execute this command(s)
+	echo "cd "$JOB_HOME >> $JOB_HOME/fluka$j-$i.sh
+	echo "$FLUPRO/flutil/rfluka -e $FLUPRO/flutil/flukadpm3 -N0 -M1 runjob"$i >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "/work/leverington/readfluka/tools/eventdat2root runjob"$i"001_eventdata"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "$HOME/eventdatroot runjob"$i"001_eventdata "$j  >> $JOB_HOME/fluka$j-$i.sh
+	
+	qsub -l os=slc6 -l hio=10 -l ujl=40 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+
+#	$FLUPRO/flutil/rfluka -e /home/leverington/software/fluka/flutil/flukadpm3 -N0 -M1 $JOB_HOME/runjob$i  #run 
+#	sleep 2
+#	eventdat2root runjob$i\001_eventdata  #convert to root file
+#	sleep 2
+#	$HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+ 	sleep 1
+    done
+    sleep 3
+done
diff --git a/pion/runanalysis.sh b/pion/runanalysis.sh
new file mode 100755
index 0000000..7399985
--- /dev/null
+++ b/pion/runanalysis.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+HOME=/work/leverington/fluka_gfortran7/project/pion
+
+
+for j in {0..5}
+do
+    JOB_HOME=$HOME/jobs$j
+    cd $JOB_HOME
+
+	     for i in 0 #26
+	     do
+		 
+		 $HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+		 sleep 0.1 
+	     done
+done
diff --git a/pion/runjobtemplate.inp b/pion/runjobtemplate.inp
new file mode 100644
index 0000000..daa29c3
--- /dev/null
+++ b/pion/runjobtemplate.inp
@@ -0,0 +1,95 @@
+* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
+TITLE
+FLUKA Course Exercise
+*
+* use names everywhere and free format for geometry
+DEFAULTS                                                              HADROTHE
+*
+* beam definitions
+BEAM      ENERGYIN   0.0       -1.7      -0.8      -0.8               TYPE
+*HI-PROPE          6.       12.
+BEAMPOS                             -0.1
+DELTARAY       1E-06       0.0       0.0  HYDROGEN  @LASTMAT          NOPRINT
+PART-THR  -0.00001  PROTON    PROTON                 0.0   
+*
+* Geometry
+* --------
+GEOBEGIN                                                              COMBNAME
+    0    0             Cylindrical Target
+*
+* Bodies
+* ------
+*
+* Blackhole to include geometry
+SPH BLK        0.0 0.0 0.0 10000.
+* Void sphere
+SPH VOID       0.0 0.0 0.0 1000.
+* Infinite cylinder
+ZCC TARG       0.0 0.0 5.
+*  planes cutting the cylinder
+XYP ZTlow      0.0
+XYP T1seg      THICK
+XYP T2seg      2.
+XYP ZThigh     10.
+*  additional plane for scoring
+XYP ZTscor     9.9999
+END
+*
+* Regions
+* -------
+*
+* Blackhole
+BLKHOLE      5 +BLK -VOID
+*
+* Target segment 1
+TARGS1       5 +TARG -ZTlow +T1seg
+* Target segment 2
+TARGS2       5 +TARG -T1seg +T2seg
+* Target segment 3
+TARGS3       5 +TARG -T2seg +ZTscor
+* Scoring region
+TARGS4       5 +TARG -ZTscor +ZThigh
+*
+* Air around target
+INAIR        5 | +VOID -TARG
+               | +VOID +ZTlow
+               | +VOID -ZThigh
+END
+*
+GEOEND
+*
+* Materials definition
+* --------------------
+MATERIAL                        0.001965                              CO2
+COMPOUND          1.    CARBON        2.    OXYGEN                    CO2
+*
+* Assign materials
+* ----------------
+ASSIGNMA    BLCKHOLE   BLKHOLE
+ASSIGNMA    POLYSTYR    TARGS1
+ASSIGNMA      VACUUM    TARGS2
+ASSIGNMA      VACUUM    TARGS3
+ASSIGNMA      VACUUM    TARGS4
+ASSIGNMA      VACUUM     INAIR
+*
+* Exercise: heavy ions
+* --------------------
+*
+*  charge spectrum of ions
+USRYIELD       1422.  HEAVYION      -68.    TARGS1    TARGS2        1.FragZ1
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS2    TARGS3        1.FragZ2
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS3    TARGS4        1.FragZ3
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+*
+*  LET in water of ions and charged particles
+USRYIELD       2223.  HEAVYION      -69.    TARGS3    TARGS4        1.LETHI
+USRYIELD         20.       0.0      200.       9.5       2.5     2703. &
+USRYIELD       2223.  ALL-CHAR      -69.    TARGS3    TARGS4        1.LETCh
+USRYIELD         20.       0.0      200.       9.5      -2.5     2703. &
+EVENTDAT        -21.                                                  eventdata
+*
+RANDOMIZ          1.
+START         10000.                                     0.0
+STOP
diff --git a/proton/#langaus.C# b/proton/#langaus.C#
new file mode 100644
index 0000000..ce0132f
--- /dev/null
+++ b/proton/#langaus.C#
@@ -0,0 +1,352 @@
+ /// \ingroup tutorial_fit
+ /// \notebook
+ /// Convoluted Landau and Gaussian Fitting Function
+ ///         (using ROOT's Landau and Gauss functions)
+ ///
+ ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
+ ///
+ ///  to execute this example, do:
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C
+ /// ~~~
+ ///
+ /// or
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C++
+ /// ~~~
+ ///
+ /// \macro_image
+ /// \macro_output
+ /// \macro_code
+ ///
+ /// \authors H.Pernegger, Markus Friedl
+ 
+ #include "TH1.h"
+ #include "TF1.h"
+ #include "TROOT.h"
+ #include "TStyle.h"
+ #include "TMath.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+ 
+ Double_t langaufun(Double_t *x, Double_t *par) {
+ 
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+ 
+       // Numeric constants
+       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+       Double_t mpshift  = -0.22278298;       // Landau maximum location
+ 
+       // Control constants
+       Double_t np = 200.0;      // number of convolution steps
+       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
+ 
+       // Variables
+       Double_t xx;
+       Double_t mpc;
+       Double_t fland;
+       Double_t sum = 0.0;
+       Double_t xlow,xupp;
+       Double_t step;
+       Double_t i;
+ 
+ 
+       // MP shift correction
+       mpc = par[1] - mpshift * par[0];
+ 
+       // Range of convolution integral
+       xlow = x[0] - sc * par[3];
+       xupp = x[0] + 2*sc * par[3];
+ 
+       step = (xupp-xlow) / np;
+ 
+       // Convolution integral of Landau and Gaussian by sum
+       for(i=1.0; i<=np/2; i++) {
+          xx = xlow + (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+ 
+          xx = xupp - (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+       }
+ 
+       return (par[2] * step * sum * invsq2pi / par[3]);
+ }
+ 
+ 
+ 
+ TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
+ {
+    // Once again, here are the Landau * Gaussian parameters:
+    //   par[0]=Width (scale) parameter of Landau density
+    //   par[1]=Most Probable (MP, location) parameter of Landau density
+    //   par[2]=Total area (integral -inf to inf, normalization constant)
+    //   par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    // Variables for langaufit call:
+    //   his             histogram to fit
+    //   fitrange[2]     lo and hi boundaries of fit range
+    //   startvalues[4]  reasonable start values for the fit
+    //   parlimitslo[4]  lower parameter limits
+    //   parlimitshi[4]  upper parameter limits
+    //   fitparams[4]    returns the final fit parameters
+    //   fiterrors[4]    returns the final fit errors
+    //   ChiSqr          returns the chi square
+    //   NDF             returns ndf
+ 
+    Int_t i;
+    Char_t FunName[100];
+ 
+    sprintf(FunName,"Fitfcn_%s",his->GetName());
+ 
+    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
+    if (ffitold) delete ffitold;
+ 
+    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
+    ffit->SetParameters(startvalues);
+    ffit->SetParNames("Width","MP","Area","GSigma");
+ 
+    for (i=0; i<4; i++) {
+       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
+    }
+ 
+    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
+ 
+    ffit->GetParameters(fitparams);    // obtain fit parameters
+    for (i=0; i<4; i++) {
+       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
+    }
+    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
+    NDF[0] = ffit->GetNDF();           // obtain ndf
+ 
+    return (ffit);              // return fit function
+ 
+ }
+ 
+ 
+ Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
+ 
+    // Seaches for the location (x value) at the maximum of the
+    // Landau-Gaussian convolute and its full width at half-maximum.
+    //
+    // The search is probably not very efficient, but it's a first try.
+ 
+    Double_t p,x,fy,fxr,fxl;
+    Double_t step;
+    Double_t l,lold;
+    Int_t i = 0;
+    Int_t MAXCALLS = 10000;
+ 
+ 
+    // Search for maximum
+ 
+    p = params[1] - 0.1 * params[0];
+    step = 0.05 * params[0];
+    lold = -2.0;
+    l    = -1.0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = langaufun(&x,params);
+ 
+       if (l < lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-1);
+ 
+    maxx = x;
+ 
+    fy = l/2;
+ 
+ 
+    // Search for right x location of fy
+ 
+    p = maxx + params[0];
+    step = params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-2);
+ 
+    fxr = x;
+ 
+ 
+    // Search for left x location of fy
+ 
+    p = maxx - 0.5 * params[0];
+    step = -params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-3);
+ 
+ 
+    fxl = x;
+ 
+    FWHM = fxr - fxl;
+    return (0);
+ }
+ 
+void langaus() {
+    // Fill Histogram
+   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
+                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
+                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
+                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
+                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
+    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
+   */
+   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
+   TH1D * hSNR;
+
+    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
+    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+   
+    
+    // hSNR = h_beamSignal_b0[13];
+    // Fitting SNR histo
+    printf("Fitting...\n");
+    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
+    
+    // Setting fit range and start values
+    Double_t fr[2];
+    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
+    Double_t chisqr;
+    Int_t    ndf;
+    TF1 *fitsnr;
+    Double_t SNRPeak, SNRFWHM;
+    char rootfilename[50] = "";
+    char saveplotname[50] = "";
+    char fout_mpv_name[50] = "";
+
+    int j = 1;
+     Double_t norm;
+     ofstream myfile, fout_mpv;
+     myfile.open ("MPVcorrection_proton0mm05.txt");
+     sprintf(fout_mpv_name, "jobs%i/plots/mpv.txt",j);
+     fout_mpv.open (fout_mpv_name);
+
+    for (int i = 0; i<26;i++){
+
+      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile *rootFile = new TFile(rootfilename,"OPEN");
+
+      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
+      norm = hSNR->GetEntries();
+      hSNR->Scale(1/norm);
+      
+      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
+      fr[1]=hSNR->GetMean() + 2*hSNR->GetRMS();
+      //fr[1] = 2.0;
+      
+      pllo[0]=0.00001; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01; 
+      plhi[0]=0.1; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.1;
+      sv[0]=0.05;
+      //sv[1]=2.1;
+      sv[2]=0.01;
+      sv[3]=0.02;
+      sv[1] = hSNR->GetMean();
+    
+      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
+
+      
+      langaupro(fp,SNRPeak,SNRFWHM);
+      graph_x[i] = beta_proton[i];
+      graph_y[i] = fp[1];
+      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
+      fout_mpv << graph_x[i] << " " << graph_y[i] << " " << graph_yerr[i] << endl;
+      printf("Fitting done\nPlotting results...\n");
+      myfile << i << " " << graph_y[i] << endl;
+      // Global style settings
+      gStyle->SetOptStat(1111);
+      gStyle->SetOptFit(111);
+      //gStyle->SetLabelSize(0.03,"x");
+      //gStyle->SetLabelSize(0.03,"y");
+      
+      // hSNR->GetXaxis()->SetRange(0,70);
+      hSNR->Draw();
+      fitsnr->Draw("lsame");
+      c1->Update();
+      sprintf(saveplotname, "jobs%i/plots/runjob%i001_h_spratio.pdf",j,i);
+      c1->SaveAs(saveplotname);
+      //c1->WaitPrimitive();
+      hSNR->Delete();
+      rootFile->Close();
+    }
+     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+    graph_1->Draw("A*");
+    c2->SaveAs("MPVcorrection_proton0mm05.pdf");
+    myfile.close();
+    fout_mpv.close();
+ }
diff --git a/proton/MPVcorrection_proton0mm05 b/proton/MPVcorrection_proton0mm05
new file mode 100644
index 0000000..e3c67c0
--- /dev/null
+++ b/proton/MPVcorrection_proton0mm05
@@ -0,0 +1,117 @@
+%!PS-Adobe-2.0
+%%Title: MPVcorrection_proton0mm05: c2
+%%Creator: ROOT Version 6.12/06
+%%CreationDate: Fri Jul 19 11:40:02 2019
+%%Orientation: Landscape
+%%DocumentNeededResources: ProcSet (FontSetInit)
+%%EndComments
+%%BeginProlog
+/s {stroke} def /l {lineto} def /m {moveto} def /t {translate} def
+/r {rotate} def /rl {roll}  def /R {repeat} def
+/d {rlineto} def /rm {rmoveto} def /gr {grestore} def /f {eofill} def
+/c {setrgbcolor} def /black {0 setgray} def /sd {setdash} def
+/cl {closepath} def /sf {scalefont setfont} def /lw {setlinewidth} def
+/box {m dup 0 exch d exch 0 d 0 exch neg d cl} def
+/NC{systemdict begin initclip end}def/C{NC box clip newpath}def
+/bl {box s} def /bf {gsave box gsave f grestore 1 lw [] 0 sd s grestore} def /Y { 0 exch d} def /X { 0 d} def 
+/K {{pop pop 0 moveto} exch kshow} bind def
+/ita {/ang 15 def gsave [1 0 ang dup sin exch cos div 1 0 0] concat} def 
+/mp {newpath /y exch def /x exch def} def
+/side {[w .77 mul w .23 mul] .385 w mul sd w 0 l currentpoint t -144 r} def
+/mr {mp x y w2 0 360 arc} def /m24 {mr s} def /m20 {mr f} def
+/mb {mp x y w2 add m w2 neg 0 d 0 w neg d w 0 d 0 w d cl} def
+/mt {mp x y w2 add m w2 neg w neg d w 0 d cl} def
+/w4 {w 4 div} def
+/w6 {w 6 div} def
+/w8 {w 8 div} def
+/m21 {mb f} def /m25 {mb s} def /m22 {mt f} def /m26{mt s} def
+/m23 {mp x y w2 sub m w2 w d w neg 0 d cl f} def
+/m27 {mp x y w2 add m w3 neg w2 neg d w3 w2 neg d w3 w2 d cl s} def
+/m28 {mp x w2 sub y w2 sub w3 add m w3 0 d  0 w3 neg d w3 0 d 0 w3 d w3 0 d  0 w3 d w3 neg 0 d 0 w3 d w3 neg 0 d 0 w3 neg d w3 neg 0 d cl s } def
+/m29 {mp gsave x w2 sub y w2 add w3 sub m currentpoint t 4 {side} repeat cl fill gr} def
+/m30 {mp gsave x w2 sub y w2 add w3 sub m currentpoint t 4 {side} repeat cl s gr} def
+/m31 {mp x y w2 sub m 0 w d x w2 sub y m w 0 d x w2 sub y w2 add m w w neg d x w2 sub y w2 sub m w w d s} def
+/m32 {mp x y w2 sub m w2 w d w neg 0 d cl s} def
+/m33 {mp x y w2 add m w3 neg w2 neg d w3 w2 neg d w3 w2 d cl f} def
+/m34 {mp x w2 sub y w2 sub w3 add m w3 0 d  0 w3 neg d w3 0 d 0 w3 d w3 0 d  0 w3 d w3 neg 0 d 0 w3 d w3 neg 0 d 0 w3 neg d w3 neg 0 d cl f } def
+/m35 {mp x y w2 add m w2 neg w2 neg d w2 w2 neg d w2 w2 d w2 neg w2 d x y w2 sub m 0 w d x w2 sub y m w 0 d s} def
+/m36 {mb x w2 sub y w2 add m w w neg d x w2 sub y w2 sub m w w d s} def
+/m37 {mp x y m w4 neg w2 d w4 neg w2 neg d w2 0 d  w4 neg w2 neg d w2 0 d w4 neg w2 d w2 0 d w4 neg w2 d w4 neg w2 neg d cl s} def
+/m38 {mp x w4 sub y w2 add m w4 neg w4 neg d 0 w2 neg d w4 w4 neg d w2 0 d w4 w4 d 0 w2 d w4 neg w4 d w2 neg 0 d x y w2 sub m 0 w d x w2 sub y m w 0 d cl s} def
+/m39 {mp x y m w4 neg w2 d w4 neg w2 neg d w2 0 d  w4 neg w2 neg d w2 0 d w4 neg w2 d w2 0 d w4 neg w2 d w4 neg w2 neg d cl f} def
+/m40 {mp x y m w4 w2 d w4 w4 neg d w2 neg w4 neg d w2 w4 neg d w4 neg w4 neg d w4 neg w2 d w4 neg w2 neg d w4 neg w4 d w2 w4 d w2 neg w4 d w4 w4 d w4 w2 neg d cl s} def
+/m41 {mp x y m w4 w2 d w4 w4 neg d w2 neg w4 neg d w2 w4 neg d w4 neg w4 neg d w4 neg w2 d w4 neg w2 neg d w4 neg w4 d w2 w4 d w2 neg w4 d w4 w4 d w4 w2 neg d cl f} def
+/m42 {mp x y w2 add m w8 neg w2 -3 4 div mul d w2 -3 4 div mul w8 neg d w2 3 4 div mul w8 neg d w8 w2 -3 4 div mul d w8 w2 3 4 div mul d w2 3 4 div mul w8 d w2 -3 4 div mul w8 d w8 neg w2 3 4 div mul d cl s} def
+/m43 {mp x y w2 add m w8 neg w2 -3 4 div mul d w2 -3 4 div mul w8 neg d w2 3 4 div mul w8 neg d w8 w2 -3 4 div mul d w8 w2 3 4 div mul d w2 3 4 div mul w8 d w2 -3 4 div mul w8 d w8 neg w2 3 4 div mul d cl f} def
+/m44 {mp x y m w6 neg w2 d w2 2 3 div mul 0 d w6 neg w2 neg d w2 w6 d 0 w2 -2 3 div mul d w2 neg w6 d w6 w2 neg d w2 -2 3 div mul 0 d w6 w2 d w2 neg w6 neg d 0 w2 2 3 div mul d w2 w6 neg d cl s} def
+/m45 {mp x y m w6 neg w2 d w2 2 3 div mul 0 d w6 neg w2 neg d w2 w6 d 0 w2 -2 3 div mul d w2 neg w6 d w6 w2 neg d w2 -2 3 div mul 0 d w6 w2 d w2 neg w6 neg d 0 w2 2 3 div mul d w2 w6 neg d cl f} def
+/m46 {mp x y w4 add m w4 neg w4 d w4 neg w4 neg d  w4 w4 neg d w4 neg w4 neg d w4 w4 neg d w4 w4 d w4 w4 neg d w4 w4 d w4 neg w4 d w4 w4 d w4 neg w4 d w4 neg w4 neg d cl s} def
+/m47 {mp x y w4 add m w4 neg w4 d w4 neg w4 neg d w4 w4 neg d w4 neg w4 neg d  w4 w4 neg d w4 w4 d w4 w4 neg d w4 w4 d w4 neg w4 d w4 w4 d w4 neg w4 d w4 neg w4 neg d cl f} def
+/m48 {mp x y w4 add m w4 neg w4 d w4 neg w4 neg d w4 w4 neg d  w4 neg w4 neg d w4 w4 neg d w4 w4 d w4 w4 neg d w4 w4 d w4 neg w4 d w4 w4 d w4 neg w4 d w4 neg w4 neg d  w4 w4 neg d w4 neg w4 neg d w4 neg w4 d w4 w4 d cl f} def
+/m49 {mp x w2 sub w3 add y w2 sub w3 add m  0 w3 neg d w3 0 d 0 w3 d w3 0 d 0 w3 d w3 neg 0 d 0 w3 d w3 neg 0 d 0 w3 neg d w3 neg 0 d 0 w3 neg d w3 0 d 0 w3 d w3 0 d 0 w3 neg d w3 neg 0 d cl f } def
+/m2 {mp x y w2 sub m 0 w d x w2 sub y m w 0 d s} def
+/m5 {mp x w2 sub y w2 sub m w w d x w2 sub y w2 add m w w neg d s} def
+%%IncludeResource: ProcSet (FontSetInit)
+%%IncludeResource: font Times-Roman
+%%IncludeResource: font Times-Italic
+%%IncludeResource: font Times-Bold
+%%IncludeResource: font Times-BoldItalic
+%%IncludeResource: font Helvetica
+%%IncludeResource: font Helvetica-Oblique
+%%IncludeResource: font Helvetica-Bold
+%%IncludeResource: font Helvetica-BoldOblique
+%%IncludeResource: font Courier
+%%IncludeResource: font Courier-Oblique
+%%IncludeResource: font Courier-Bold
+%%IncludeResource: font Courier-BoldOblique
+%%IncludeResource: font Symbol
+%%IncludeResource: font ZapfDingbats
+/reEncode {exch findfont dup length dict begin {1 index /FID eq  {pop pop} {def} ifelse } forall /Encoding exch def currentdict end dup /FontName get exch definefont pop } def [/Times-Bold /Times-Italic /Times-BoldItalic /Helvetica /Helvetica-Oblique
+ /Helvetica-Bold /Helvetica-BoldOblique /Courier /Courier-Oblique /Courier-Bold /Courier-BoldOblique /Times-Roman /AvantGarde-Book /AvantGarde-BookOblique /AvantGarde-Demi /AvantGarde-DemiOblique /Bookman-Demi /Bookman-DemiItalic /Bookman-Light
+ /Bookman-LightItalic /Helvetica-Narrow /Helvetica-Narrow-Bold /Helvetica-Narrow-BoldOblique /Helvetica-Narrow-Oblique /NewCenturySchlbk-Roman /NewCenturySchlbk-Bold /NewCenturySchlbk-BoldItalic /NewCenturySchlbk-Italic /Palatino-Bold
+ /Palatino-BoldItalic /Palatino-Italic /Palatino-Roman ] {ISOLatin1Encoding reEncode } forall/Zone {/iy exch def /ix exch def  ix 1 sub  3144 mul  1 iy sub  2224 mul t} def
+%%EndProlog
+%%BeginSetup
+%%EndSetup
+newpath  gsave  90 r 0 -594 t  28 20 t .25 .25 scale  gsave 
+%%Page: 1 1
+ gsave  gsave 
+ 1 1 Zone
+ gsave  0 0 t black[  ] 0 sd 3 lw 1 1 1 c 2948 2117 0 0 bf black 1 1 1 c 2358 1693 295 212 bf black 2358 1693 295 212 bl 295 212 m 2358 X s 426 262 m -50 Y s 503 237 m -25 Y s 579 237 m -25 Y s 656 237 m -25 Y s 733 237 m -25 Y s 809 262 m -50 Y s
+ 886 237 m -25 Y s 963 237 m -25 Y s 1039 237 m -25 Y s 1116 237 m -25 Y s 1193 262 m -50 Y s 1269 237 m -25 Y s 1346 237 m -25 Y s 1423 237 m -25 Y s 1499 237 m -25 Y s 1576 262 m -50 Y s 1653 237 m -25 Y s 1729 237 m -25 Y s 1806 237 m -25 Y s 1883
+ 237 m -25 Y s 1959 262 m -50 Y s 2036 237 m -25 Y s 2113 237 m -25 Y s 2189 237 m -25 Y s 2266 237 m -25 Y s 2343 262 m -50 Y s 426 262 m -50 Y s 349 237 m -25 Y s 2343 262 m -50 Y s 2420 237 m -25 Y s 2496 237 m -25 Y s 2573 237 m -25 Y s 2650 237
+ m -25 Y s
+ gsave  2948 2117 0 0 C 376.816 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.3) show NC gr 
+ gsave  2948 2117 0 0 C 738.855 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.35) show NC gr 
+ gsave  2948 2117 0 0 C 1141.53 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.4) show NC gr 
+ gsave  2948 2117 0 0 C 1503.57 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.45) show NC gr 
+ gsave  2948 2117 0 0 C 1909.94 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.5) show NC gr 
+ gsave  2948 2117 0 0 C 2271.98 144.077 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.55) show NC gr  295 212 m 1693 Y s 366 272 m -71 X s 330 313 m -35 X s 330 354 m -35 X s 330 395 m -35 X s 366 436 m -71 X s 330 477 m -35 X s 330 518 m -35 X s 330
+ 559 m -35 X s 366 600 m -71 X s 330 642 m -35 X s 330 683 m -35 X s 330 724 m -35 X s 366 765 m -71 X s 330 806 m -35 X s 330 847 m -35 X s 330 888 m -35 X s 366 929 m -71 X s 330 970 m -35 X s 330 1011 m -35 X s 330 1052 m -35 X s 366 1093 m -71 X
+ s 330 1134 m -35 X s 330 1175 m -35 X s 330 1217 m -35 X s 366 1258 m -71 X s 330 1299 m -35 X s 330 1340 m -35 X s 330 1381 m -35 X s 366 1422 m -71 X s 330 1463 m -35 X s 330 1504 m -35 X s 330 1545 m -35 X s 366 1586 m -71 X s 330 1627 m -35 X s
+ 330 1668 m -35 X s 330 1709 m -35 X s 366 1750 m -71 X s 366 272 m -71 X s 330 231 m -35 X s 366 1750 m -71 X s 330 1792 m -35 X s 330 1833 m -35 X s 330 1874 m -35 X s
+ gsave  2948 2117 0 0 C 181.019 247.516 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.8) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 413.759 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.82) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 576.307 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.84) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 742.549 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.86) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 905.097 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.88) show NC gr 
+ gsave  2948 2117 0 0 C 181.019 1067.65 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.9) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1233.89 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.92) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1396.44 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.94) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1562.68 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.96) show NC gr 
+ gsave  2948 2117 0 0 C 140.382 1725.23 t 0 r /Helvetica findfont 70.1912 sf 0 0 m (0.98) show NC gr  1 1 1 c black
+ gsave  2948 2117 0 0 C 1333.63 2013.38 t 0 r /Helvetica findfont 103.44 sf 0 0 m (Graph) show NC gr  1 1 1 c black /w 30 def /w2 {w 2 div} def /w3 {w 3 div} def 491 1667 740 1435 933 1191 1091 1078 1225 1003 1342 933 1445 892 1539 830 1625 771 1703
+ 736 1777 705 1845 685 1908 660 1969 628 2026 607 2080 581 2132 545 2183 557 2231 533 2277 499 2340 490 2400 468 2457 434 23 { m31} R 491 1667 m 97 Y s 484 1764 m 15 X s 491 1667 m -97 Y s 484 1570 m 15 X s 740 1435 m 95 Y s 733 1530 m 15 X s 740
+ 1435 m -94 Y s 733 1341 m 15 X s 933 1191 m 91 Y s 926 1282 m 15 X s 933 1191 m -91 Y s 926 1100 m 15 X s 1091 1078 m 89 Y s 1084 1167 m 14 X s 1091 1078 m -89 Y s 1084 989 m 14 X s 1225 1003 m 89 Y s 1218 1092 m 14 X s 1225 1003 m -89 Y s 1218 914
+ m 14 X s 1342 933 m 87 Y s 1334 1020 m 15 X s 1342 933 m -88 Y s 1334 845 m 15 X s 1445 892 m 87 Y s 1438 979 m 15 X s 1445 892 m -87 Y s 1438 805 m 15 X s 1539 830 m 87 Y s 1532 917 m 15 X s 1539 830 m -86 Y s 1532 744 m 15 X s 1625 771 m 86 Y s
+ 1617 857 m 15 X s 1625 771 m -85 Y s 1617 686 m 15 X s 1703 736 m 85 Y s 1696 821 m 15 X s 1703 736 m -85 Y s 1696 651 m 15 X s 1777 705 m 84 Y s 1769 789 m 15 X s 1777 705 m -85 Y s 1769 620 m 15 X s 1845 685 m 84 Y s 1837 769 m 15 X s 1845 685 m
+ -85 Y s 1837 600 m 15 X s 1908 660 m 84 Y s 1901 744 m 15 X s 1908 660 m -85 Y s 1901 575 m 15 X s 1969 628 m 84 Y s 1961 712 m 15 X s 1969 628 m -83 Y s 1961 545 m 15 X s 2026 607 m 84 Y s 2018 691 m 15 X s 2026 607 m -83 Y s 2018 524 m 15 X s 2080
+ 581 m 84 Y s 2073 665 m 15 X s 2080 581 m -83 Y s 2073 498 m 15 X s 2132 545 m 83 Y s 2125 628 m 15 X s 2132 545 m -83 Y s 2125 462 m 15 X s 2183 557 m 82 Y s 2175 639 m 15 X s 2183 557 m -83 Y s 2175 474 m 15 X s 2231 533 m 83 Y s 2223 616 m 15 X s
+ 2231 533 m -82 Y s 2223 451 m 15 X s 2277 499 m 83 Y s 2269 582 m 15 X s 2277 499 m -82 Y s 2269 417 m 15 X s 2340 490 m 82 Y s 2332 572 m 15 X s 2340 490 m -82 Y s 2332 408 m 15 X s 2400 468 m 82 Y s 2392 550 m 15 X s 2400 468 m -82 Y s 2392 386 m
+ 15 X s 2457 434 m 82 Y s 2449 516 m 15 X s 2457 434 m -81 Y s 2449 353 m 15 X s 1 1 1 c black
+ gsave  2948 2117 0 0 C 1333.63 2013.38 t 0 r /Helvetica findfont 103.44 sf 0 0 m (Graph) show NC gr  gr showpage
+ gr 
+%%Trailer
+%%Pages:  1
+ gr  gr  gr 
+%%EOF
diff --git a/proton/MPVcorrection_proton0mm05.pdf b/proton/MPVcorrection_proton0mm05.pdf
new file mode 100644
index 0000000..7ce993a
Binary files /dev/null and b/proton/MPVcorrection_proton0mm05.pdf differ
diff --git a/proton/MPVcorrection_proton0mm05.txt b/proton/MPVcorrection_proton0mm05.txt
new file mode 100644
index 0000000..bf08015
--- /dev/null
+++ b/proton/MPVcorrection_proton0mm05.txt
@@ -0,0 +1,26 @@
+0 0.939261
+1 0.923705
+2 0.908651
+3 0.897928
+4 0.888853
+5 0.881592
+6 0.876172
+7 0.869721
+8 0.862654
+9 0.858657
+10 0.854603
+11 0.851932
+12 0.849187
+13 0.845096
+14 0.843032
+15 0.839206
+16 0.834533
+17 0.835875
+18 0.833985
+19 0.82905
+20 0.828067
+21 0.825095
+22 0.821344
+23 0.819369
+24 0.817124
+25 0.816329
diff --git a/proton/Makefile b/proton/Makefile
new file mode 100644
index 0000000..c16461b
--- /dev/null
+++ b/proton/Makefile
@@ -0,0 +1,38 @@
+
+ROOTCFLAGS   := $(shell root-config --cflags)
+ROOTLIBS     := $(shell root-config --libs)
+ROOTGLIBS    := $(shell root-config --glibs)
+
+GSLCFLAGS   := $(shell gsl-config --cflags)
+GSLLIBS     := $(shell gsl-config --libs)
+GSLGLIBS    := $(shell gsl-config --glibs)
+
+
+LDFLAGS       = -O
+LIBS         += $(ROOTLIBS) $(GSLLIBS) -L/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/lib
+CFLAGS       += $(ROOTCFLAGS) $(GSLCFLAGS)  -I/cvmfs/lhcb.cern.ch/lib/lcg/releases/tbb/2018_U1-d3621/x86_64-slc6-gcc7-opt/include
+
+#  Not sure why Minuit isn't being included -- put in by hand
+#
+LIBS         += -lMinuit -ltbb
+
+
+all: eventdatroot
+
+
+
+eventdatroot: eventdatroot.o
+	g++ -o eventdatroot  eventdatroot.o $(LDFLAGS) $(LIBS)
+
+%.o: %.c
+	g++ ${CFLAGS} -c  -g $<
+
+
+test:
+	@echo $(ROOTCFLAGS) 
+	@echo $(LDFLAGS) 
+	@echo $(LIBS)
+
+clean:	
+	-rm -f  *.o
+
diff --git a/proton/eventdatroot b/proton/eventdatroot
new file mode 100755
index 0000000..2473e30
Binary files /dev/null and b/proton/eventdatroot differ
diff --git a/proton/eventdatroot.c b/proton/eventdatroot.c
new file mode 100644
index 0000000..a13394e
--- /dev/null
+++ b/proton/eventdatroot.c
@@ -0,0 +1,135 @@
+#define eventdatroot_cxx
+#include "eventdatroot.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TH1.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+//#include <TROOT.h>
+#include <TF1.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+
+
+using namespace std;
+
+void eventdatroot::Loop(Float_t scale)
+{
+  //   In a ROOT session, you can do:
+  //      root> .L eventdatroot.C
+  //      root> eventdatroot t
+  //      root> t.GetEntry(12); // Fill t data members with entry number 12
+  //      root> t.Show();       // Show values of entry 12
+  //      root> t.Show(16);     // Read and show values of entry 16
+  //      root> t.Loop();       // Loop on all entries
+  //
+
+  //     This is the loop skeleton where:
+  //    jentry is the global entry number in the chain
+  //    ientry is the entry number in the current Tree
+  //  Note that the argument to GetEntry must be:
+  //    jentry for TChain::GetEntry
+  //    ientry for TTree::GetEntry and TBranch::GetEntry
+  //
+  //       To read only selected branches, Insert statements like:
+  // METHOD1:
+  //    fChain->SetBranchStatus("*",0);  // disable all branches
+  //    fChain->SetBranchStatus("branchname",1);  // activate branchname
+  // METHOD2: replace line
+  //    fChain->GetEntry(jentry);       //read all branches
+  //by  b_branchname->GetEntry(ientry); //read only this branch
+
+  Float_t ratio_electron, ratio_primary, etotal, spratio;
+  TF1 * f_sp_ps = new TF1("f_sp_ps","[0]*pow(x,[1])+[2]", 50, 250); //stopping power of protons in polystyrene [MeV cm2 /g]
+  f_sp_ps->SetParameters(361.936, -0.892255, 1.19133);  //protons between 50 and 220 MeV
+
+  
+  if (fChain == 0) return;
+
+  Long64_t nentries = fChain->GetEntriesFast();
+
+  Long64_t nbytes = 0, nb = 0;
+  for (Long64_t jentry=0; jentry<nentries;jentry++) {
+    Long64_t ientry = LoadTree(jentry);
+    if (ientry < 0) break;
+    nb = fChain->GetEntry(jentry);   nbytes += nb;
+    // if (Cut(ientry) < 0) continue;
+    //    cout << jentry << endl;
+
+    h_endist_0->Fill(ENDIST[0]); // edep by ionisation
+    h_endist_1->Fill(ENDIST[1]); // edep by pi0, e-, e+ and #gamma
+    h_endist_2->Fill(ENDIST[2]); // edep by nuclear recoils and heavy fragments
+    h_endist_3->Fill(ENDIST[3]); // edep by part. #lt threshold
+    h_endist_4->Fill(ENDIST[4]); // E leaving the syste
+    h_endist_5->Fill(ENDIST[5]); // E carried by discarded particles
+    h_endist_6->Fill(ENDIST[6]); // resid, excit. E after evap.
+    h_endist_7->Fill(ENDIST[7]); // edep by low-energy neutrons
+    h_endist_8->Fill(ENDIST[8]); // E of part. #gt time limit
+    h_endist_9->Fill(ENDIST[9]); // E lost in endothermic nuclear reactions
+    h_endist_10->Fill(ENDIST[10]); // E lost in endothermic low-E n-reaction
+    h_endist_11->Fill(ENDIST[11]); // missing E
+
+    etotal = 0.;
+    for (unsigned int i = 0; i<3; i++ ){
+      etotal += ENDIST[i];
+    }
+
+    ratio_electron = ENDIST[1] / etotal;
+    ratio_primary = ENDIST[0] / etotal;
+
+    h_ratio_e->Fill(ratio_electron);
+    h_ratio_p->Fill(ratio_primary);
+    h_sum_ep->Fill(ENDIST[1]+ENDIST[0]);
+
+    h_let_ep->Fill( (ENDIST[1]+ENDIST[0])/scale );
+    h_let_e->Fill( (ENDIST[1])/scale );
+    h_let_p->Fill( (ENDIST[0])/scale );
+
+    spratio = ( ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06) / f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) ); //ratio of deposited energy to PSTAR stopping power (MeV/cm/density) / MeVcm2/g
+    h_spratio->Fill(spratio);
+
+    //   cout <<  ( (ENDIST[0]+ENDIST[1])/scale*1000/1.06)  << " " << f_sp_ps->Eval( (DATA_ENERGY[0]+ENDIST[4])*1000 ) << " " << spratio << endl;
+
+  }
+}
+
+
+
+int main(int argc, const char* argv[])
+{
+  Int_t thickkey = atoi(argv[2]);
+  Float_t THICKNESS[6] = {0.025, 0.050, 0.100, 0.150, 0.500, 1.000}; //#cm
+  // Float_t THICKNESS[5] = {0.050, 0.100, 0.150, 1.000, 10.00}; //#cm
+
+  ///open the input and output root files
+  char  finname[50];
+  sprintf(finname, "%s.root",argv[1]);
+  TFile * f = new TFile(finname,"OPEN");
+
+  char  foutname[50];
+  sprintf(foutname, "%s_out.root",argv[1]);
+  cout << foutname << endl;
+  TFile * fout = new TFile(foutname, "RECREATE");
+
+  //run the analysis loop above
+  eventdatroot t;
+  TTree * tree;
+  f->GetObject("EVENTDAT",tree);  
+  t.Init(tree);
+  t.Loop( THICKNESS[thickkey] );
+  t.Closefile(fout);
+
+  return 1;
+
+}
diff --git a/proton/eventdatroot.h b/proton/eventdatroot.h
new file mode 100644
index 0000000..0026047
--- /dev/null
+++ b/proton/eventdatroot.h
@@ -0,0 +1,199 @@
+//////////////////////////////////////////////////////////
+// This class has been automatically generated on
+// Fri Mar  8 13:05:30 2019 by ROOT version 6.12/06
+// from TTree EVENTDAT/FLUKA Course Exercise	DATE:  3/ 4/19,  TIME: 11:54:46
+// found on file: jobs/runjob1001_eventdata.root
+//////////////////////////////////////////////////////////
+
+#ifndef eventdatroot_h
+#define eventdatroot_h
+
+#include <TROOT.h>
+#include <TChain.h>
+#include <TFile.h>
+#include <TH1.h>
+
+// Header file for the classes stored in the TTree if any.
+
+class eventdatroot {
+public :
+   TTree          *fChain;   //!pointer to the analyzed TTree or TChain
+   Int_t           fCurrent; //!current Tree number in a TChain
+
+// Fixed size dimensions of array or collections stored in the TTree if any.
+
+   // Declaration of leaf types
+   Float_t         DATA_ALL_PART[6];
+   Float_t         DATA_BEAMPART[6];
+   Float_t         DATA_ENERGY[6];
+   Float_t         DATA_id211[6];
+   Int_t           SEEDS[4];
+   Float_t         ENDIST[12];
+
+   // List of branches
+   TBranch        *b_DATA;   //!
+   TBranch        *b_seed;   //!
+   TBranch        *b_endist;   //!
+
+      
+   TH1F * h_endist_0;
+   TH1F * h_endist_1;
+   TH1F * h_endist_2;
+   TH1F * h_endist_3;
+   TH1F * h_endist_4;
+   TH1F * h_endist_5;
+   TH1F * h_endist_6;
+   TH1F * h_endist_7;
+   TH1F * h_endist_8;
+   TH1F * h_endist_9;
+   TH1F * h_endist_10;
+   TH1F * h_endist_11;
+
+   TH1F * h_ratio_e;
+   TH1F * h_ratio_p;
+   TH1F * h_sum_ep;
+   
+   TH1F * h_let_ep;
+   TH1F * h_let_e;
+   TH1F * h_let_p;
+
+   TH1F * h_spratio;
+
+   eventdatroot(TTree *tree=0);
+   virtual ~eventdatroot();
+   virtual Int_t    Cut(Long64_t entry);
+   virtual Int_t    GetEntry(Long64_t entry);
+   virtual Long64_t LoadTree(Long64_t entry);
+   virtual void     Init(TTree *tree);
+   virtual void     Loop(Float_t scale);
+   virtual Bool_t   Notify();
+   virtual void     Show(Long64_t entry = -1);
+   virtual Int_t    Closefile(TFile * file);
+};
+
+#endif
+
+#ifdef eventdatroot_cxx
+eventdatroot::eventdatroot(TTree *tree) : fChain(0) 
+{
+// if parameter tree is not specified (or zero), connect the file
+// used to generate this class and read the Tree.
+  //  if (tree == 0) {
+     //   TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject("jobs/runjob1001_eventdata.root");
+      //     if (!f || !f->IsOpen()) {
+		//      f = new TFile("jobs/runjob1001_eventdata.root");
+	 //     }
+      //    f->GetObject("EVENTDAT",tree);
+
+      //   }
+   //   Init(tree);
+}
+
+eventdatroot::~eventdatroot()
+{
+   if (!fChain) return;
+   delete fChain->GetCurrentFile();
+}
+
+Int_t eventdatroot::GetEntry(Long64_t entry)
+{
+// Read contents of entry.
+   if (!fChain) return 0;
+   return fChain->GetEntry(entry);
+}
+Long64_t eventdatroot::LoadTree(Long64_t entry)
+{
+// Set the environment to read one entry
+   if (!fChain) return -5;
+   Long64_t centry = fChain->LoadTree(entry);
+   if (centry < 0) return centry;
+   if (fChain->GetTreeNumber() != fCurrent) {
+      fCurrent = fChain->GetTreeNumber();
+      Notify();
+   }
+   return centry;
+}
+
+void eventdatroot::Init(TTree *tree)
+{
+   // The Init() function is called when the selector needs to initialize
+   // a new tree or chain. Typically here the branch addresses and branch
+   // pointers of the tree will be set.
+   // It is normally not necessary to make changes to the generated
+   // code, but the routine can be extended by the user if needed.
+   // Init() will be called many times when running on PROOF
+   // (once per file to be processed).
+
+   // Set branch addresses and branch pointers
+   if (!tree) return;
+   fChain = tree;
+   fCurrent = -1;
+   fChain->SetMakeClass(1);
+
+   fChain->SetBranchAddress("DATA", DATA_ALL_PART, &b_DATA);
+   fChain->SetBranchAddress("SEEDS", SEEDS, &b_seed);
+   fChain->SetBranchAddress("ENDIST", ENDIST, &b_endist);
+
+   h_endist_0 = new TH1F("h_endist_0","edep by ionisation",200,0.0,0.004);
+   h_endist_1 = new TH1F("h_endist_1","edep by pi0, e-, e+ and #gamma",200,0.0,0.005);
+   h_endist_2 = new TH1F("h_endist_2","edep by nuclear recoils and heavy fragments",200,0.0,0.001);
+   h_endist_3 = new TH1F("h_endist_3","edep by part. #lt threshold",200,0.0,0.001);
+   h_endist_4 = new TH1F("h_endist_4","E leaving the system",200,0.0,0.1);
+   h_endist_5 = new TH1F("h_endist_5","E carried by discarded particles",200,0.0,0.1);
+   h_endist_6 = new TH1F("h_endist_6","resid, excit. E after evap.",200,0.0,0.1);
+   h_endist_7 = new TH1F("h_endist_7","edep by low-energy neutrons",200,0.0,0.1);
+   h_endist_8 = new TH1F("h_endist_8","E of part. #gt time limit",200,0.0,0.1);
+   h_endist_9 = new TH1F("h_endist_9","E lost in endothermic nuclear reactions",200,0.0,0.1);
+   h_endist_10 = new TH1F("h_endist_10","E lost in endothermic low-E n-reactions",200,0.0,0.1);
+   h_endist_11 = new TH1F("h_endist_11","missing E",200,0.0,0.01);
+
+   h_ratio_e = new TH1F("h_ratio_e","fraction of edep by pi0, e-, e+ and #gamma  ",200,0.0,1.0);
+   h_ratio_p = new TH1F("h_ratio_p","fraction of edep by ionisation",200,0.0,1.0);
+   h_sum_ep = new TH1F("h_sum_ep","sum of edep by ion, pi0, e-, e+ and #gamma,",200,0.0,0.003);
+
+   h_let_ep = new TH1F("h_let_ep","LET/Z^{2} (GeV/cm) by ion, pi0, e-, e+ and #gamma,",200,0.0,0.05);
+   h_let_e = new TH1F("h_let_e","LET/Z^{2} (GeV/cm) by pi0, e-, e+ and #gamma,",200,0.0,0.02);
+   h_let_p = new TH1F("h_let_p","LET/Z^{2} (GeV/cm) by ion",200,0.0,0.02);
+
+   h_spratio = new TH1F("h_spratio","total LET/SP (MeV/cm/#rho)/SP(MeVcm2/g)",200,0.5,3.0);
+
+   
+   Notify();
+}
+
+Bool_t eventdatroot::Notify()
+{
+   // The Notify() function is called when a new file is opened. This
+   // can be either for a new TTree in a TChain or when when a new TTree
+   // is started when using PROOF. It is normally not necessary to make changes
+   // to the generated code, but the routine can be extended by the
+   // user if needed. The return value is currently not used.
+
+   return kTRUE;
+}
+
+void eventdatroot::Show(Long64_t entry)
+{
+// Print contents of entry.
+// If entry is not specified, print current entry
+   if (!fChain) return;
+   fChain->Show(entry);
+}
+Int_t eventdatroot::Cut(Long64_t entry)
+{
+// This function may be called from Loop.
+// returns  1 if entry is accepted.
+// returns -1 otherwise.
+   return 1;
+}
+
+Int_t eventdatroot::Closefile(TFile * file)
+{
+  file->cd();
+  file->Write();
+  file->Close();
+
+  return 1;
+}
+
+#endif // #ifdef eventdatroot_cxx
diff --git a/proton/langaus.C b/proton/langaus.C
new file mode 100644
index 0000000..ce0132f
--- /dev/null
+++ b/proton/langaus.C
@@ -0,0 +1,352 @@
+ /// \ingroup tutorial_fit
+ /// \notebook
+ /// Convoluted Landau and Gaussian Fitting Function
+ ///         (using ROOT's Landau and Gauss functions)
+ ///
+ ///  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
+ ///
+ ///  to execute this example, do:
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C
+ /// ~~~
+ ///
+ /// or
+ ///
+ /// ~~~{.cpp}
+ ///  root > .x langaus.C++
+ /// ~~~
+ ///
+ /// \macro_image
+ /// \macro_output
+ /// \macro_code
+ ///
+ /// \authors H.Pernegger, Markus Friedl
+ 
+ #include "TH1.h"
+ #include "TF1.h"
+ #include "TROOT.h"
+ #include "TStyle.h"
+ #include "TMath.h"
+#include <TH2.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <iostream>
+#include <fstream>
+#include <iomanip> 
+#include <cstdlib>
+#include <cmath>
+#include <string.h>
+#include <TLorentzVector.h>
+#include <vector>
+#include <TApplication.h>
+#include <TAxis.h>
+#include <stdio.h>
+#include <stdlib.h>
+#include <TFile.h>
+#include <TTree.h>
+ 
+ Double_t langaufun(Double_t *x, Double_t *par) {
+ 
+    //Fit parameters:
+    //par[0]=Width (scale) parameter of Landau density
+    //par[1]=Most Probable (MP, location) parameter of Landau density
+    //par[2]=Total area (integral -inf to inf, normalization constant)
+    //par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    //In the Landau distribution (represented by the CERNLIB approximation),
+    //the maximum is located at x=-0.22278298 with the location parameter=0.
+    //This shift is corrected within this function, so that the actual
+    //maximum is identical to the MP parameter.
+ 
+       // Numeric constants
+       Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
+       Double_t mpshift  = -0.22278298;       // Landau maximum location
+ 
+       // Control constants
+       Double_t np = 200.0;      // number of convolution steps
+       Double_t sc =   4.0;      // convolution extends to +-sc Gaussian sigmas
+ 
+       // Variables
+       Double_t xx;
+       Double_t mpc;
+       Double_t fland;
+       Double_t sum = 0.0;
+       Double_t xlow,xupp;
+       Double_t step;
+       Double_t i;
+ 
+ 
+       // MP shift correction
+       mpc = par[1] - mpshift * par[0];
+ 
+       // Range of convolution integral
+       xlow = x[0] - sc * par[3];
+       xupp = x[0] + 2*sc * par[3];
+ 
+       step = (xupp-xlow) / np;
+ 
+       // Convolution integral of Landau and Gaussian by sum
+       for(i=1.0; i<=np/2; i++) {
+          xx = xlow + (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+ 
+          xx = xupp - (i-.5) * step;
+          fland = TMath::Landau(xx,mpc,par[0]) / par[0];
+          sum += fland * TMath::Gaus(x[0],xx,par[3]);
+       }
+ 
+       return (par[2] * step * sum * invsq2pi / par[3]);
+ }
+ 
+ 
+ 
+ TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
+ {
+    // Once again, here are the Landau * Gaussian parameters:
+    //   par[0]=Width (scale) parameter of Landau density
+    //   par[1]=Most Probable (MP, location) parameter of Landau density
+    //   par[2]=Total area (integral -inf to inf, normalization constant)
+    //   par[3]=Width (sigma) of convoluted Gaussian function
+    //
+    // Variables for langaufit call:
+    //   his             histogram to fit
+    //   fitrange[2]     lo and hi boundaries of fit range
+    //   startvalues[4]  reasonable start values for the fit
+    //   parlimitslo[4]  lower parameter limits
+    //   parlimitshi[4]  upper parameter limits
+    //   fitparams[4]    returns the final fit parameters
+    //   fiterrors[4]    returns the final fit errors
+    //   ChiSqr          returns the chi square
+    //   NDF             returns ndf
+ 
+    Int_t i;
+    Char_t FunName[100];
+ 
+    sprintf(FunName,"Fitfcn_%s",his->GetName());
+ 
+    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
+    if (ffitold) delete ffitold;
+ 
+    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
+    ffit->SetParameters(startvalues);
+    ffit->SetParNames("Width","MP","Area","GSigma");
+ 
+    for (i=0; i<4; i++) {
+       ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
+    }
+ 
+    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot
+ 
+    ffit->GetParameters(fitparams);    // obtain fit parameters
+    for (i=0; i<4; i++) {
+       fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
+    }
+    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
+    NDF[0] = ffit->GetNDF();           // obtain ndf
+ 
+    return (ffit);              // return fit function
+ 
+ }
+ 
+ 
+ Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {
+ 
+    // Seaches for the location (x value) at the maximum of the
+    // Landau-Gaussian convolute and its full width at half-maximum.
+    //
+    // The search is probably not very efficient, but it's a first try.
+ 
+    Double_t p,x,fy,fxr,fxl;
+    Double_t step;
+    Double_t l,lold;
+    Int_t i = 0;
+    Int_t MAXCALLS = 10000;
+ 
+ 
+    // Search for maximum
+ 
+    p = params[1] - 0.1 * params[0];
+    step = 0.05 * params[0];
+    lold = -2.0;
+    l    = -1.0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = langaufun(&x,params);
+ 
+       if (l < lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-1);
+ 
+    maxx = x;
+ 
+    fy = l/2;
+ 
+ 
+    // Search for right x location of fy
+ 
+    p = maxx + params[0];
+    step = params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-2);
+ 
+    fxr = x;
+ 
+ 
+    // Search for left x location of fy
+ 
+    p = maxx - 0.5 * params[0];
+    step = -params[0];
+    lold = -2.0;
+    l    = -1e300;
+    i    = 0;
+ 
+    while ( (l != lold) && (i < MAXCALLS) ) {
+       i++;
+ 
+       lold = l;
+       x = p + step;
+       l = TMath::Abs(langaufun(&x,params) - fy);
+ 
+       if (l > lold)
+          step = -step/10;
+ 
+       p += step;
+    }
+ 
+    if (i == MAXCALLS)
+       return (-3);
+ 
+ 
+    fxl = x;
+ 
+    FWHM = fxr - fxl;
+    return (0);
+ }
+ 
+void langaus() {
+    // Fill Histogram
+   /*   Int_t data[100] = {10,20,50,3,10,5,2,6,11,18,18,55,90,141,255,323,454,563,681,
+                     737,821,796,832,720,637,558,519,460,357,291,279,241,212,
+                     153,164,139,106,95,91,76,80,80,59,58,51,30,49,23,35,28,23,
+                     22,27,27,24,20,16,17,14,20,12,12,13,10,17,7,6,12,6,12,4,
+                     9,9,10,3,4,5,2,4,1,5,5,1,7,1,6,3,3,3,4,5,4,4,2,2,7,2,4};
+    TH1F *hSNR = new TH1F("snr","Signal-to-noise",400,0,400);
+   */
+   // for (Int_t i=0; i<100; i++) hSNR->Fill(i,data[i]);
+   TH1D * hSNR;
+
+    Double_t graph_x[30], graph_y[30], graph_xerr[30], graph_yerr[30];
+    Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+   
+    
+    // hSNR = h_beamSignal_b0[13];
+    // Fitting SNR histo
+    printf("Fitting...\n");
+    TCanvas * c1 = new TCanvas("c1","c1", 800, 600);
+    
+    // Setting fit range and start values
+    Double_t fr[2];
+    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
+    Double_t chisqr;
+    Int_t    ndf;
+    TF1 *fitsnr;
+    Double_t SNRPeak, SNRFWHM;
+    char rootfilename[50] = "";
+    char saveplotname[50] = "";
+    char fout_mpv_name[50] = "";
+
+    int j = 1;
+     Double_t norm;
+     ofstream myfile, fout_mpv;
+     myfile.open ("MPVcorrection_proton0mm05.txt");
+     sprintf(fout_mpv_name, "jobs%i/plots/mpv.txt",j);
+     fout_mpv.open (fout_mpv_name);
+
+    for (int i = 0; i<26;i++){
+
+      sprintf(rootfilename, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile *rootFile = new TFile(rootfilename,"OPEN");
+
+      TH1D * hSNR = (TH1D*)rootFile->Get("h_spratio");
+      norm = hSNR->GetEntries();
+      hSNR->Scale(1/norm);
+      
+      fr[0]=hSNR->GetMean() - 2*hSNR->GetRMS();
+      fr[1]=hSNR->GetMean() + 2*hSNR->GetRMS();
+      //fr[1] = 2.0;
+      
+      pllo[0]=0.00001; pllo[1]=0.50; pllo[2]=0.001; pllo[3]=0.01; 
+      plhi[0]=0.1; plhi[1]=1.5; plhi[2]=0.055; plhi[3]=0.1;
+      sv[0]=0.05;
+      //sv[1]=2.1;
+      sv[2]=0.01;
+      sv[3]=0.02;
+      sv[1] = hSNR->GetMean();
+    
+      fitsnr = langaufit(hSNR,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);
+
+      
+      langaupro(fp,SNRPeak,SNRFWHM);
+      graph_x[i] = beta_proton[i];
+      graph_y[i] = fp[1];
+      graph_yerr[i] = sqrt(fitsnr->GetParError(1)*fitsnr->GetParError(1)+0.012*0.012*graph_y[i]*graph_y[i]);
+      fout_mpv << graph_x[i] << " " << graph_y[i] << " " << graph_yerr[i] << endl;
+      printf("Fitting done\nPlotting results...\n");
+      myfile << i << " " << graph_y[i] << endl;
+      // Global style settings
+      gStyle->SetOptStat(1111);
+      gStyle->SetOptFit(111);
+      //gStyle->SetLabelSize(0.03,"x");
+      //gStyle->SetLabelSize(0.03,"y");
+      
+      // hSNR->GetXaxis()->SetRange(0,70);
+      hSNR->Draw();
+      fitsnr->Draw("lsame");
+      c1->Update();
+      sprintf(saveplotname, "jobs%i/plots/runjob%i001_h_spratio.pdf",j,i);
+      c1->SaveAs(saveplotname);
+      //c1->WaitPrimitive();
+      hSNR->Delete();
+      rootFile->Close();
+    }
+     TGraphErrors * graph_1 = new TGraphErrors(23,graph_x, graph_y, 0, graph_yerr);
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+    graph_1->Draw("A*");
+    c2->SaveAs("MPVcorrection_proton0mm05.pdf");
+    myfile.close();
+    fout_mpv.close();
+ }
diff --git a/proton/plots.C b/proton/plots.C
new file mode 100644
index 0000000..12658ab
--- /dev/null
+++ b/proton/plots.C
@@ -0,0 +1,65 @@
+{
+
+
+   Double_t beta_proton[30] = {0.308525262, 0.340993523, 0.366173485, 0.386743776, 0.404197708, 0.4194131, 0.43294541, 0.445179695, 0.456345494, 0.466616582, 0.476154213, 0.48502264, 0.493348242, 0.501186212, 0.50863865, 0.515744144, 0.522562549, 0.52911069, 0.535379917, 0.541397728, 0.549575745, 0.557428612, 0.564849395, 0.571867977, 0.578541117, 0.584900169};
+    Double_t beta_helium[30] = {0.316661966, 0.34727202, 0.371222186, 0.391037227, 0.408018455, 0.422922098, 0.436235455, 0.44827542, 0.459299095, 0.469441244, 0.478845524, 0.487649369, 0.495886825, 0.503656244, 0.510990953, 0.517959457, 0.52457247, 0.530888884, 0.536925849, 0.54269899, 0.550671248, 0.55845178, 0.565814653, 0.572798702, 0.579448698, 0.585785313};
+    Double_t beta_carbon[30] = {0.407931067, 0.448122448, 0.479020432, 0.504000457, 0.524971648, 0.543076553, 0.559041917, 0.573329576, 0.586254563, 0.598061846, 0.608917559, 0.618952311, 0.62829287, 0.637039726, 0.645286945, 0.65311609, 0.660570813, 0.667689852, 0.67446932, 0.680943928, 0.689677353, 0.698000799, 0.70580765, 0.71315081, 0.720086739, 0.726650602};
+    Double_t beta_oxygen[30] = {0.43638582, 0.479000679, 0.51134888, 0.537325331, 0.559061572, 0.57764689, 0.594123989, 0.608730698, 0.621997639, 0.63402408, 0.645050809, 0.655226738, 0.664724227, 0.673475951, 0.681810969, 0.689681788, 0.697243281, 0.704219104, 0.710968918, 0.717442538, 0.726111033, 0.734356548, 0.742073831, 0.749383281, 0.756156282, 0.762562424};
+  
+  TGraph * gr_letratio = new TGraph();
+
+    char  finname[50] = "";
+    char fout_mean_name[50] = "";
+    char fout_median_name[50] = "";
+    char fout_mpv_name[50] = "";
+
+    ofstream fout_mean, fout_median;
+
+    int j = 1;
+      sprintf(fout_mean_name, "jobs%i/plots/mean.txt",j);
+      sprintf(fout_median_name, "jobs%i/plots/median.txt",j);
+      sprintf(fout_mpv_name, "jobs%i/plots/mpv.txt",j);
+
+    
+    fout_median.open (fout_median_name);
+    fout_mean.open (fout_mean_name);
+
+    
+    Double_t x, q;
+    q = 0.5; // 0.5 for "median"
+    
+    
+    //  for (int j = 0; j<5;j++){
+    for (int i = 0; i<26;i++){
+
+      sprintf(finname, "jobs%i/runjob%i001_eventdata_out.root",j,i);
+      TFile * f = new TFile(finname,"OPEN");
+      TH1F * h = (TH1F*)f->Get("h_spratio");
+      cout << i << " " << h->GetMean() << " +/- " << h->GetRMS() << endl;
+      fout_mean  << beta_proton[i] << " " << h->GetMean() << " " << 0.012*h->GetMean() << endl;
+      h->ComputeIntegral(); // just a precaution
+      h->GetQuantiles(1, &x, &q);
+      std::cout << "median = " << x << std::endl;
+      fout_median  << beta_proton[i] << " " << x << " " << 0.012*x << endl;
+
+      
+      h->Delete();
+      
+      f->Close();
+      
+    }
+    //  }
+    TGraphErrors * graph_1 = new TGraphErrors(fout_mean_name,"%lg %lg %lg" );
+    TGraphErrors * graph_2 = new TGraphErrors(fout_median_name,"%lg %lg %lg" );
+    TGraphErrors * graph_3 = new TGraphErrors(fout_mpv_name,"%lg %lg %lg" );
+
+    TCanvas * c2 = new TCanvas("c2","c2", 800, 600);
+    c2->cd();
+    graph_3->Draw("A*");
+    graph_1->Draw("*");
+    graph_2->Draw("*");
+    fout_mean.close();
+    fout_median.close();
+
+}
+
diff --git a/proton/runJob.sh b/proton/runJob.sh
new file mode 100755
index 0000000..73ea524
--- /dev/null
+++ b/proton/runJob.sh
@@ -0,0 +1,63 @@
+#!/bin/bash
+
+TYPE="PROTON"
+energy=(0.04812 0.05982 0.07003 0.07917 0.08753 0.09530 0.10261 0.10956 0.11620 0.12257 0.12872 0.13465 0.14041 0.14601 0.15150 0.15689 0.16221 0.16746 0.17262 0.17770 0.18481 0.19187 0.19876 0.20548 0.21206 0.21851); #GeV/u, protons
+#energy=(0.08883 0.11058 0.12979 0.14713 0.16309 0.17801 0.19213 0.20560 0.21852 0.23098 0.24303 0.25471 0.26608 0.27719 0.28810 0.29887 0.30952 0.32007 0.33048 0.34077 0.35522 0.36964 0.38378 0.39766 0.41132 0.42477); #GeV/u, carbon
+#energy=(0.05057 0.06186 0.07173 0.08064 0.08885 0.09652 0.10376 0.11064 0.11723 0.12355 0.12964 0.13555 0.14127 0.14684 0.15226 0.15756 0.16273 0.16780 0.17277 0.17764 0.18456 0.19154 0.19836 0.20503 0.21157 0.21798); #GeV/u, helium
+
+#THICKNESS=(0.025 0.050 0.075 0.100 1.25 0.150 0.2500 0.500) #cm
+THICKNESS=(0.025 0.050 0.100 0.150 0.500 1.000 0.125) #cm
+
+HOME=/work/leverington/fluka_gfortran7/project/proton
+ZED=2.
+ATN=4.
+
+for j in {6..6}
+do 
+    mkdir -p $HOME/jobs$j #make the directory if it doesn't exist
+    JOB_HOME=$HOME/jobs$j 
+    rm -r  $JOB_HOME/* #clean up the directory
+    cd $JOB_HOME
+    
+    
+    for i in {0..25} #26
+    do
+	# make a copy of the template
+	cp $HOME/runjobtemplate.inp  $JOB_HOME/runjob$i.inp
+	touch $JOB_HOME/runjob$i.inp
+	
+	#replace template placeholders with values
+	sed -i 's/ENERGYIN/'"-${energy[i]}"'/g' $JOB_HOME/runjob$i.inp #the negative sign indicates energy, not momentum
+	sed -i 's/TYPE/'"$TYPE"'/g' $JOB_HOME/runjob$i.inp  #beam particle type
+	sed -i 's/THICK/'"${THICKNESS[j]}"'/g' $JOB_HOME/runjob$i.inp  #thickness of polystyrene layer
+	sed -i 's/ZED/'"$ZED"'/g' $JOB_HOME/runjob$i.inp  #atomic charge
+	sed -i 's/ATN/'"$ATN"'/g' $JOB_HOME/runjob$i.inp  #atomic number
+
+	echo "Job: $j-$i   Ion type: $TYPE    Energy: ${energy[i]} GeV/u    Plane thickness: ${THICKNESS[j]} cm"
+	echo ""
+	
+	touch $JOB_HOME/fluka$j-$i.sh  #script to submit to the queue
+	echo "#!/bin/bash" >> $JOB_HOME/fluka$j-$i.sh
+	echo ". /local/env.sh" >> $JOB_HOME/fluka$j-$i.sh
+#	echo ". /cvmfs/lhcb.cern.ch/lib/lhcb/LBSCRIPTS/LBSCRIPTS_v9r2p2/InstallArea/scripts/LbLogin.sh -c x86_64-slc6-gcc7-opt" >> $JOB_HOME/fluka$j-$i.sh
+	echo "source /home/lhcb/leverington/setup.sh" >> $JOB_HOME/fluka$j-$i.sh
+        
+        #execute this command(s)
+	echo "cd "$JOB_HOME >> $JOB_HOME/fluka$j-$i.sh
+	echo "$FLUPRO/flutil/rfluka -e $FLUPRO/flutil/flukadpm3 -N0 -M1 runjob"$i >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "/work/leverington/readfluka/tools/eventdat2root runjob"$i"001_eventdata"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "sleep 2"  >> $JOB_HOME/fluka$j-$i.sh
+	echo "$HOME/eventdatroot runjob"$i"001_eventdata "$j  >> $JOB_HOME/fluka$j-$i.sh
+	
+	qsub -l os=slc6 -l hio=10 -l ujl=40 -cwd -j yes $JOB_HOME/fluka$j-$i.sh
+
+#	$FLUPRO/flutil/rfluka -e /home/leverington/software/fluka/flutil/flukadpm3 -N0 -M1 $JOB_HOME/runjob$i  #run 
+#	sleep 2
+#	eventdat2root runjob$i\001_eventdata  #convert to root file
+#	sleep 2
+#	$HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+ 	sleep 1
+    done
+    sleep 30
+done
diff --git a/proton/runanalysis.sh b/proton/runanalysis.sh
new file mode 100755
index 0000000..56f2beb
--- /dev/null
+++ b/proton/runanalysis.sh
@@ -0,0 +1,17 @@
+#!/bin/bash
+
+HOME=/work/leverington/fluka_gfortran7/project/proton
+
+
+for j in {0..5}
+do
+    JOB_HOME=$HOME/jobs$j
+    cd $JOB_HOME
+
+	     for i in {0..25} #26
+	     do
+		 
+		 $HOME/eventdatroot runjob$i\001_eventdata $j #run analysis code
+		 sleep 0.1 
+	     done
+done
diff --git a/proton/runjobtemplate.inp b/proton/runjobtemplate.inp
new file mode 100644
index 0000000..1b2c6b8
--- /dev/null
+++ b/proton/runjobtemplate.inp
@@ -0,0 +1,95 @@
+* ..+....1....+....2....+....3....+....4....+....5....+....6....+....7...
+TITLE
+FLUKA Course Exercise
+*
+* use names everywhere and free format for geometry
+DEFAULTS                                                              HADROTHE
+*
+* beam definitions
+BEAM      ENERGYIN  0.0       -1.7      -0.8      -0.8                TYPE
+*HI-PROPE          6.       12.
+BEAMPOS                             -0.1
+DELTARAY       1E-06       0.0       0.0  HYDROGEN  @LASTMAT          NOPRINT
+PART-THR  -0.00001  PROTON    PROTON                 0.0   
+*
+* Geometry
+* --------
+GEOBEGIN                                                              COMBNAME
+    0    0             Cylindrical Target
+*
+* Bodies
+* ------
+*
+* Blackhole to include geometry
+SPH BLK        0.0 0.0 0.0 10000.
+* Void sphere
+SPH VOID       0.0 0.0 0.0 1000.
+* Infinite cylinder
+ZCC TARG       0.0 0.0 5.
+*  planes cutting the cylinder
+XYP ZTlow      0.0
+XYP T1seg      THICK
+XYP T2seg      2.
+XYP ZThigh     10.
+*  additional plane for scoring
+XYP ZTscor     9.9999
+END
+*
+* Regions
+* -------
+*
+* Blackhole
+BLKHOLE      5 +BLK -VOID
+*
+* Target segment 1
+TARGS1       5 +TARG -ZTlow +T1seg
+* Target segment 2
+TARGS2       5 +TARG -T1seg +T2seg
+* Target segment 3
+TARGS3       5 +TARG -T2seg +ZTscor
+* Scoring region
+TARGS4       5 +TARG -ZTscor +ZThigh
+*
+* Air around target
+INAIR        5 | +VOID -TARG
+               | +VOID +ZTlow
+               | +VOID -ZThigh
+END
+*
+GEOEND
+*
+* Materials definition
+* --------------------
+MATERIAL                        0.001965                              CO2
+COMPOUND          1.    CARBON        2.    OXYGEN                    CO2
+*
+* Assign materials
+* ----------------
+ASSIGNMA    BLCKHOLE   BLKHOLE
+ASSIGNMA    POLYSTYR    TARGS1
+ASSIGNMA      VACUUM    TARGS2
+ASSIGNMA      VACUUM    TARGS3
+ASSIGNMA      VACUUM    TARGS4
+ASSIGNMA      VACUUM     INAIR
+*
+* Exercise: heavy ions
+* --------------------
+*
+*  charge spectrum of ions
+USRYIELD       1422.  HEAVYION      -68.    TARGS1    TARGS2        1.FragZ1
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS2    TARGS3        1.FragZ2
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+USRYIELD       1422.  HEAVYION      -68.    TARGS3    TARGS4        1.FragZ3
+USRYIELD         9.5       2.5        7.       90.       0.0        3. &
+*
+*  LET in water of ions and charged particles
+USRYIELD       2223.  HEAVYION      -69.    TARGS3    TARGS4        1.LETHI
+USRYIELD         20.       0.0      200.       9.5       2.5     2703. &
+USRYIELD       2223.  ALL-CHAR      -69.    TARGS3    TARGS4        1.LETCh
+USRYIELD         20.       0.0      200.       9.5      -2.5     2703. &
+EVENTDAT        -21.                                                  eventdata
+*
+RANDOMIZ          1.
+START         10000.                                     0.0
+STOP