leverington/fluka_gfortran7_hit
1
1
Fork 0
Code Issues Pull Requests repo.project_board Releases Wiki Activity

first push

Browse Source
This commit is contained in:
Blake Leverington 2021-08-30 16:15:26 +02:00
parent 061348dedc
commit 090051977b
58 changed files with 5444 additions and 0 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

7
.gitignore vendored
View File

@ -66,3 +66,10 @@
*.out
*.app
# my stuff
**/job*
**/*.gz
**/*.tgz
**/*~
**/*.root

26
carbon/MPVcorrection_carbon0mm05.txt Normal file
View File

@ -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

38
carbon/Makefile Normal file
View File

@ -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

BIN
carbon/carbon0mm5_SPratio.pdf Normal file
View File

Binary file not shown.

BIN
carbon/eventdatroot Executable file
View File

Binary file not shown.

137
carbon/eventdatroot.c Normal file
View File

@ -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;
}

199
carbon/eventdatroot.h Normal file
View File

@ -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

348
carbon/langaus.C Normal file
View File

@ -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();
}

32
carbon/plots.C Normal file
View File

@ -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();
}
// }
}

64
carbon/runJob.sh Executable file
View File

@ -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

17
carbon/runanalysis.sh Executable file
View File

@ -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

95
carbon/runjobtemplate.inp Normal file
View File

@ -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

BIN
helium/MPVcorrection_helium0mm05.pdf Normal file
View File

Binary file not shown.

26
helium/MPVcorrection_helium0mm05.txt Normal file
View File

@ -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

38
helium/Makefile Normal file
View File

@ -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

BIN
helium/eventdatroot Executable file
View File

Binary file not shown.

143
helium/eventdatroot.c Normal file
View File

@ -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;
}

199
helium/eventdatroot.h Normal file
View File

@ -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

355
helium/langaus.C Normal file
View File

@ -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();
}

63
helium/plots.C Normal file
View File

@ -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();
}

63
helium/runJob.sh Executable file
View File

@ -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

17
helium/runanalysis.sh Executable file
View File

@ -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

95
helium/runjobtemplate.inp Normal file
View File

@ -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

25
oxygen/MPVcorrection_oxygen0mm05.txt Normal file
View File

@ -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

38
oxygen/Makefile Normal file
View File

@ -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

BIN
oxygen/eventdatroot Executable file
View File

Binary file not shown.

139
oxygen/eventdatroot.c Normal file
View File

@ -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;
}

199
oxygen/eventdatroot.h Normal file
View File

@ -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

348
oxygen/langaus.C Normal file
View File

@ -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();
}

BIN
oxygen/oxygen0mm5_SPratio.pdf Normal file
View File

Binary file not shown.

32
oxygen/plots.C Normal file
View File

@ -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();
}
// }
}

65
oxygen/runJob.sh Executable file
View File

@ -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

17
oxygen/runanalysis.sh Executable file
View File

@ -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

95
oxygen/runjobtemplate.inp Normal file
View File

@ -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

116
pion/MPVcorrection_proton0mm05 Normal file
View File

@ -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

26
pion/MPVcorrection_proton0mm05.txt Normal file
View File

@ -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

38
pion/Makefile Normal file
View File

@ -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

BIN
pion/eventdatroot Executable file
View File

Binary file not shown.

135
pion/eventdatroot.c Normal file
View File

@ -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;
}

199
pion/eventdatroot.h Normal file
View File

@ -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

343
pion/langaus.C Normal file
View File

@ -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();
}

32
pion/plots.C Normal file
View File

@ -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();
}
// }
}

64
pion/runJob.sh Executable file
View File

@ -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

17
pion/runanalysis.sh Executable file
View File

@ -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

95
pion/runjobtemplate.inp Normal file
View File

@ -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

352
proton/#langaus.C# Normal file
View File

@ -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();
}

117
proton/MPVcorrection_proton0mm05 Normal file
View File

@ -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

BIN
proton/MPVcorrection_proton0mm05.pdf Normal file
View File

Binary file not shown.

26
proton/MPVcorrection_proton0mm05.txt Normal file
View File

@ -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

38
proton/Makefile Normal file
View File

@ -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

BIN
proton/eventdatroot Executable file
View File

Binary file not shown.

135
proton/eventdatroot.c Normal file
View File

@ -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;
}

199
proton/eventdatroot.h Normal file
View File

@ -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

352
proton/langaus.C Normal file
View File

@ -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();
}

65
proton/plots.C Normal file
View File

@ -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();
}

63
proton/runJob.sh Executable file
View File

@ -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

17
proton/runanalysis.sh Executable file
View File

@ -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

95
proton/runjobtemplate.inp Normal file
View File

@ -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