LMU/musrsim
5
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

27.1.2011 Kamil Sedlak

Browse Source 
1) correction of volume that is assigned to optical photons
   (odet_ID) - now it is the volume from postStepPoint instead of
   preStepPoint
2) new variable added (nOptPhot) - number of optical photons
   generated in the event
This commit is contained in:
sedlak
2011-01-27 14:36:03 +00:00
parent 543a0e667b
commit 7a8ea16943
7 changed files with 137 additions and 92 deletions
Download Patch File Download Diff File Expand all files Collapse all files

BIN
doc/musrSim.pdf
View File

Binary file not shown.

49
doc/musrSim.tex
View File

@ -688,8 +688,16 @@ Three special volumes ``Target, M0, M1 and M2''.
Normaly the simulation of ``detected signal'' stops at the level of the deposited energy in
a sensitive volume (e.g.\ in a scintilator tile). However, in some special cases, one would
like to know how the light is propagated through the scintillators. In such case simulation
of optical photons is possible. It will, however, significantly (in some cases dramatically)
reduce the speed of the simulation. The output of the optical photon simulation is stored
of optical photons is possible. Note that the optical photon in Geant are treaded as a class
of particles distinct from higher energy gamma particles -- and there is no smooth transition
between the two. Some additional material properties
of an active detector and of the detector surface have to be defined for optical photons.
We strongly recommend the users willing to simulate optical photons to read the
chapter ``Optical Photon Processes'' in~\cite{geantUserManual} to understand the rest of
this chapter.
The simulation of optical photons will significantly (in some cases dramatically)
reduce the speed of the program. The output of the optical photon simulation is stored
in variables starting with the string ``odet\_''. The user has to specify several parameters
in order to simulate optical photons:
%
@ -700,21 +708,46 @@ in order to simulate optical photons:
will be ignored internally in musrSim, and the user therefore does not have to
comment out other parameters connected with optical photons in the macro file.
\item{\bf /musr/command materialPropertiesTable \emph{MPT\_name} \emph{property} \emph{n} \emph{val(1)} ... \emph{val(n)} \emph{E(1)} ... \emph{E(n)} \\
\item{\bf /musr/command materialPropertiesTable \emph{MPT\_name} \emph{property} \emph{n}
\emph{E(1)} ... \emph{E(n)} \emph{val(1)} ... \emph{val(n)} } \\
Defines some optical \emph{property} of a given material (e.g.\ absorption lenght,
refractive index, scintillation yield, ...} to a material property table.
refractive index, scintillation yield, ...) to a material property table.
\emph{MPT\_name} stands for the material property table name. The table has
\emph{n} different values \emph{val(1)} ... \emph{val(n)} for the same number
of optical photon energies \emph{E(1)} ... \emph{E(n)} expressed in MeV.
If \emph{n}=0, the \emph{property} is called ``constant property''.
Possible \emph{property} keywords are: ABSLENGTH, RINDEX, FASTCOMPONENT, SLOWCOMPONENT,
SCINTILLATIONYIELD, and constant \emph{properties} are SCINTILLATIONYIELD, RESOLUTIONSCALE,
FASTTIMECONSTANT, SLOWTIMECONSTANT, and YIELDRATIO.
Some of the possible \emph{property} keywords are: ABSLENGTH, RINDEX, FASTCOMPONENT, SLOWCOMPONENT,
SCINTILLATIONYIELD, and some of the constant \emph{properties} are SCINTILLATIONYIELD, RESOLUTIONSCALE,
FASTTIMECONSTANT, SLOWTIMECONSTANT, and YIELDRATIO. See other \emph{property} keywords
in chapter ``Optical Photon Processes'' in~\cite{geantUserManual}.\\
\begin{description}
\item{\bf SCINTILLATIONYIELD} ... nr. of photons per 1\,MeV of deposited energy.
\item{\bf RESOLUTIONSCALE} ... intrinsic resolution -- normally 1, larger than
1 for crystals with impurities, smaller than 1 when Fano factor plays a role.
\item{\bf YIELDRATIO} ... relative strength of the fast component as a fraction
of total scintillation yeald.
\end{description}
\item {\bf /musr/command setMaterialPropertiesTable \emph{MPT\_name} \emph{material\_name}} \\
Assigns a material property table defined by ``/musr/command materialPropertiesTable''
to a given material.
\item {\bf /musr/command opticalSurface \emph{surface\_name} \emph{phys\_volName1} \emph{phys\_volName2} \emph{surfaceType} \emph{surfaceFinish} \emph{surfaceModel} \emph{MPT\_name}}\\
Defines ``border surface'' (G4LogicalBorderSurface, see subsection ``Boundary Process''
of chapter ``Optical Photon Processes'' in ~\cite{geantUserManual}).
\begin{description}
\item{\bf surface\_name} name of the newly defined surface.
\item{\bf phys\_volName1, phys\_volName2} names of the physics volume in question (it is an ordered pair!).
\item{\bf surfaceType} one of the following: dielectric\_dielectric, dielectric\_metal, dielectric\_LUT, firsov, x\_ray.
\item{\bf surfaceFinish} surface finish properties, such as: polished, ground, etchedair, ...
\item{\bf surfaceModel} one of the following: glisur, unified, LUT.
\item{\bf MPT\_name} optional ``material properties table name'' which should be assigned to the surface.
E.g.\ REFLECTIVITY or EFFICIENCY for a given surface may be assigned this way.
\end{description}
One has to assign a non-zero EFFICIENCY and a REFLECTIVITY smaller than 1 to a boundary surface
between the scintillator and sensitive device (e.g.\ an APD).
\end{description}
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
@ -1097,6 +1130,8 @@ in~\cite{Aktas:2004px}.
Nucl. Inst. and Meth. in Phys. Res. A 389 (1997) 81.%-86.
See also http://root.cern.ch/.
\bibitem{geantUserManual} Geant4 User Manual
\bibitem{turtle}
K.L.~Brown, Ch.~Iselin, D.C.~Carey, {\it``Decay Turtle''}, CERN 74-2 (1974). \\
U.~Rohrer, {\it ``Compendium of Decay Turtle Enhancements''},

