Skip to content

A
AculUtils
Commit 5016cc63 authored by Vratislav Chudoba's avatar Vratislav Chudoba
Browse files
Options

Added class TELoss for calculation of energy losses.

parent 67b93633
Expand all Hide whitespace changes
Inline Side-by-side
Showing with 862 additions and 1990 deletions
AculData/AculDataCint.cpp deleted 100644 → 0
View file @ 67b93633
This diff is collapsed.
AculData/AculDataCint.h deleted 100644 → 0
View file @ 67b93633
/********************************************************************
* /home/dariak/Utilities/AculData/AculDataCint.h
* CAUTION: DON'T CHANGE THIS FILE. THIS FILE IS AUTOMATICALLY GENERATED
* FROM HEADER FILES LISTED IN G__setup_cpp_environmentXXX().
* CHANGE THOSE HEADER FILES AND REGENERATE THIS FILE.
********************************************************************/
#ifdef __CINT__
#error /home/dariak/Utilities/AculData/AculDataCint.h/C is only for compilation. Abort cint.
#endif
#include <stddef.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#define G__ANSIHEADER
#define G__DICTIONARY
#define G__PRIVATE_GVALUE
#include "G__ci.h"
#include "FastAllocString.h"
extern "C" {
extern void G__cpp_setup_tagtableAculDataCint();
extern void G__cpp_setup_inheritanceAculDataCint();
extern void G__cpp_setup_typetableAculDataCint();
extern void G__cpp_setup_memvarAculDataCint();
extern void G__cpp_setup_globalAculDataCint();
extern void G__cpp_setup_memfuncAculDataCint();
extern void G__cpp_setup_funcAculDataCint();
extern void G__set_cpp_environmentAculDataCint();
}
#include "TObject.h"
#include "TMemberInspector.h"
#include "/home/dariak/Utilities/AculData/AculCalibration.h"
#include "/home/dariak/Utilities/AculData/ConfigDictionary.h"
#include <algorithm>
namespace std { }
using namespace std;
#ifndef G__MEMFUNCBODY
#endif
extern G__linked_taginfo G__AculDataCintLN_TClass;
extern G__linked_taginfo G__AculDataCintLN_TBuffer;
extern G__linked_taginfo G__AculDataCintLN_TMemberInspector;
extern G__linked_taginfo G__AculDataCintLN_TObject;
extern G__linked_taginfo G__AculDataCintLN_string;
extern G__linked_taginfo G__AculDataCintLN_vectorlEROOTcLcLTSchemaHelpercOallocatorlEROOTcLcLTSchemaHelpergRsPgR;
extern G__linked_taginfo G__AculDataCintLN_reverse_iteratorlEvectorlEROOTcLcLTSchemaHelpercOallocatorlEROOTcLcLTSchemaHelpergRsPgRcLcLiteratorgR;
extern G__linked_taginfo G__AculDataCintLN_TObjArray;
extern G__linked_taginfo G__AculDataCintLN_vectorlETVirtualArraymUcOallocatorlETVirtualArraymUgRsPgR;
extern G__linked_taginfo G__AculDataCintLN_reverse_iteratorlEvectorlETVirtualArraymUcOallocatorlETVirtualArraymUgRsPgRcLcLiteratorgR;
extern G__linked_taginfo G__AculDataCintLN_iteratorlEbidirectional_iterator_tagcOTObjectmUcOlongcOconstsPTObjectmUmUcOconstsPTObjectmUaNgR;
extern G__linked_taginfo G__AculDataCintLN_TFile;
extern G__linked_taginfo G__AculDataCintLN_maplEstringcOTObjArraymUcOlesslEstringgRcOallocatorlEpairlEconstsPstringcOTObjArraymUgRsPgRsPgR;
extern G__linked_taginfo G__AculDataCintLN_TH1;
extern G__linked_taginfo G__AculDataCintLN_TVectorTlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TVectorTlEdoublegR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTBaselEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTBaselEdoublegR;
extern G__linked_taginfo G__AculDataCintLN_TCanvas;
extern G__linked_taginfo G__AculDataCintLN_THStack;
extern G__linked_taginfo G__AculDataCintLN_AculCalibration;
extern G__linked_taginfo G__AculDataCintLN_TElementActionTlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TElementPosActionTlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTRow_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTRowlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTDiag_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTColumn_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTFlat_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTSub_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTSparseRow_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTSparseDiag_constlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTColumnlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTDiaglEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTFlatlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTSublEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTSparseRowlEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_TMatrixTSparseDiaglEfloatgR;
extern G__linked_taginfo G__AculDataCintLN_ConfigDictionary;
extern G__linked_taginfo G__AculDataCintLN_maplEstringcOstringcOlesslEstringgRcOallocatorlEpairlEconstsPstringcOstringgRsPgRsPgR;
extern G__linked_taginfo G__AculDataCintLN_maplEstringcOstringcOlesslEstringgRcOallocatorlEpairlEconstsPstringcOstringgRsPgRsPgRcLcLiterator;
/* STUB derived class for protected member access */
TELoss/ELOSS.f90 0 → 100755
View file @ 5016cc63
!======================================================================
!== Energy loss in material
!==
!== File: ELOSS.f90
!== Date: 08/12/2000
!======================================================================
! ELOSS.f90
!
! FUNCTIONS/SUBROUTINES exported from ELOSS.dll:
! ELOSS - subroutine
!
subroutine ELOSS(NEL,ZEL,AEL,WEL,DEN,ZP,AP,NE,ETAB,RE,ZW,AW)
! Expose subroutine ELOSS to users of this DLL
!
!DEC$ ATTRIBUTES DLLEXPORT::ELOSS
IMPLICIT REAL*8(A-H,O-Z)
INTERFACE
SUBROUTINE rde(E, R, R1, R2, IRE)
REAL*8 E, R, R1, R2
INTEGER IRE
END SUBROUTINE rde
END INTERFACE
! Variables
! Input:
INTEGER NEL, NE
REAL*8 ZEL(NEL),AEL(NEL),WEL(NEL), DEN,ZP,AP, ETAB(NE)
! Output:
REAL*8 RE(NE),ZW,AW
! Local
COMMON adj(5),as(5),z(5),f(5),um(5),oz,ro,azm,aj,zion,aion,&
zmed,amed,jc
! Body of ELOSS
zion = ZP
aion = AP
jc = NEL
ro = DEN
!== two parameters with question
ire = 1
oz = 1.
DO i = 1, NEL
as(i) = AEL(i)
z(i) = ZEL(i)
um(i) = WEL(i)
IF(z(i) .LE. 13.) THEN
adj(i) = 12.*z(i) + 7.
ELSE
adj(i) = 9.76*z(i) + 58.8 / z(i)**0.19
END IF
END DO
uma=0.
zam=0.
umed=0.
zme=0.
aj=0.
DO i = 1, jc
uma=uma+um(i)*as(i)
umed=umed+um(i)
zme=zme+um(i)*z(i)
END DO
amed=uma/umed
zmed=zme/umed
DO i = 1, jc
f(i)=um(i)*as(i)/uma
END DO
DO i = 1, jc
zam=zam+f(i)*z(i)/as(i)
END DO
azm=1./zam
aj=0.
DO i = 1, jc
aj=aj+f(i)*z(i)*log(adj(i))/as(i)
END DO
aj=aj*azm
aj=exp(aj)
ZW = zmed
AW = amed
DO j = 1, NE
!== result [mg/cm^2]
CALL rde(ETAB(j), RE(j), rel1, rel, ire)
END DO
end subroutine ELOSS
SUBROUTINE rde(e,range,rel1,rel,ix)
IMPLICIT REAL*8(A-H,O-Z)
INTERFACE
REAL*8 FUNCTION c(x)
REAL*8 x
END FUNCTION c
REAL*8 FUNCTION c1(x)
REAL*8 x
END FUNCTION c1
END INTERFACE
! calculates range and de/dx for compounds
dimension ala(3),ala1(3)
dimension a(3,3),a1(4,4),a2(4,4),b(2,3),cc(5)
common adj(5),as(5),z(5),f(5),um(5),oz,ro,azm,aj,zion,aion,&
zmed,amed,jc
dimension alaj1(4),altau1(4)
data a/ -.75265, .073736, .040556, 2.5398,-.312,.018664, &
-.24598, .11548, -.0099661/
data a1/-8.0155, 0.36916, -1.4307e-02, 3.4718e-03, &
1.8371, -1.452e-02, -3.0142e-02, 2.3603e-03, &
0.045233, -9.5873e-04, 7.1303e-03, -6.8538e-04, &
-5.9898e-03,-5.2315e-04, -3.3802e-04, 3.9405e-05/
data a2/-8.725, 0.8309, -0.13396, 0.012625, &
1.8797, 0.11139, -0.064808, 0.0054043, &
0.74192, -0.528805, 0.1264232, -0.00934142, &
0.752005, -0.5558937, 0.12843119, -0.009306336/
data b/ 0.422377e-06, 3.858019e-09, 0.0304043e-06, -0.1667989e-09, &
-0.00038106e-06, 0.00157955e-09/
if(e.le.0.002)then
range=0.
rel1=0.
rel=0.
return
endif
! this a flag
ibis=-1
en = e/aion
tau = en*1.008
altau = log(tau)
alaj = log(aj)
alaj1(1)= 1.
altau1(1)= 1.
DO kk = 2,4
alaj1(kk)=alaj**(kk-1)
altau1(kk)=altau**(kk-1)
END DO
t = tau/938.256
tt = 2.*t + t*t
beta= sqrt(tt)/(1.+t)
s1=0.
DO i=1,4
DO j=1,4
s1 = s1 + a2(j,i)*alaj1(j)*altau1(i)
END DO
END DO
ala(1) = azm*exp(s1)
s2=0.
DO i=2,4
DO j=1,4
s2 = s2 + (i-1)*a2(j,i)*alaj1(j)*altau1(i-1)
END DO
END DO
ala1(1)=ala(1)*s2/tau
s1=0.
DO i=1,4
DO j=1,4
s1 = s1 + a1(j,i)*alaj1(j)*altau1(i)
END DO
END DO
ala(3) = azm*exp(s1)
s2=0.
DO i=2,4
DO j=1,4
s2 = s2 + (i-1)*a1(j,i)*alaj1(j)*altau1(i-1)
END DO
END DO
ala1(3)=ala(3)*s2/tau
s1=0.
DO i=1,3
DO j=1,3
s1 = s1 + a(j,i)*alaj1(j)*altau1(i)
END DO
END DO
ala(2)=azm*exp(s1)/1000.
s2=0.
DO i=2,3
DO j=1,3
s2 = s2 + (i-1)*a(j,i)*alaj1(j)*altau1(i-1)
END DO
END DO
ala1(2)=ala(2)*s2/tau
25 continue
coefa=3.
coefb=1.
alaa=(ala(1)*(tanh(coefa*(.98 - en))+1.)/2. &
+ ala(2)*(tanh(coefa*(en - .98))+1.)/2.) &
* (tanh(coefb*(8.0 - en))+1.)/2. &
+ ala(3)*(tanh(coefb*(en - 8.))+1.)/2.
alim1=0.
alim2=0.
if(coefa*(en-.98).lt.85)then
alim1=1.008/cosh(coefa*(.98-en))/cosh(coefa*(.98-en))
endif
if(coefb*(en-8.).lt.85)then
alim2=1.008/cosh(coefb*(8.-en))/cosh(coefb*(8.-en))
endif
alaa1=(ala1(1)*(tanh(coefa*(.98-en))+1.)/2.+ &
ala1(2)*(tanh(coefa*(en-.98))+1.)/2.)* &
(tanh(coefb*(8.-en))+1.)/2.+ &
ala1(3)*(tanh(coefb*(en-8.))+1.)/2.+ &
coefa/2.*(tanh(coefb*(8.-en))+1.)/2.* &
(ala(2)*alim1-ala(1)*alim1)+ &
coefb/2.*(ala(3)*alim2- &
(ala(1)*(tanh(coefa*(.98-en))+1.)/2.+ &
ala(2)*(tanh(coefa*(en-.98))+1.)/2.)*alim2)
hi=137.*beta/zion
bz=(31.8+3.86*(aj**.625))*azm*.000001
bz=bz*(zion**2.666667)*c(hi)
bz1=(4.357+.5288*(aj**.625))*azm*.001
bz1=bz1*(zion**1.666667)*c1(hi)
bep=beta*beta
rel1=zion*zion/(alaa1+bz1*((1.-bep)**1.5)/931.141/beta)
rel1=rel1/1000.
range=(alaa+bz)*aion/1.008/zion/zion
range=range*1000.
! Atention!! this version do not work correctly for ix.ne.1
if(ix.ne.1)return
ank=.153574*ro/azm
z23=zion**.666667
abet=beta*125./z23
zef=zion*(1.-exp(-abet))
zef2=zef*zef
om=1022000.*bep/(1.-bep)
cbet=0.
DO k=1,jc
cc(k)=0.
DO i=1,2
DO j=1,3
cc(k)=cc(k)+b(i,j)*((1./bep-1.)**j)*((adj(k)**(i+1)))
END DO
END DO
cbet=cbet+f(k)*cc(k)/as(k)
END DO
cbet=cbet*azm
del1=ank*zef2*(log(om/oz)-bep)/bep
del1=del1/ro
del2=2.*ank*zef2*(log(oz/aj)-cbet)/bep
del2=del2/ro
rel2=rel1-del1
rel3=rel1-rel2+del2
if(del1)8,8,9
8 rel=rel1
goto 12
9 if(del1+del2-rel2)10,10,11
10 rel=rel3
goto 12
11 if(del1.lt.rel1)then
rel=rel2
else
rel=rel1
endif
12 return
END
REAL*8 FUNCTION c(x)
REAL*8 x
IF(x .LE. 0.2) THEN
c = -0.00006 + (0.05252 + 0.12857*x)*x
ELSE IF(x .LE. 2.) THEN
c = -0.00185 + (0.07355 + (0.07172 - 0.02723*x)*x)*x
ELSE IF(x .LE. 3.) THEN
c = -0.0793 + (0.3323 - (0.1234 - 0.0153*x)*x)*x
ELSE
c = 0.22
END IF
END
REAL*8 FUNCTION c1(x)
REAL*8 x
IF(x .LE. 0.2) THEN
c1 = 0.05252+.25714*x
ELSE IF(x .LE. 2.0) THEN
c1 = 0.07355 + (0.14344 - 0.08169*x)*x
ELSE IF(x .LE. 3.0) THEN
c1 = 0.3323 - (0.2468 - 0.0459*x)*x
ELSE
c1 = 0.
END IF
END
TELoss/TELoss.cpp 0 → 100755
View file @ 5016cc63
This diff is collapsed.
TELoss/TELoss.h 0 → 100755
View file @ 5016cc63
#pragma once
#include <TROOT.h>
#include <TObject.h>
#include <TSpline.h>
#include <TArrayD.h>
#include <TMath.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <iostream>
using std::cout;
using std::endl;
class TELoss : public TObject
{
private:
int mel;
int nel;
TArrayD zel, ael, wel;
double den;
double zp,ap;
int ne;
TArrayD etab, rtab, detab; //tables: energy, range, deltaE
// TArrayD vtab;
double zw,aw;
double a, b; //linear interpolation coeficients, y = a*x + b
int bsearch(int ntab, double *xtab, double x);
double aitken(int ntab, double *xtab, double *ytab, double x);
double aitken3(int ntab, double *xtab, double *ytab, double x);
Double_t linear(int ntab, Double_t *xtab, Double_t *ytab, double x);
Double_t linear(int ntab, Double_t *xtab, Double_t *ytab, Int_t i0, double x);
public:
TELoss(void);
TELoss(int _mel, double _den);
virtual ~TELoss(void);
void SetEL(int _mel, double _den);
void SetDen(double _den) {den = _den;}
int AddEL(double _zel, double _ael, double _wel=1.);
void SetZP(double _zp, double _ap);
void SetEtab(int _ne, double e2); //from 0 MeV to e2 MeV, _ne elements
void SetDeltaEtab(double r); //pro kazdou tlostku bude zvlastni tabulka, nebo taky ne
double GetZ() { return zp; };
double GetA() { return ap; };
double GetE(double e0, double r);
double GetE0dE(double de, double r); //de in MeV, r in microns
double GetE0dE(double de); //de in MeV
double GetE0(double e, double r);
double GetEold(double e0, double r);
double GetE0old(double e, double r);
double GetR(double e0, double e);
double GetRold(double e0, double e);
void PrintRE();
void PrintdEE();
void PrintREtoFile();
void PrintdEEtoFile(const char* outfile = "outputdEE.txt");
ClassDef(TELoss, 1)
};
TELoss/TELoss.mk 0 → 100755
View file @ 5016cc63
################################################################################
# TELoss input with some variables
################################################################################
TELOSSLIBS := -lCore -lCint -lMathCore -lMatrix -lHist -lgfortran
# Add inputs and outputs from these tool invocations to the build variables
TELOSSCPP_SRCS += \
$(TELOSS)/TELoss.cpp \
$(TELOSS)/TELossCint.cpp
TELOSSOBJS += \
$(TELOSS)/ELOSS.o \
$(TELOSS)/TELoss.o \
$(TELOSS)/TELossCint.o
TELOSSCPP_DEPS += \
$(TELOSS)/TELoss.d \
$(TELOSS)/TELossCint.d
TELoss/linkdef.h 0 → 100755
View file @ 5016cc63
#ifdef __CINT__
#pragma link off all globals;
#pragma link off all classes;
#pragma link off all functions;
#pragma link C++ class TELoss;
#endif
\ No newline at end of file
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment