diff --git a/doc/musrSim.pdf b/doc/musrSim.pdf index b31e50d..572dda6 100644 Binary files a/doc/musrSim.pdf and b/doc/musrSim.pdf differ diff --git a/doc/musrSim.tex b/doc/musrSim.tex index 4a87d50..9a1bf4a 100644 --- a/doc/musrSim.tex +++ b/doc/musrSim.tex @@ -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''}, diff --git a/include/musrRootOutput.hh b/include/musrRootOutput.hh index af7a7cc..6857a33 100644 --- a/include/musrRootOutput.hh +++ b/include/musrRootOutput.hh @@ -100,17 +100,18 @@ class musrRootOutput { G4cout<<" numberOfGeneratedEvents = "<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:"<SessionStart(); delete session; } diff --git a/src/musrDetectorConstruction.cc b/src/musrDetectorConstruction.cc index 49d671d..a8007d1 100644 --- a/src/musrDetectorConstruction.cc +++ b/src/musrDetectorConstruction.cc @@ -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="< update musrDetectorConstruction.cc file!)"<< G4endl; G4cout << " ==> S T O P F O R C E D"< update musrDetectorConstruction.cc file!)"<< G4endl; G4cout << " ==> S T O P F O R C E D"< update musrDetectorConstruction.cc file!)"<< G4endl; G4cout << " ==> S T O P F O R C E D"<SetFinish(OpticalFinish); optSurfTMP->SetModel(OpticalModel); // Assign the "Material properties table" if required by the user: - G4cout<<"materialPropertiesTableName="<second; } optSurfTMP->SetMaterialPropertiesTable(MPT_tmp); - G4cout<GetMaterialPropertiesTable()->DumpTable(); + // optSurfTMP->GetMaterialPropertiesTable()->DumpTable(); } G4cout<<"Optical surface \""< 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,¶meterValue); + double parameterValue; + sscanf(&line[0],"%*s %*s %*s %s %lf",parameterName,¶meterValue); 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!"<0) {myUserLimits->SetMaxAllowedStep(ustepMax*mm); G4cout<<"ustepMax = "<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="<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; isecond; } myMaterial->SetMaterialPropertiesTable(MPT_tmp); - G4cout<GetMaterialPropertiesTable()->DumpTable(); + // G4cout<GetMaterialPropertiesTable()->DumpTable(); } } } diff --git a/src/musrRootOutput.cc b/src/musrRootOutput.cc index 4012aa6..79adde8 100644 --- a/src/musrRootOutput.cc +++ b/src/musrRootOutput.cc @@ -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; diff --git a/src/musrScintSD.cc b/src/musrScintSD.cc index 7fe042b..ec14fd3 100644 --- a/src/musrScintSD.cc +++ b/src/musrScintSD.cc @@ -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="<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 ="<