diff --git a/geant4/LEMuSR/macro.c b/geant4/LEMuSR/macro.c
new file mode 100644
index 0000000..3a7d0e5
--- /dev/null
+++ b/geant4/LEMuSR/macro.c
@@ -0,0 +1,597 @@
+gStyle->SetPalette(1);
+gStyle->SetErrorX(0);
+gStyle->SetOptStat(0);
+
+TFile *trans = TFile::Open("SDInnerScint.root");
+TTree *tree= (TTree*)trans->Get("tree");
+
+TFile *foil = TFile::Open("SDCryoE0.root");
+TTree *plane= (TTree*)foil->Get("tree");
+
+
+
+TCanvas *ec1;
+TCanvas *ec2;
+
+
+double meanMU, meanmuplus, ie ;
+
+  double  implE;
+
+double Bfield=0;
+
+TGraph *g1, *g2;
+
+
+void esquisse(){
+gStyle->SetPalette(1);
+  ec1= new TCanvas();
+  ec2= new TCanvas();
+
+  ec1->SetGrid();
+  ec1->Divide(2,2);
+
+  TString poscut="positionz==14";
+
+  ec1->cd(1);
+  plane->Draw("positionz>>muPOS(1500,14,14.00003)");
+
+  ec1->cd(2);
+  plane->Draw("muon",poscut);
+
+  ec1->cd(3);
+  plane->Draw("globaltime>>hTOF(150)",poscut,"E1");
+
+  ec1->cd(4);
+  plane->Draw("kenergy+edeposit>>hKE(150)",poscut,"E1");
+
+
+
+  ec2->SetGrid();
+  ec2->Divide(2,2);
+  
+  TFile *foil = TFile::Open("SDCryoE0.root");
+  TTree *plane= (TTree*)foil->Get("tree");
+  
+  ec2->cd(1);
+  plane->Draw("positiony:positionx>>mxy(100,100)",poscut+"&&muon==1","zcol");
+  
+  ec2->cd(2);
+  plane->Draw("positiony:positionx>>Mxy(100,100)",poscut+"&&muonium==1","zcol");
+  
+  ec2->cd(3);
+  plane->Draw("globaltime>>mpg(100)",poscut+"&&muon==1","E1");
+  
+  ec2->cd(4);
+  plane->Draw("globaltime>>mg(100)",poscut+"&&muonium==1","E1");
+ 
+  
+    
+}
+
+
+void moyennes(int ID=-1)
+{
+
+gStyle->SetPalette(1);
+gStyle->SetErrorX(0);
+gStyle->SetOptStat(11);
+
+TCanvas *ec1;
+TCanvas *ec2;
+
+
+double meanMU, meanmuplus, ie ;
+
+  double  implE;
+
+  TString cut = "positionz==14";
+
+  if(ID!=-1)
+    {
+      cut=cut+"&&runID=="+FS(ID);
+    }
+
+  plane->Draw("globaltime>>hmuplus(1000,0,500)",cut+"&&muon==1");
+  meanmuplus = hmuplus->GetMean();
+
+  plane->Draw("globaltime>>hMU(1000,0,500)",cut+"&&muonium==1");
+  meanMU = hMU->GetMean();
+   
+  //  plane->Draw("Init_Energy>>hIE");
+  plane->Draw("kenergy+edeposit>>hIE(500,0,5)",cut+"&&muon==1");
+
+   ie=hIE->GetMean();
+
+
+  ofstream outf;
+
+  outf.open("means.dat");
+  
+  if(outf.is_open())
+    {
+      std::cout<<meanmuplus << " " << meanMU<<" " << ie <<std::endl;
+       outf<<meanmuplus << " " << meanMU<<" " << ie <<std::endl;
+      outf.close();
+      std::cout<<"Output file means.dat created."<<std::endl;
+    }
+  else
+    {
+      cout << "Error opening file";
+      return 0;
+    }
+  
+}
+
+
+void mtime(TString particle="mu+",TString grids="ON", TString guards="ON"){
+  
+ gDirectory->DeleteAll();
+
+ int ID;
+
+  std::cout<<"Run? (-1 for all): " <<std::endl;
+  std::cin>>ID;
+  std::cout<<"Magnetic Field [Gauss]?" <<std::endl;
+  std::cin>>Bfield;
+  //  std::cout<<"Implantation Energy?" <<std::endl;
+  //  std::cin>>implE;
+
+  int b_f=1;
+  if (Bfield==0)b_f=0;
+
+#define GMUON 0.01355
+#define TAU_MUON 2.197 
+#define ASYM 0.28
+  Int_t i,j;
+  double r;
+  char  func_str[256];
+  char  text_str[256];
+ 
+ GetValues();
+
+  
+  sprintf(func_str, "exp(-x/%f)", TAU_MUON);
+  TF1 *f1 = new TF1("f1", func_str, 0., 20.);
+  //sprintf(func_str, "exp(-x/%f)", TAU_MUON);
+  TF1 *f2 = new TF1("f2", func_str, 0., 20.);
+  TF1 *fitdecay = new TF1("fitdecay", TFdecayfit, meanmuplus, 10000., 2);
+  TF1 *fitasym1  = new TF1("fitasym1", TFasymuonfit, meanmuplus, 10000., 3);
+  TF1 *fitasym2  = new TF1("fitasym2", TFasymuonfit, meanmuplus, 10000., 3);
+
+  /*----------- initialise fit hforw and hback ----------------*/
+  fitdecay->SetParNames("N(0)","MuonLifeTime");
+  fitdecay->SetParameters(0, 450.);
+  fitdecay->SetParameters(1, 2.2);
+
+  // convert frequence from Hz 1/s to 1/ns (10^-9 Hz)
+  double fz = 13.*b_f; 
+  double ampl = 0.2; 
+  
+  if(particle=="Mu")
+    {
+      fz = -1.39e3*b_f;
+      ampl = 0.16;
+    }
+  fitasym1->SetParNames("Amplitude","Frequency [kHz]","Offset/Pi");
+  fitasym1->SetParameters(ampl,fz,0*b_f);
+  fitasym1->SetLineColor(46);
+
+
+  fitasym2->SetParNames("Amplitude","Frequency [kHz]","Offset/Pi");
+  fitasym2->SetParameters(ampl,fz,0.5*b_f);
+  fitasym2->SetLineColor(50);
+
+
+  //===================================================================
+  //===================================================================
+
+ /*---------------------------------------------------------*/
+ TH1F *hRight = new TH1F("hRight", "hRight",100, 0., 2400.);
+ TH1F *hLeft = new TH1F("hLeft", "hLeft",100, 0., 2400.);
+ TH1F *LRdiff = new TH1F("LRdiff", "LRdiff",100, 0., 2400.);
+ TH1F *LRsum  = new TH1F("LRsum",  "LRsum",100, 0., 2400.);
+ TH1F *LRasym = new TH1F("LRasym", "LRasym",100, 0., 2400.);
+ /*---------------------------------------------------------*/
+ TH1F *hTop = new TH1F("hTop", "hTop",100, 0., 2400.);
+ TH1F *hBottom = new TH1F("hBottom", "hBottom",100, 0., 2400.);
+ TH1F *BTdiff = new TH1F("BTdiff", "BTdiff",100, 0., 2400.);
+ TH1F *BTsum  = new TH1F("BTsum",  "BTsum",100, 0., 2400.);
+ TH1F *BTasym = new TH1F("BTasym", "BTasym",100, 0., 2400.);
+ /*---------------------------------------------------------*/
+ if(particle=="Mu")
+   { 
+ hRight->SetBins( 100, 0.,  meanMU+200.);
+ hLeft->SetBins(  100, 0.,  meanMU+200.);
+ LRdiff->SetBins( 100, 0.,  meanMU+200.);
+ LRsum->SetBins(  100, 0.,  meanMU+200.);
+ LRasym->SetBins( 100, 0.,  meanMU+200.);
+ /*---------------------------------------------------------*/
+ hTop->SetBins(100, 0., 500.);
+ hBottom->SetBins(100, 0., 500.);
+ BTdiff->SetBins(100, 0., 500.);
+ BTsum->SetBins(100, 0., 500.);
+ BTasym->SetBins(100, 0., 500.);
+   }
+ /*---------------------------------------------------------*/
+   
+
+
+ hLeft->Sumw2();
+ hRight->Sumw2();
+ hTop->Sumw2();
+ hBottom->Sumw2();
+ LRdiff->Sumw2();
+ LRsum->Sumw2();
+ LRasym->Sumw2();
+ BTdiff->Sumw2();
+ BTsum->Sumw2();
+ BTasym->Sumw2();
+
+  //===================================================================
+  //===================================================================
+  
+  //  moyennes();
+  TString cut2;
+  cut2 = "";
+  if(ID==-1)
+    {  
+      cut2 = "";
+    }
+  else
+    {
+      cut2 = "positron.runID=="+FS(ID)+"&&";
+    }
+
+  tree->Draw("positron.globaltime>>hLeft",cut2+"positron.Left==1","");
+  tree->Draw("positron.globaltime>>hRight",cut2+"positron.Right==1","");
+  tree->Draw("positron.globaltime>>hBottom",cut2+"positron.Bottom==1","");
+  tree->Draw("positron.globaltime>>hTop",cut2+"positron.Top==1","");
+   
+ 
+  LRsum->Add(hLeft,hRight,1,1);
+  LRdiff->Add(hLeft,hRight,1,-1);
+  LRasym->Divide(LRdiff,LRsum);
+
+  hLeft->SetLineColor(1);
+  hTop->SetLineColor(2);
+  hRight->SetLineColor(3);
+  hBottom->SetLineColor(4);
+
+
+ 
+
+  //=============================================================================
+  //                              LEFT-RIGHT
+  //=============================================================================
+
+
+
+
+  c1=new TCanvas("c1","LR Counts");
+  c1->SetGrid();
+  c1->SetFillColor(0);
+  c1->SetFrameFillColor(0);
+  c1->SetFrameBorderMode(0);
+  
+  hLeft->SetTitle("Counts in Left and Right Scintillators");
+  hLeft->GetYaxis()->SetTitle("Counts (arbitrary scale)");
+  hLeft->GetXaxis()->SetTitle("Time [ns]");
+  hLeft->SetMarkerStyle(3);
+  hLeft->SetMarkerColor(4);
+  hLeft->SetLineColor(4);
+
+  hLeft->Draw("E1lp");
+
+
+
+  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
+  pt->SetName("title");
+  pt->SetTextSize(0.032);
+   pt->SetBorderSize(0);
+  pt->SetFillColor(0);
+  text = pt->AddText("");
+  pt->Draw();
+
+  
+  hRight->SetMarkerStyle(26);
+  hRight->SetMarkerColor(2);
+  hRight->SetLineColor(2);
+  hRight->Draw("E1lpsame");
+
+  TOF_plot();
+  
+  leg = new TLegend(0.2,0.75,0.43,0.88);
+  leg->SetFillColor(0);
+  leg->SetBorderSize(1);
+  
+  leg->AddEntry(hLeft,"Left ","p");
+  leg->AddEntry(hRight,"Right ","p");
+    leg->AddEntry(g1,"#mu+ TOF ("+FS(meanmuplus)+" ns)","l");
+    leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
+  leg->Draw();
+  
+
+
+
+  //++++++++++++++++++++++++++++++++++
+  //++++++++++++++++++++++++++++++++++
+
+  c3=new TCanvas("c3","LR Asymmetry");
+  c3->SetGrid();
+  c3->SetFillColor(0);
+  c3->SetFrameFillColor(0);
+  c3->SetFrameBorderMode(0);
+  TString title = "Left/Right Asymmetry for "+FS(implE)+" KeV impl. energy. (B="+FS(Bfield)+"G)";
+
+
+  LRasym->SetTitle(title);
+  LRasym->GetXaxis()->SetTitle("Time [ns]");
+  LRasym->GetYaxis()->SetTitle("Asymmetry");
+  LRasym->GetYaxis()->SetRangeUser(-1,1);
+  LRasym->SetMarkerStyle(2);
+  LRasym->SetMarkerColor(38);
+  LRasym->SetLineColor(38);
+
+  LRasym->Fit("fitasym1");
+  //  gStyle->SetOptStat(0);
+  gStyle->SetOptFit(111);  
+  LRasym->Draw("e1");
+
+ 
+  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
+  pt->SetName("title");
+  pt->SetBorderSize(0);
+  pt->SetTextSize(0.032);
+   pt->SetFillColor(0);
+  text = pt->AddText("");
+  pt->Draw();
+
+
+
+  
+  TOF_plot();
+
+  
+
+  leg = new TLegend(0.2,0.75,0.43,0.88);
+  leg->SetFillColor(0);
+  leg->SetBorderSize(1);
+
+  leg->AddEntry(LRasym,"L-R Asym","p");
+  leg->AddEntry(g1,"#mu+ TOF ("+FS(meanmuplus)+" ns)","l");
+    leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
+  leg->Draw();
+
+
+  //=============================================================================
+  //                             Bottom-Top
+  //=============================================================================
+gStyle->SetPalette(1);
+gStyle->SetOptStat(11);
+gStyle->SetErrorX(0);
+
+
+
+  BTsum->Add(hBottom,hTop,1,1);
+  BTdiff->Add(hBottom,hTop,1,-1);
+  BTasym->Divide(BTdiff,BTsum);
+  BTasym->SetMarkerStyle(3);
+
+
+  c4=new TCanvas("c4","BT Counts");
+  c4->SetGrid();
+  c4->SetFillColor(0);
+   c4->SetFillColor(0);
+  c4->SetFrameFillColor(0);
+  c4->SetFrameBorderMode(0);
+  hTop->SetTitle("Counts in Top and Bottom Scintillators");
+
+  hTop->SetMarkerStyle(26);
+  hTop->SetMarkerColor(2);
+  hTop->SetLineColor(2);
+ hTop->GetYaxis()->SetTitle("Counts (arbitrary scale)");
+  hTop->GetXaxis()->SetTitle("Time [ns]");
+  hTop->Draw("lp");
+
+
+  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
+  pt->SetName("title");
+  pt->SetBorderSize(0);
+  pt->SetTextSize(0.032);
+   pt->SetFillColor(0);
+  text = pt->AddText("");
+  pt->Draw();
+
+
+
+  hBottom->SetMarkerStyle(3);
+  hBottom->SetMarkerColor(4);
+  hBottom->SetLineColor(4);
+  hBottom->Draw("lpsame");
+
+  TOF_plot();
+
+  leg = new TLegend(0.2,0.75,0.43,0.88);
+  leg->SetFillColor(0);
+  leg->SetBorderSize(1);
+  leg->AddEntry(hBottom,"Bottom ","p");
+  leg->AddEntry(hTop,"Top ","p");
+  leg->AddEntry(g1,"#mu+ TOF("+FS(meanmuplus)+" ns)","l");
+  leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
+  leg->Draw();
+  
+
+
+  //++++++++++++++++++++++++++++++++++
+  //++++++++++++++++++++++++++++++++++
+
+
+  c6=new TCanvas("c6","BT Asymmetry");
+  c6->SetGrid();
+  c6->SetFillColor(0);
+  c6->SetFrameFillColor(0);
+  c6->SetFrameBorderMode(0);
+
+
+
+  title = "Bottom/Top Asymmetry for "+FS(implE)+" KeV impl. energy. (B="+FS(Bfield)+"G)";
+  BTasym->SetTitle(title);
+  BTasym->GetXaxis()->SetTitle("Time [ns]");
+  BTasym->GetYaxis()->SetTitle("Asymmetry");
+  BTasym->GetYaxis()->SetRangeUser(-1,1);
+  BTasym->SetMarkerColor(31);
+  BTasym->SetLineColor(31);
+
+
+  //  BTasym->GetXaxis()->SetRangeUser(0,5000);
+  BTasym->Fit("fitasym2");
+  gStyle->SetOptStat(11);
+  gStyle->SetOptFit(111);  
+
+  BTasym->Draw("p");
+  pt = new TPaveText(0.01,0.91,0.6,0.99,"blNDC");
+  pt->SetName("title");
+  pt->SetTextSize(0.032);
+  pt->SetBorderSize(0);
+  pt->SetFillColor(0);
+  text = pt->AddText("");
+  pt->Draw();
+
+  TOF_plot();
+
+  leg = new TLegend(0.2,0.75,0.43,0.88);
+  leg->SetFillColor(0);
+  leg->SetBorderSize(1);
+
+  leg->AddEntry(BTasym,"B-T Asym","p");
+    leg->AddEntry(g1,"#mu+ TOF("+FS(meanmuplus)+" ns)","l");
+    leg->AddEntry(g2,"Mu  TOF("+FS(meanMU)+" ns)","l");
+ leg->Draw();
+  
+
+  //  c6->SetEditable(0);
+
+
+
+}
+
+
+
+
+TString flt2string(const float& number)
+{
+  ostringstream oss;
+  oss << number;
+  return oss.str();
+}
+
+TString FS(const float& number)
+{
+  ostringstream oss;
+  oss << number;
+  return oss.str();
+}
+
+saveall(  TString save_ans="")
+{
+
+ 
+      c1->SaveAs("LRcounts"+save_ans+".eps");
+      c3->SaveAs("LRasym"+save_ans+".eps");
+      c4->SaveAs("BTcounts"+save_ans+".eps");
+      c6->SaveAs("BTasym"+save_ans+".eps");
+      /* c1->SaveAs("LRcounts"+FS(implE)+save_ans+".gif");
+      c3->SaveAs("LRasym"+FS(implE)+save_ans+".gif");
+      c4->SaveAs("BTcounts"+FS(implE)+save_ans+".gif");
+      c6->SaveAs("BTasym"+FS(implE)+save_ans+".gif");*/
+
+
+}
+
+void TOF_plot()
+{
+    Double_t a[2],b[2];
+  
+  a[0]=meanmuplus;a[1]=meanmuplus;
+  b[0]=-1;b[1]=2500;
+  g1= new TGraph(2,a,b);
+  a[0]=meanMU;a[1]=meanMU;
+  b[0]=-1;b[1]=2500;
+  g2= new TGraph(2,a,b);
+  g1->SetLineColor(1);
+  g1->SetLineWidth(2);
+  g1->Draw("");
+  g2->SetLineColor(2);
+  g2->SetLineWidth(2);
+  g2->Draw("");
+}
+
+
+void GetValues()
+{
+ 
+
+
+ ifstream inf;
+  inf.open("means.dat");
+  while(!inf.is_open())
+    {
+      std::cout<<"No 'means.dat' file.\nAttempt to create it."<<std::endl;
+      
+      moyennes();
+      inf.open("means.dat");
+    }
+
+  inf >> meanmuplus >> meanMU >> implE;// ie;
+  inf.close();
+
+}
+
+
+
+Double_t MuonLifeFunc(Double_t *x, Double_t *par) {
+  // two fit parameters: the normalisation
+  //                     the muon life time
+  return par[0]*TMath::Exp(-x[0]/par[1]);
+}
+
+Double_t TFasymuonfit(Double_t *x, Double_t *par) {
+  // one fit parameter: the asymmetry
+  Double_t fitval, arg;
+  arg=-x[0]*Bfield*6.28*par[1]*1.e-6+3.14*par[2];
+  fitval = par[0]*TMath::Cos(arg);
+  return fitval;
+}
+
+
+Double_t TFdecayfit(Double_t *x, Double_t *par) {
+ // 
+ Double_t fitval;
+ fitval = MuonLifeFunc(x, par);
+ return fitval;
+}
+
+
+format()
+{
+
+  for(Int_t i==0; i<3; i++)
+    {
+      CASYM->cd(i)->SetFrameBorderMode(0);
+      CASYM->cd(i)->SetFrameBorderSize(0);
+      CASYM->cd(i)->SetGrid();
+      CASYM->cd(i)->SetEditable(1);
+   }
+
+
+format(TCanvas cv)
+{
+      cv->SetFrameBorderMode(0);
+      cv->SetFrameBorderSize(0);
+      cv->SetGrid();
+      cv->SetEditable(0);
+   }
+
+}
+