25
include/musrRootOutput.hh
View File

@ -100,17 +100,18 @@ class musrRootOutput {
G4cout<<" numberOfGeneratedEvents = "<<GeantParametersD[7]<<G4endl;
}
void SetPolInTarget(G4ThreeVector pol) {muTargetPolX=pol.x(); muTargetPolY=pol.y(); muTargetPolZ=pol.z();};
void SetTimeInTarget(G4double time) {muTargetTime = time/microsecond;};
void SetMomentumInTarget(G4ThreeVector mom) {muTargetMomX=(mom.x())/MeV; muTargetMomY=(mom.y())/MeV; muTargetMomZ=(mom.z())/MeV;};
void SetPolInM0(G4ThreeVector pol) {muM0PolX=pol.x(); muM0PolY=pol.y(); muM0PolZ=pol.z();};
void SetTimeInM0(G4double time) {muM0Time = time/microsecond;};
void SetPolInM1(G4ThreeVector pol) {muM1PolX=pol.x(); muM1PolY=pol.y(); muM1PolZ=pol.z();};
void SetTimeInM1(G4double time) {muM1Time = time/microsecond;};
void SetPolInM2(G4ThreeVector pol) {muM2PolX=pol.x(); muM2PolY=pol.y(); muM2PolZ=pol.z();};
void SetTimeInM2(G4double time) {muM2Time = time/microsecond;};
void SetInitialPositronMomentum(G4ThreeVector mom) {posIniMomx=mom.x(); posIniMomy=mom.y(); posIniMomz=mom.z();};
void SetDecayTime(G4double time) {muDecayTime=time/microsecond;};
void SetPolInTarget(G4ThreeVector pol) {muTargetPolX=pol.x(); muTargetPolY=pol.y(); muTargetPolZ=pol.z();}
void SetTimeInTarget(G4double time) {muTargetTime = time/microsecond;}
void SetMomentumInTarget(G4ThreeVector mom) {muTargetMomX=(mom.x())/MeV; muTargetMomY=(mom.y())/MeV; muTargetMomZ=(mom.z())/MeV;}
void SetPolInM0(G4ThreeVector pol) {muM0PolX=pol.x(); muM0PolY=pol.y(); muM0PolZ=pol.z();}
void SetTimeInM0(G4double time) {muM0Time = time/microsecond;}
void SetPolInM1(G4ThreeVector pol) {muM1PolX=pol.x(); muM1PolY=pol.y(); muM1PolZ=pol.z();}
void SetTimeInM1(G4double time) {muM1Time = time/microsecond;}
void SetPolInM2(G4ThreeVector pol) {muM2PolX=pol.x(); muM2PolY=pol.y(); muM2PolZ=pol.z();}
void SetTimeInM2(G4double time) {muM2Time = time/microsecond;}
void SetInitialPositronMomentum(G4ThreeVector mom) {posIniMomx=mom.x(); posIniMomy=mom.y(); posIniMomz=mom.z();}
void SetNOptPhot(G4int value) {nOptPhot=value;}
void SetDecayTime(G4double time) {muDecayTime=time/microsecond;}
void SetNrFieldNomVal(G4int n) {nFieldNomVal = n;}
void SetFieldNomVal(G4int i, G4double value);
G4int GetNrOfVolumes() {return det_nMax;}
@ -174,6 +175,7 @@ class musrRootOutput {
static G4bool store_posIniMomX;
static G4bool store_posIniMomY;
static G4bool store_posIniMomZ;
static G4bool store_nOptPhot;
static G4bool store_det_ID;
static G4bool store_det_edep;
static G4bool store_det_edep_el;
@ -259,6 +261,7 @@ class musrRootOutput {
Double_t muDecayPosX, muDecayPosY, muDecayPosZ;
Double_t muDecayTime;
Double_t posIniMomx, posIniMomy, posIniMomz;
Int_t nOptPhot;
public:
static const Int_t maxNFieldnNominalValues=30;

6
musrSim.cc
View File

@ -3,7 +3,7 @@
#include "musrPrimaryGeneratorAction.hh"
#include "musrRunAction.hh"
#include "musrEventAction.hh"
//#include "StackingAction.hh"
#include "musrStackingAction.hh"
#include "musrSteppingAction.hh"
#include "musrSteppingVerbose.hh"
@ -97,7 +97,8 @@ int main(int argc,char** argv) {
runManager->SetUserAction(new musrPrimaryGeneratorAction(musrdetector));
runManager->SetUserAction(new musrRunAction);
runManager->SetUserAction(new musrEventAction);
// runManager->SetUserAction(new StackingAction);
// Initiate musrStackingAction only if optical photons are required (otherwise not needed)
if (musrParameters::boolG4OpticalPhotons) runManager->SetUserAction(new musrStackingAction);
runManager->SetUserAction(new musrSteppingAction);
//Initialize G4 kernel
@ -136,6 +137,7 @@ int main(int argc,char** argv) {
#else
session = new G4UIterminal();
#endif
G4cout<<"Go to idle state now:"<<G4endl;
session->SessionStart();
delete session;
}

141
src/musrDetectorConstruction.cc
View File

@ -153,8 +153,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
// Define rotation matrix, which might be later on used for some volumes
if (strcmp(tmpString1,"rotation")==0){
char matrixName[100]="Unset";
float pp1, pp2, pp3, pp4=0;
sscanf(&line[0],"%*s %*s %s %g %g %g %g",matrixName,&pp1,&pp2,&pp3,&pp4);
double pp1, pp2, pp3, pp4=0;
sscanf(&line[0],"%*s %*s %s %lf %lf %lf %lf",matrixName,&pp1,&pp2,&pp3,&pp4);
if (pp4==0) {
G4RotationMatrix* pRotMatrix = new G4RotationMatrix(pp1*deg,pp2*deg,pp3*deg); // pp1=phi, pp2=theta, pp3=psi
pointerToRotationMatrix[matrixName]=pRotMatrix;
@ -167,8 +167,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
else if (strcmp(tmpString1,"construct")==0){
float x1=0,x2=0,x3=0,x4=0,x5=0,x6=0,x7=0,x8=0,x9=0,x10=0,x11,x12;
float posx,posy,posz;
double x1=0,x2=0,x3=0,x4=0,x5=0,x6=0,x7=0,x8=0,x9=0,x10=0,x11,x12;
double posx,posy,posz;
char name[100];
char mothersName[100];
char material[100];
@ -181,49 +181,51 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
char actualFieldName[100]="nofield";
if (strcmp(tmpString2,"tubs")==0){
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
solid = new G4Tubs(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg);
}
else if (strcmp(tmpString2,"cons")==0){
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
solid = new G4Cons(solidName,x1*mm,x2*mm,x3*mm,x4*mm,x5*mm,x6*deg,x7*deg);
}
else if (strcmp(tmpString2,"box")==0){
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %s %g %g %g %s %s",
// sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %s %lf %lf %lf %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
// G4cout<<"DEBUG KAMIL: xxx1="<<xxx1<<", xxx2="<<xxx2<<", xxx3="<<xxx3<<G4endl;
solid = new G4Box(solidName,x1*mm,x2*mm,x3*mm);
}
else if ((strcmp(tmpString2,"trd")==0)||(strcmp(tmpString2,"trd90y")==0)) {
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
solid = new G4Trd(solidName,x1*mm,x2*mm,x3*mm,x4*mm,x5*mm);
}
else if (strcmp(tmpString2,"sphere")==0){
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
solid = new G4Sphere(solidName,x1*mm,x2*mm,x3*deg,x4*deg,x5*deg,x6*deg);
}
else if (strcmp(tmpString2,"para")==0){ // NOT YET TESTED
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
solid = new G4Para(solidName,x1*mm,x2*mm,x3*mm,x4*deg,x5*deg,x6*deg);
}
else if (strcmp(tmpString2,"cylpart")==0){ // Volume introduced by Pavel Bakule on 12 May 2009
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -236,7 +238,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
// Create a U-profile geometry. x1, x2, x3 define the outer dimensions of the U-profile (as a box),
// x4 is the wall thickness of the U-profile. The centre of the U-profile
// is in the centre of the box defined by x1,x2,x3.
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -252,7 +254,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
else if (strcmp(tmpString2,"alcSupportPlate")==0){
// Create an ALC holder geometry: x1 half-width of the holder (as a box),
// x2 half-height of the spacer
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -289,7 +291,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
// | |
// ----------
// First 5 parameters as for the outer tube, the 6th parameter is the depth of the hole.
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,&x6,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -316,7 +318,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
// | |
// ----------
// First 5 parameters as for the outer tube, the 6th parameter is the depth of the hole.
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -330,7 +332,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(tmpString2,"tubsbox")==0){
// Create a tube, from which center a box is cut out. x1=box half-width; x2,x3,x4,x5 define the tube.
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -343,9 +345,9 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
// Create a box, from which a box is cut out. x1, x2, x3 = box half-widths;
// x4,x5,x6,x7 define the tube, x8, x9 and x10 are the distances between the tube and box centres.
// x11, x12 are the half-withs of the ractangular opening in the collimator
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %g %g %g %g %g %g %g %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %lf %s",
name,&x1,&x2,&x3,&x4,&x5,&x6,&x7,&x8,&x9,&x10,&x11,&x12,material);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %g %g %g %s %s %s %d %s",&posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*g %*s %lf %lf %lf %s %s %s %d %s",&posx,&posy,&posz,mothersName,rotMatrix,sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
G4Box* solidDetBox = new G4Box("SolidDetBox",x1*mm,x2*mm,x3*mm);
G4Box* solidHole = new G4Box("SolidDetBox",x11*mm,x2*mm+0.1,x12*mm);
@ -358,7 +360,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
else if (strcmp(tmpString2,"tubsboxsegm")==0){
// Create a volume that looks like an intersection of tube and box.
char orientation[100];
sscanf(&line[0],"%*s %*s %*s %s %s %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,orientation,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d %s",sensitiveDet,&volumeID,actualFieldName);
solidName+=name;
@ -381,7 +383,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
else if (strcmp(tmpString2,"GPDsampleHolderA")==0){
// First part of the GPD sample holder, where posx, posy, posz = centre of the whole (long) tube
// (=111.25mm below the centre of the holes)
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,&x4,&x5,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID);
solidName+=name;
@ -405,7 +407,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(tmpString2,"GPDmHolder")==0){
// Light guide that holds m0 in its position
sscanf(&line[0],"%*s %*s %*s %s %g %g %g %s %g %g %g %s %s",
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf %s %lf %lf %lf %s %s",
name,&x1,&x2,&x3,material,&posx,&posy,&posz,mothersName,rotMatrix);
sscanf(&line[0],"%*s %*s %*s %*s %*g %*g %*g %*s %*g %*g %*g %*s %*s %s %d",sensitiveDet,&volumeID);
solidName+=name;
@ -641,18 +643,21 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
G4OpticalSurfaceFinish OpticalFinish = OpticalFinishMap[finish];
if ((OpticalType==0)&&(strcmp(type,"dielectric_metal")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical type \""<<type<<"\" not found!"<<G4endl;
G4cout<<" (Maybe it exists in Geant4 but was not added to musrSim ==> update musrDetectorConstruction.cc file!)"<< G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
exit(1);
}
if ((OpticalModel==0)&&(strcmp(model,"glisur")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical surface model \""<<model<<"\" not found!"<<G4endl;
G4cout<<" (Maybe it exists in Geant4 but was not added to musrSim ==> update musrDetectorConstruction.cc file!)"<< G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
exit(1);
}
if ((OpticalFinish==0)&&(strcmp(finish,"polished")!=0)) {
G4cout<<"ERROR! musrDetectorConstruction::Construct(): Optical surface finish \""<<finish<<"\" not found!"<<G4endl;
G4cout<<" (Maybe it exists in Geant4 but was not added to musrSim ==> update musrDetectorConstruction.cc file!)"<< G4endl;
G4cout << " ==> S T O P F O R C E D"<<G4endl;
G4cout<<" "<<line;
exit(1);
@ -663,7 +668,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
optSurfTMP->SetFinish(OpticalFinish);
optSurfTMP->SetModel(OpticalModel);
// Assign the "Material properties table" if required by the user:
G4cout<<"materialPropertiesTableName="<<materialPropertiesTableName<<G4endl;
// G4cout<<"materialPropertiesTableName="<<materialPropertiesTableName<<G4endl;
if (strcmp(materialPropertiesTableName,"Undefined")!=0) {
G4MaterialPropertiesTable* MPT_tmp=NULL;
itMPT = materialPropertiesTableMap.find(materialPropertiesTableName);
@ -679,8 +684,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
MPT_tmp = itMPT->second;
}
optSurfTMP->SetMaterialPropertiesTable(MPT_tmp);
G4cout<<optSurfTMP<<G4endl;
optSurfTMP->GetMaterialPropertiesTable()->DumpTable();
// optSurfTMP->GetMaterialPropertiesTable()->DumpTable();
}
G4cout<<"Optical surface \""<<optSurfaceName<<"\" created. OpticalType="<<OpticalType
<<" OpticalFinish="<<OpticalFinish<<" OpticalModel="<<OpticalModel<<G4endl;
@ -695,7 +699,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
musrErrorMessage::GetInstance()->musrError(FATAL,eMessage,false);
}
char varName[100];
float fVarValue;
double fVarValue;
int iVarValue;
sscanf(&line[0],"%*s %*s %s",varName);
if (strcmp(varName,"minNrOfDetectedPhotons")==0) {
@ -703,12 +707,12 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
myMusrScintSD -> Set_OPSA_minNrOfDetectedPhotons(iVarValue);
}
else if (strcmp(varName,"signalSeparationTime")==0) {
sscanf(&line[0],"%*s %*s %*s %g",&fVarValue);
sscanf(&line[0],"%*s %*s %*s %lf",&fVarValue);
myMusrScintSD -> Set_OPSA_SignalSeparationTime(fVarValue*nanosecond);
}
else if (strcmp(varName,"photonFractions")==0) {
float a, b, c, d, e;
sscanf(&line[0],"%*s %*s %*s %g %g %g %g %g",&a, &b, &c, &d, &e);
double a, b, c, d, e;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %lf %lf",&a, &b, &c, &d, &e);
myMusrScintSD -> Set_OPSA_frac(a,b,c,d,e);
}
else if (strcmp(varName,"eventsForOPSAhistos")==0) {
@ -718,8 +722,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(varName,"OPSAhist")==0) {
int nBins;
float min, max;
sscanf(&line[0],"%*s %*s %*s %d %g %g",&nBins,&min,&max);
double min, max;
sscanf(&line[0],"%*s %*s %*s %d %lf %lf",&nBins,&min,&max);
myMusrScintSD -> SetOPSAhistoBinning(nBins,min*nanosecond,max*nanosecond);
}
}
@ -734,18 +738,18 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
(void)F04GlobalField::getObject();
char typeOfField[100]="Unset";
float pp1=0; float pp2=0; float pp3=0;
sscanf(&line[0],"%*s %*s %*s %g %g %g %s",&pp1,&pp2,&pp3,typeOfField);
double pp1=0, pp2=0, pp3=0;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %s",&pp1,&pp2,&pp3,typeOfField);
G4ThreeVector position = G4ThreeVector(pp1,pp2,pp3);
float fieldValue=0.000000001;
float fieldValueFinal=0;
double fieldValue=0.000000001;
double fieldValueFinal=0;
int fieldNrOfSteps=0;
// if (strcmp(tmpString2,"magnetic")==0){
if (strcmp(typeOfField,"fromfile")==0) {
char fieldInputFileName[100];
char fieldTableType[100];
char logicalVolumeName[100];
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %s %s %g %g %d",fieldTableType,fieldInputFileName,logicalVolumeName,
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %s %s %lf %lf %d",fieldTableType,fieldInputFileName,logicalVolumeName,
&fieldValue,&fieldValueFinal,&fieldNrOfSteps);
// Find out the logical volume, to which the field will be placed:
G4LogicalVolume* logVol = FindLogicalVolume(logicalVolumeName);
@ -770,19 +774,17 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(typeOfField,"uniform")==0) {
float fieldValue[6];
double fieldValue[6];
char logicalVolumeName[100];
float positionX, positionY, positionZ;
double positionX, positionY, positionZ;
G4double half_x = pp1; // halfwith of the box with the field
G4double half_y = pp2; // halfhight of the box with the field
G4double half_z = pp3; // halflenght of the box with the field
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %g %g %g %s %g %g %g %g %g %g",
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %lf %lf %lf %s %lf %lf %lf %lf %lf %lf",
&positionX, &positionY, &positionZ, logicalVolumeName,
&fieldValue[0],&fieldValue[1],&fieldValue[2],&fieldValue[3],&fieldValue[4],&fieldValue[5]);
G4ThreeVector position = G4ThreeVector(positionX,positionY,positionZ);
//old sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %s %g %g %g %g %g %g", ////// %g %d",
//old logicalVolumeName, &fieldValue[0],&fieldValue[1],&fieldValue[2],&fieldValue[3],&fieldValue[4],&fieldValue[5]);
G4LogicalVolume* logVol = FindLogicalVolume(logicalVolumeName);
if (logVol==NULL) {
sprintf(eMessage,"musrDetectorConstruction.cc::Construct(): GLOBAL FIELD (uniform): Logical volume \"%s\" not found.", logicalVolumeName);
@ -797,10 +799,10 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(typeOfField,"quadrupole")==0) {
float halfLength, fieldRadius, gradientValue, gradientValueFinal, fringeFactor;
double halfLength, fieldRadius, gradientValue, gradientValueFinal, fringeFactor;
int gradientNrOfSteps = 0;
char logicalVolumeName[100];
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %g %g %g %s %g %g %d",&halfLength,&fieldRadius,&fringeFactor,logicalVolumeName,
sscanf(&line[0],"%*s %*s %*s %*g %*g %*g %*s %lf %lf %lf %s %lf %lf %d",&halfLength,&fieldRadius,&fringeFactor,logicalVolumeName,
&gradientValue,&gradientValueFinal,&gradientNrOfSteps);
G4LogicalVolume* logVol = FindLogicalVolume(logicalVolumeName);
if (logVol==NULL) {
@ -831,8 +833,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else {
char parameterName[100];
float parameterValue;
sscanf(&line[0],"%*s %*s %*s %s %g",parameterName,&parameterValue);
double parameterValue;
sscanf(&line[0],"%*s %*s %*s %s %lf",parameterName,&parameterValue);
if (strcmp(parameterName,"SetDeltaIntersection")==0){ fieldMgr->SetDeltaIntersection(parameterValue*mm); }
else if (strcmp(parameterName,"SetDeltaOneStep")==0){ fieldMgr->SetDeltaOneStep(parameterValue*mm); }
else if (strcmp(parameterName,"SetMinimumEpsilonStep")==0){ fieldMgr->SetMinimumEpsilonStep(parameterValue); }
@ -872,8 +874,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if ((strcmp(tmpString2,"printFieldValueAtPoint")==0)||(strcmp(tmpString2,"printFieldDerivativeAtPoint")==0)) { // Print the fieldvalue at the given point
float p0, p1, p2;
sscanf(&line[0],"%*s %*s %*s %g %g %g",&p0,&p1,&p2);
double p0, p1, p2;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf",&p0,&p1,&p2);
if (F04GlobalField::Exists()) {
F04GlobalField* myGlobalField = F04GlobalField::getObject();
if (myGlobalField) {
@ -901,13 +903,13 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
<< tmpString2 <<"\" not found!"<<G4endl<<"S T O P F O R C E D"<<G4endl;
ReportGeometryProblem(line);
}
// float maxStep, ;
float ustepMax = -1;
float utrakMax = -1;
float utimeMax = -1;
float uekinMin = -1;
float urangMin = -1;
sscanf(&line[0],"%*s %*s %*s %g %g %g %g %g", &ustepMax, &utrakMax, &utimeMax, &uekinMin, &urangMin);
double ustepMax = -1;
double utrakMax = -1;
double utimeMax = -1;
double uekinMin = -1;
double urangMin = -1;
sscanf(&line[0],"%*s %*s %*s %lf %lf %lf %lf %lf", &ustepMax, &utrakMax, &utimeMax, &uekinMin, &urangMin);
G4UserLimits* myUserLimits = new G4UserLimits();
G4cout<<"musrDetectorConstruction.cc: G4UserLimits in "<<tmpString2<<": ";
if (ustepMax>0) {myUserLimits->SetMaxAllowedStep(ustepMax*mm); G4cout<<"ustepMax = "<<ustepMax<<" mm, ";}
@ -969,14 +971,14 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
else if (strcmp(tmpString1,"signalSeparationTime")==0){
float timeSeparation;
sscanf(&line[0],"%*s %*s %g",&timeSeparation);
double timeSeparation;
sscanf(&line[0],"%*s %*s %lf",&timeSeparation);
musrParameters::signalSeparationTime = timeSeparation*nanosecond;
}
else if (strcmp(tmpString1,"maximumRunTimeAllowed")==0){ // in seconds
float timeMax;
sscanf(&line[0],"%*s %*s %g",&timeMax);
double timeMax;
sscanf(&line[0],"%*s %*s %lf",&timeMax);
musrEventAction::maximumRunTimeAllowed = timeMax;
}
@ -1031,8 +1033,8 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
}
else if (strcmp(tmpString2,"setProductionCut")==0) {
char charRegionName[100];
float fGammaCut=0, fElectronCut=0, fPositronCut=0;
sscanf(&line[0],"%*s %*s %*s %s %g %g %g",charRegionName,&fGammaCut,&fElectronCut,&fPositronCut);
double fGammaCut=0, fElectronCut=0, fPositronCut=0;
sscanf(&line[0],"%*s %*s %*s %s %lf %lf %lf",charRegionName,&fGammaCut,&fElectronCut,&fPositronCut);
G4String regionName = charRegionName;
G4Region* myRegion = G4RegionStore::GetInstance()->GetRegion(regionName,false);
if( myRegion == NULL ) { // G4Region does not exist
@ -1110,6 +1112,7 @@ G4VPhysicalVolume* musrDetectorConstruction::Construct() {
if (strcmp(tmpString2,"posIniMomX")==0) {musrRootOutput::store_posIniMomX = false;}
if (strcmp(tmpString2,"posIniMomY")==0) {musrRootOutput::store_posIniMomY = false;}
if (strcmp(tmpString2,"posIniMomZ")==0) {musrRootOutput::store_posIniMomZ = false;}
if (strcmp(tmpString2,"nOptPhot")==0) {musrRootOutput::store_nOptPhot = false;}
if (strcmp(tmpString2,"fieldNomVal")==0) {musrRootOutput::store_fieldNomVal = false;}
if (strcmp(tmpString2,"det_ID")==0) {musrRootOutput::store_det_ID = false;}
if (strcmp(tmpString2,"det_edep")==0) {musrRootOutput::store_det_edep = false;}
@ -1336,30 +1339,28 @@ void musrDetectorConstruction::DefineMaterials()
char propertyName[100];
int nEntries;
sscanf(&line[0],"%*s %*s %*s %s %d",propertyName,&nEntries);
std::cout<<" Optical Material Def: MPT_tmp="<<MPT_tmp<<", materialPropertiesTableName="<<materialPropertiesTableName
<<", propertyName="<<propertyName<<", nEntries="<<nEntries<<std::endl;
// G4cout<<" Optical Material Def: MPT_tmp="<<MPT_tmp<<", materialPropertiesTableName="<<materialPropertiesTableName
// <<", propertyName="<<propertyName<<", nEntries="<<nEntries<<G4endl;
if (nEntries==0) { // AddConstProperty
float value;
sscanf(&line[0],"%*s %*s %*s %*s %*d %g",&value);
double value;
sscanf(&line[0],"%*s %*s %*s %*s %*d %lf",&value);
MPT_tmp->AddConstProperty(propertyName,value);
}
else { // AddProperty
char* pch = strstr(line,propertyName)+strlen(propertyName);
float value;
double value;
G4double photonEnergyArray[100];
G4double valueArray[100];
char dummy[100];
sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
for (int i=0; i<nEntries; i++) {
sscanf(pch,"%g",&value);
// G4cout<<" DDDDD var1="<<value<<" &pch="<<&pch<<G4endl;
sscanf(pch,"%lf",&value);
photonEnergyArray[i]=value;
sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
}
for (int i=0; i<nEntries; i++) {
sscanf(pch,"%g",&value);
// G4cout<<" DDDDD var2="<<value<<" &pch="<<&pch<<G4endl;
sscanf(pch,"%lf",&value);
valueArray[i]=value;
sscanf(pch,"%s",dummy); char* pch2=strstr(pch,dummy)+strlen(dummy); pch = pch2;
}
@ -1388,8 +1389,8 @@ void musrDetectorConstruction::DefineMaterials()
MPT_tmp = itMPT->second;
}
myMaterial->SetMaterialPropertiesTable(MPT_tmp);
G4cout<<myMaterial<<G4endl;
myMaterial->GetMaterialPropertiesTable()->DumpTable();
// G4cout<<myMaterial<<G4endl;
// myMaterial->GetMaterialPropertiesTable()->DumpTable();
}
}
}

3
src/musrRootOutput.cc
View File

@ -110,6 +110,7 @@ G4bool musrRootOutput::store_muM2PolZ = false;
G4bool musrRootOutput::store_posIniMomX = true;
G4bool musrRootOutput::store_posIniMomY = true;
G4bool musrRootOutput::store_posIniMomZ = true;
G4bool musrRootOutput::store_nOptPhot = true;
G4bool musrRootOutput::store_det_ID = true;
G4bool musrRootOutput::store_det_edep = true;
G4bool musrRootOutput::store_det_edep_el = true;
@ -214,6 +215,7 @@ void musrRootOutput::BeginOfRunAction() {
if (store_posIniMomX) {rootTree->Branch("posIniMomX",&posIniMomx,"posIniMomX/D");}
if (store_posIniMomY) {rootTree->Branch("posIniMomY",&posIniMomy,"posIniMomY/D");}
if (store_posIniMomZ) {rootTree->Branch("posIniMomZ",&posIniMomz,"posIniMomZ/D");}
if (store_nOptPhot) {rootTree->Branch("nOptPhot",&nOptPhot,"nOptPhot/I");}
// if (store_globalTime) {rootTree->Branch("globalTime",&globalTime,"globalTime/D");}
// if (store_fieldValue) {rootTree->Branch("fieldValue",&fieldValue,"fieldValue/D");}
if (store_fieldNomVal) {
@ -368,6 +370,7 @@ void musrRootOutput::ClearAllRootVariables() {
muDecayPosX=-1000;muDecayPosY=-1000;muDecayPosZ=-1000;
muDecayTime=-1000;
posIniMomx=-1000;posIniMomy=-1000;posIniMomz=-1000;
nOptPhot=0;
BxIntegral = -1000; ByIntegral = -1000; BzIntegral = -1000;
BzIntegral1 = -1000; BzIntegral2 = -1000; BzIntegral3 = -1000;
det_n=0;

5
src/musrScintSD.cc
View File

@ -212,9 +212,10 @@ void musrScintSD::ProcessOpticalPhoton(G4Step* aStep) {
//Prior to Geant4.6.0-p1 this would not have been enough to check
if(aStep->GetPostStepPoint()->GetStepStatus()==fGeomBoundary){
// G4cout<<" boundaryStatus="<<boundaryStatus<<" ";
G4String actualVolume = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
// G4String actualVolume = aStep->GetTrack()->GetVolume()->GetLogicalVolume()->GetName();
G4String actualVolume = aStep->GetPostStepPoint()->GetPhysicalVolume()->GetLogicalVolume()->GetName();
Int_t detID = myRootOutput->ConvertVolumeToID(actualVolume);
// G4cout<<" detID ="<<detID<<" actualVolume="<<actualVolume;
// G4cout<<" detID ="<<detID<<" actualVolume="<<actualVolume<<G4endl;
optHitMapType::iterator iter = optHitMap.find(detID);
if (iter==optHitMap.end()) { // optHitDetectorMapType does not exist ==> create it
optHitDetectorMapType* optHitDetectorMapTmp = new optHitDetectorMapType;